In [ ]:
import matplotlib.pyplot as plt
%matplotlib notebook
import numpy as np
import pandas as pd
import pickle
import xmeos
from xmeos import models
from xmeos import datamod
CONSTS = models.CONSTS
import copy
from mpltools import annotation
from collections import OrderedDict
In [ ]:
In [ ]:
analysis_file = 'data/analysis.pkl'
with open(analysis_file, 'rb') as f:
analysis = pickle.load(f)
In [ ]:
In [ ]:
In [ ]:
datamodel = analysis['datamodel']
data = datamodel['data']
eos_mod = datamodel['eos_mod']
In [ ]:
T0 = eos_mod.get_refstate()['T0']
In [ ]:
In [ ]:
In [ ]:
Nsamp = 300
plt.rc('text', usetex=True)
lbl_siz = 16-2
par_lbl_siz = lbl_siz-2
par_col = np.array([43,71,20])/255.
par_col = np.array([78,128,0])/255.
f, ax = plt.subplots(2, 2, sharex='col')
f.set_size_inches(8,6)
# f.subplots_adjust(hspace=0.05,wspace=0.25,left=.1,right=.97,top=.95)
f.subplots_adjust(hspace=0.03,wspace=0.18,left=.05,right=.97,top=.97,bottom=.05)
# plt.draw()
# Set labels
ax[0,0].set_ylabel(r'$P$')
ax[1,1].set_ylabel(r'$E$')
ax[0,1].set_ylabel(r'$C_{P/V}$')
ax[1,0].set_ylabel(r'$T$')
ax[1,0].set_xlabel(r'$V$')
ax[1,1].set_xlabel(r'$T$')
ax[0,0].set_yticklabels([])
ax[1,0].set_yticklabels([])
ax[0,1].set_yticklabels([])
ax[1,1].set_yticklabels([])
ax[1,0].set_xticklabels([])
ax[1,1].set_xticklabels([])
Vbnds=[0.4,1.2]
V0 = eos_mod.get_params()['V0']
Vmod = np.linspace(Vbnds[0],Vbnds[1],Nsamp)*V0
Vdat = V0*np.array([0.5,.6,.7,.8,.9,1.0])
def make_PV_diagram(Vmod, Vdat, eos_mod, ax, delT=4000):
V0 = eos_mod.get_params()['V0']
T0 = eos_mod.get_refstate()['T0']
Pmod = eos_mod.press(Vmod, T0)
Pmod_hi = eos_mod.press(Vmod, T0+delT)
Pdat = eos_mod.press(Vdat, T0)
Pdat_hi = eos_mod.press(Vdat, T0+delT)
ax.plot(Vmod/V0,Pmod,'k-', Vmod/V0,Pmod_hi,'r-')
ax.plot(Vdat/V0,Pdat,'ko', mew=2,markeredgecolor='k')
ax.plot(Vdat/V0,Pdat_hi,'ro', mew=2,markeredgecolor='r')
ax.set_ylim(-2,190)
ax.set_xlim(0.49,1.11)
ax.text(.65,17-12,r'$P(V,T_0)$',horizontalalignment='right',fontsize=lbl_siz)
ax.text(.65,85,r"1. Isothermal Compression "+"\n"+r"$\;\;\;\;\; \{V_0,K_0,K'_0\}$",
horizontalalignment='left', fontsize=par_lbl_siz,color=par_col,
bbox=dict(lw=2,boxstyle="round", fc=[.95,1,0.95], ec=par_col))
ax.text(.65,60,r'$P(V,T)$',horizontalalignment='left',color='r',fontsize=lbl_siz)
ax.text(.57,170+5,r"$\rm{(I)\,Isothermal\,Compression:}$"+
r" $P(V,T_0)$"+"\n"+r"$\rm{(II)\,Thermal\,Pressure\,Deriv:}$"+"\n"+
r"$\;\;\;\;\;\;\;\;\left.\frac{dP}{dT}\right|_V = \left.\frac{dS}{dV}\right|_T$",horizontalalignment='left',
verticalalignment='top', fontsize=par_lbl_siz,
bbox=dict(lw=2,boxstyle="square", fc=".9", ec='k'))
ax.text(0.95,30,r"$+\Delta T$",verticalalignment='top',
horizontalalignment='left',fontsize=par_lbl_siz+6,color='r')
for (iV,iP0,iPT) in zip(Vdat/V0,Pdat,Pdat_hi):
idy = iPT-iP0
ax.annotate("", xy=(iV, iPT), xycoords='data',
xytext=(iV, iP0), textcoords='data',
arrowprops=dict( facecolor='r',shrink=0.0,
width=3.,headwidth=6.,headlength=3.0,edgecolor='k'))
pass
def make_TV_diagram(Vmod, Vdat, eos_mod, ax, delT=4000):
V0 = eos_mod.get_params()['V0']
T0 = eos_mod.get_refstate()['T0']
ax.plot(Vmod/V0, eos_mod.ref_temp_adiabat(Vmod)/T0,'k-')
ax.plot(Vdat/V0, eos_mod.ref_temp_adiabat( Vdat)/T0,'ko')
ax.set_xlim(0.49,1.11)
ax.set_ylim(2400/T0,4700/T0)
ax.text(.52,1.27-.22-.05,r'$T_{0S}(V)$',horizontalalignment='left',fontsize=lbl_siz)
# ax_a[1,0].text(.52,1.18-.22,r"$p_2= \{\gamma_0,\gamma_0'\}$",color=par_col,
# horizontalalignment='left',fontsize=par_lbl_siz)
ax.text(.72,1.35,r"2. Adiabatic Reference: "+"\n"+r"$\;\;\;\;\; \{\gamma_0,\gamma_0'\}$",color=par_col,
horizontalalignment='left',fontsize=par_lbl_siz,
bbox=dict(lw=2,boxstyle="round", fc=[.95,1,0.95], ec=par_col))
ax.text(0.58,1.75,r"$\rm{(IIIb)\,Adiabatic\,Temperature}$"+
r"$\rm{\,Deriv:}$"+"\n"+ r" $\;\;\;\;\;\;\left.\frac{dT}{dV}\right|_{S} = -\left.\frac{dS}{dV}\right|_{T} \; / \; \left.\frac{dS}{dT}\right|_{V} $",horizontalalignment='left',
verticalalignment='top', fontsize=par_lbl_siz,
bbox=dict(lw=2,boxstyle="square", fc=".9", ec='k'))
annotation.slope_marker((.721, 1.21), -1.15, labels=(r'$dT$',r'$dV$'),
ax=ax)
pass
make_PV_diagram(Vmod, Vdat, eos_mod, ax[0,0], delT=4000)
make_TV_diagram(Vmod, Vdat, eos_mod, ax[1,0], delT=4000)
In [ ]:
In [ ]:
In [ ]:
def plot_rtpress_model_fig():
plt.rc('text', usetex=True)
lbl_siz = 16-2
par_lbl_siz = lbl_siz-2
par_col = np.array([43,71,20])/255.
par_col = np.array([78,128,0])/255.
f, ax_a = plt.subplots(2, 2, sharex='col')
f.set_size_inches(8,6)
# f.subplots_adjust(hspace=0.05,wspace=0.25,left=.1,right=.97,top=.95)
f.subplots_adjust(hspace=0.03,wspace=0.18,left=.05,right=.97,top=.97,bottom=.05)
plt.draw()
# ax_P = plt.subplot(221)
# ax_E = plt.subplot(222)
# ax_T = plt.subplot(223)
# ax_Cpv = plt.subplot(224)
# Fix ticks
# ax_P.set_xticklabels([])
# ax_E.set_xticklabels([])
# Set labels
ax_a[0,0].set_ylabel(r'$P$')
ax_a[1,1].set_ylabel(r'$E$')
ax_a[0,1].set_ylabel(r'$C_{P/V}$')
ax_a[1,0].set_ylabel(r'$T$')
ax_a[1,0].set_xlabel(r'$V$')
ax_a[1,1].set_xlabel(r'$T$')
ax_a[0,0].set_yticklabels([])
ax_a[1,0].set_yticklabels([])
ax_a[0,1].set_yticklabels([])
ax_a[1,1].set_yticklabels([])
ax_a[1,0].set_xticklabels([])
ax_a[1,1].set_xticklabels([])
plt.draw()
Vbnds=[0.4,1.2]
eos_d = datamod_d['eos_d']
full_mod = eos_d['modtype_d']['FullMod']
gamma_mod = eos_d['modtype_d']['GammaMod']
V0 = eos_d['param_d']['V0']
Vgrid_a = np.linspace(Vbnds[0],Vbnds[1],Nsamp)*V0
Tgrid_a = np.array([2500,7000])
E_grid_a, P_grid_a, dPdT_grid_a = eval_model( Vgrid_a, Tgrid_a, eos_d )
Vdat_a = V0*np.array([0.5,.6,.7,.8,.9,1.0])
Pdat_a = full_mod.press(Vdat_a,Tgrid_a[0],eos_d)
Pdat_hi_a = full_mod.press(Vdat_a,Tgrid_a[1],eos_d)
# P-V plot
ax_a[0,0].plot(Vgrid_a/V0,P_grid_a[0],'k-',
Vgrid_a/V0,P_grid_a[1],'r-')
ax_a[0,0].plot(Vdat_a/V0,Pdat_a,'ko',mew=2,markeredgecolor='k')
ax_a[0,0].plot(Vdat_a/V0,Pdat_hi_a,'ro',mew=2,markeredgecolor='r')
ax_a[0,0].text(.65,17-12,r'$P(V,T_0)$',horizontalalignment='right',fontsize=lbl_siz)
# ax_a[0,0].text(.705,2,r"$p_1= \{V_0,K_0,K'_0\}$",horizontalalignment='right',
# fontsize=par_lbl_siz,color=par_col)
ax_a[0,0].text(.65,85,r"1. Isothermal Compression "+"\n"+r"$\;\;\;\;\; \{V_0,K_0,K'_0\}$",
horizontalalignment='left', fontsize=par_lbl_siz,color=par_col,
bbox=dict(lw=2,boxstyle="round", fc=[.95,1,0.95], ec=par_col))
ax_a[0,0].text(.65,60,r'$P(V,T)$',horizontalalignment='left',color='r',fontsize=lbl_siz)
ax_a[0,0].text(.57,170+5,r"$\rm{(I)\,Isothermal\,Compression:}$"+
r" $P(V,T_0)$"+"\n"+r"$\rm{(II)\,Thermal\,Pressure\,Deriv:}$"+"\n"+
r"$\;\;\;\;\;\;\;\;\left.\frac{dP}{dT}\right|_V = \left.\frac{dS}{dV}\right|_T$",horizontalalignment='left',
verticalalignment='top', fontsize=par_lbl_siz,
bbox=dict(lw=2,boxstyle="square", fc=".9", ec='k'))
# ax_a[0,0].text(.53,150+5,r"$\rm{(II)\,Thermal\,Pressure\,Deriv:}$"+"\n"+
# r"$\;\;\;\;\;\;\left.\frac{dP}{dT}\right|_V = \left.\frac{dS}{dV}\right|_T$",horizontalalignment='left',
# verticalalignment='top', fontsize=par_lbl_siz)
# ax_a[0,0].text(.6,140,r"$\rm{(II)\,Thermal\,Pressure}$"+"\n"+
# r"$\;\;P(V,T_0)$",horizontalalignment='left',
# verticalalignment='top', fontsize=par_lbl_siz)
ax_a[0,0].text(0.95,30,r"$+\Delta T$",verticalalignment='top',
horizontalalignment='left',fontsize=par_lbl_siz+6,color='r')
ax_a[0,0].set_xlim(0.49,1.11)
ax_a[0,0].set_ylim(-2,190)
for (iV,iP0,iPT) in zip(Vdat_a/V0,Pdat_a,Pdat_hi_a):
idy = iPT-iP0
ax_a[0,0].annotate("", xy=(iV, iPT), xycoords='data',
xytext=(iV, iP0), textcoords='data',
arrowprops=dict( facecolor='r',shrink=0.0,
width=3.,headwidth=6.,headlength=3.0,edgecolor='k'))
plt.draw()
T0 = eos_d['param_d']['T0']
# T-V plot
ax_a[1,0].plot(Vgrid_a/V0, gamma_mod.temp( Vgrid_a,T0,eos_d)/T0,'k-')
ax_a[1,0].plot(Vdat_a/V0, gamma_mod.temp( Vdat_a,T0,eos_d)/T0,'ko')
ax_a[1,0].set_xlim(0.49,1.11)
ax_a[1,0].set_ylim(2400/T0,4700/T0)
ax_a[1,0].text(.52,1.27-.22-.05,r'$T_{0S}(V)$',horizontalalignment='left',fontsize=lbl_siz)
# ax_a[1,0].text(.52,1.18-.22,r"$p_2= \{\gamma_0,\gamma_0'\}$",color=par_col,
# horizontalalignment='left',fontsize=par_lbl_siz)
ax_a[1,0].text(.72,1.35,r"2. Adiabatic Reference: "+"\n"+r"$\;\;\;\;\; \{\gamma_0,\gamma_0'\}$",color=par_col,
horizontalalignment='left',fontsize=par_lbl_siz,
bbox=dict(lw=2,boxstyle="round", fc=[.95,1,0.95], ec=par_col))
ax_a[1,0].text(0.58,1.75,r"$\rm{(IIIb)\,Adiabatic\,Temperature}$"+
r"$\rm{\,Deriv:}$"+"\n"+ r" $\;\;\;\;\;\;\left.\frac{dT}{dV}\right|_{S} = -\left.\frac{dS}{dV}\right|_{T} \; / \; \left.\frac{dS}{dT}\right|_{V} $",horizontalalignment='left',
verticalalignment='top', fontsize=par_lbl_siz,
bbox=dict(lw=2,boxstyle="square", fc=".9", ec='k'))
annotation.slope_marker((.721, 1.21), -1.15, labels=(r'$dT$',r'$dV$'),
ax=ax_a[1,0])
# E-T plot
E0 = eos_d['param_d']['E0']
Tdat_a = T0*np.array([.6,1.0,1.4,1.8,2.2,2.6,3.0])
Edat_a = full_mod.energy(V0,Tdat_a,eos_d)
Tgrid_a = T0*np.linspace(.5,3.2,1001)
Egrid_a = full_mod.energy(V0,Tgrid_a,eos_d)
ax_a[1,1].plot(Tdat_a/T0, Edat_a-E0,'ko')
ax_a[1,1].plot(Tgrid_a/T0, Egrid_a-E0,'k-')
ax_a[1,1].set_xlim(.9,3.1)
# ax_a[1,1].set_ylim(-.3,1.7)
ax_a[1,1].set_ylim(-.4,2.0)
ax_a[1,1].text(1.5,0.4-.65,r'$E(V_0,T)$',horizontalalignment='right',fontsize=lbl_siz)
ax_a[1,1].text(1.2,1.5+.3,r"$\rm{(IIIc)\,Thermal\,Energy\,Deriv:}$"+
"\n"+r"$\;\;\;\;\;\;\left.\frac{dE}{dT}\right|_V = \left. T \, \frac{dS}{dT}\right|_V$",horizontalalignment='left',
verticalalignment='top', fontsize=par_lbl_siz,
bbox=dict(lw=2,boxstyle="square", fc=".9", ec='k'))
annotation.slope_marker((1.5, 0.35), 0.85, labels=(r'$dE$',r'$dT$'),
ax=ax_a[1,1])
# Cv-T plot
Cv_dat_a = full_mod.heat_capacity(V0,Tdat_a,eos_d)/eos_d['const_d']['kboltz']
Cv_dat_hiP_a = full_mod.heat_capacity(0.5*V0,Tdat_a,eos_d)/eos_d['const_d']['kboltz']
Cv_grid_a = full_mod.heat_capacity(V0,Tgrid_a,eos_d)/eos_d['const_d']['kboltz']
Cv_grid_hiP_a = full_mod.heat_capacity(0.5*V0,Tgrid_a,eos_d)/eos_d['const_d']['kboltz']
ax_a[0,1].plot(Tgrid_a/T0, Cv_grid_a,'k-', Tdat_a/T0,Cv_dat_a,'ko')
ax_a[0,1].plot(Tgrid_a/T0, Cv_grid_hiP_a,'b-')
ax_a[0,1].set_xlim(.9,3.1)
#ax_a[0,1].set_ylim(2.9,4.4)
ax_a[0,1].set_ylim(3.3,5.1)
# ax_a[0,1].text(0.75,4.15,r'$C_{P/V}(T)$',horizontalalignment='left',fontsize=lbl_siz)
ax_a[0,1].text(1.5,3.35-.26+.3,r'$C_{P/V}(T)$',horizontalalignment='right',fontsize=lbl_siz)
# ax_a[0,1].text(1.05,3.2-.23,r"$p_3= \{b_0,...,b_n\}$",color=par_col,
# horizontalalignment='left',fontsize=par_lbl_siz)
ax_a[0,1].text(1.7,3.8+.6,r"3. RT Thermal Coeffs: "+"\n"+r"$\;\;\;\;\; \{b_0,...,b_n\}$",color=par_col,
horizontalalignment='left',fontsize=par_lbl_siz,
bbox=dict(lw=2,boxstyle="round", fc=[.95,1,0.95], ec=par_col))
# ax_a[0,1].text(2.5,3.03,"low press",verticalalignment='top',
# horizontalalignment='left',fontsize=par_lbl_siz,color='k')
# ax_a[0,1].text(2.5,3.47,"high press",verticalalignment='top',
# horizontalalignment='left',fontsize=par_lbl_siz,color='b')
ax_a[0,1].text(2.6,3.45+.45,r"$+\Delta P$",verticalalignment='top',
horizontalalignment='left',fontsize=par_lbl_siz+6,color='b')
ax_a[0,1].text(1.3,4.3+.6,r"$\rm{(IIIa)\,Heat\,Capacity:}$"+
r" $C_{P/V} = T\,\frac{dS}{dT}$",horizontalalignment='left',
verticalalignment='top', fontsize=par_lbl_siz,
bbox=dict(lw=2,boxstyle="square", fc=".9", ec='k'))
for (iT,iCv0,iCvT) in zip(Tdat_a/T0,Cv_dat_a,Cv_dat_hiP_a):
ax_a[0,1].annotate("", xy=(iT, iCvT), xycoords='data',
xytext=(iT, iCv0), textcoords='data',
arrowprops=dict( facecolor='b',shrink=0.0,
width=3.,headwidth=6.,headlength=3.0,edgecolor='k'))
plt.draw()
plt.draw()
return ax_a
In [ ]:
# Make RTpress model construction plot
ax_a = plot_rtpress_model_fig()
# plt.savefig('figs/MgSiO3-rtpress-modelfit-Spera2011.eps')
In [ ]: