from matplotlib.pylab import *
import numpy
import sys
import os

fig_size = [800/72.27 ,700/72.27]
params = {'backend': 'ps', 'axes.labelsize': 25, 'font.size': 25, 'legend.fontsize': 14,
          'xtick.labelsize': 25, 'ytick.labelsize': 25, 
          'xtick.major.size': 15,'ytick.major.size': 15,
          'xtick.minor.size': 8,'ytick.minor.size': 8,
          'text.usetex': False, 'figure.figsize': fig_size, 'xtick.major.pad': 11,
          'figure.subplot.top': 0.94,'figure.subplot.left': 0.105,'figure.subplot.bottom': 0.095,'figure.subplot.right': 0.978}

rcParams.update(params)
a=subplot(111)

title_ = "First-Principles Equation of State Database";
if (len(sys.argv)>1):
   title_ = str(sys.argv[1])
title(title_)

scriptName = sys.argv[0][:-3]; # [:3] removes ".py" for file name 
try:
   os.remove(scriptName+'.pdf')
   os.remove(scriptName+'.png')
except:
   pass

d  = numpy.loadtxt('FPEOS_isochores_ideal_Debye.txt',usecols=(1,3,7,25,27,29,31)); # index,rho,T,1-P/PIdeal,1-E/EIdeal,1-PDebye/PIdeal,1-EDebye/EIdeal
dp = numpy.loadtxt('FPEOS_isochor_points.txt',      usecols=(1,3,7,25,27,29,31)); # index,rho,T,1-P/PIdeal,1-E/EIdeal,1-PDebye/PIdeal,1-EDebye/EIdeal

iMin = min(d[:,0])
iMax = max(d[:,0])
#print(iMin,iMax)
iMin = min(dp[:,0])
iMax = max(dp[:,0])
#print(iMin,iMax)

i = (d[:,0]==iMax)
ii = np.where(i==True)[0][0]
D = (4,0.5,0.5,0.5)
rho_str  =  "{:6.3f}".format(  d[ii,1] )   # 3 decimals for rho
plot(d[i,2],d[i,4],'r-',  linewidth=4.0, color='darkgreen', dashes=D,                                   zorder=10);
plot(d[i,2],d[i,6],'k--', linewidth=1.5, dashes=D, label=r'$\rho$[g/cc]='+rho_str+' Debye plasma model',zorder=0 );

i = (dp[:,0]==iMax)
ii = np.where(i==True)[0][0]
rho_str  =  "{:6.3f}".format(  dp[ii,1] )   # 3 decimals for rho
plot(dp[i,2],dp[i,4],       'D', markersize=12, mec='darkgreen', mfc='lightgreen',                                                                                      zorder=20);
plot(dp[i,2],dp[i,4]+1000.0,'D', markersize=12, mec='darkgreen', mfc='lightgreen', linewidth=4.0, color='darkgreen', dashes=D, label=r'$\rho$[g/cc]='+rho_str+' FPEOS', zorder=20);


i = (d[:,0]==(iMax+iMin)//2)
ii = np.where(i==True)[0][0]
D = (8,1,1,1,1,1,1,1,1,1)
rho_str  =  "{:6.3f}".format(  d[ii,1] )   # 3 decimals for rho
plot(d[i,2],d[i,4],'r-',  linewidth=4.0, color='b', dashes=D,                                           zorder=10);
plot(d[i,2],d[i,6],'k--', linewidth=1.5, dashes=D, label=r'$\rho$[g/cc]='+rho_str+' Debye plasma model',zorder=0 );

i = (dp[:,0]==(iMax+iMin)//2)
ii = np.where(i==True)[0][0]
rho_str  =  "{:6.3f}".format(  dp[ii,1] )   # 3 decimals for rho
plot(dp[i,2],dp[i,4],       'o', markersize=16, mec='b', mfc='lightblue',                                                                              zorder=20);
plot(dp[i,2],dp[i,4]+1000.0,'o', markersize=16, mec='b', mfc='lightblue', linewidth=4.0, color='b', dashes=D, label=r'$\rho$[g/cc]='+rho_str+' FPEOS', zorder=20);


i = (d[:,0]==(iMax+iMin)//4)
ii = np.where(i==True)[0][0]
D = ()
rho_str  =  "{:6.3f}".format(  d[ii,1] )   # 3 decimals for rho
plot(d[i,2],d[i,4],'r-',  linewidth=4.0, color='r', dashes=D,                                           zorder=10);
plot(d[i,2],d[i,6],'k--', linewidth=1.5, dashes=D, label=r'$\rho$[g/cc]='+rho_str+' Debye plasma model',zorder=0 );

i = (dp[:,0]==(iMax+iMin)//4)
ii = np.where(i==True)[0][0]
rho_str  =  "{:6.3f}".format(  dp[ii,1] )   # 3 decimals for rho
plot(dp[i,2],dp[i,4],       's', markersize=16, mec='r', mfc='w',                                                                                    zorder=20);
plot(dp[i,2],dp[i,4]+1000.0,'s', markersize=16, mec='r', mfc='w', linewidth=4.0, color='r', dashes=D, label=r'$\rho$[g/cc]='+rho_str+' FPEOS', zorder=20);

#############################################################################

a.set_ylabel(r'(E-E$_0$)/E$_0$')
#xlabel(r'Density (g$\,$cm$^{-3}$)')
#xlabel(r'Compression ratio $\rho / \rho_0$')
a.set_xlabel(r'Temperature (K)')
#ylabel(r'Pressure (GPa)')

a.set_xscale('log')
a.yaxis.set_minor_locator(MultipleLocator(0.01))

a.xaxis.set_ticks_position('both')
a.get_xaxis().set_tick_params(which='both', direction='in')
a.yaxis.set_ticks_position('both')
a.get_yaxis().set_tick_params(which='both', direction='in')

a.legend(loc='lower right',frameon=False,handlelength=4.5,ncol=1,handletextpad=0.4,numpoints=1)

a.set_xlim(1e6,1e9)
plot([1e6,1e9],[0,0],'--', color='grey',linewidth=1);#print DFT
a.set_ylim(-0.25,0.01)

savefig(scriptName+'.pdf', bbox_inches='tight')
savefig(scriptName+'.png', bbox_inches='tight', dpi=200)