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': 20,
          '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': 9,
          'figure.subplot.top': 0.94,'figure.subplot.left': 0.105,'figure.subplot.bottom': 0.095,'figure.subplot.right': 0.978}
#'figure.subplot.bottom': 0.097

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" from file name 
try:
   os.remove(scriptName+'.pdf')
   os.remove(scriptName+'.png')
except:
   pass

##############################################################################################################################
hug = numpy.loadtxt('FPEOS_mixture_Hugoniot.txt',usecols=(1,3,5,7,9,13)); # T,rho/rho0,rho0,rho,P,E
rho0 = hug[0,2]

cx = 1
cy = 4
l  = 5
plot(hug[:,cx],hug[:,cy],'-', linewidth=l, markersize=0,color='red',dashes=(16/l,1/l,1/l,1/l,1/l,1/l),mec='r',mew=2, mfc='white', label='Hugoniot curve',zorder=10);

##############################################################################################################################
d   = numpy.loadtxt('FPEOS_mixture_isotherms.txt'      ,usecols=(1,3,7,9,11)); # index,rho,T,P,E

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

for i in range(iMin,iMax+1):
    ii = (d[:,0]==i)
    if (np.where(ii==True)[0].shape[0]>0):
        iii = np.where(ii==True)[0][0]

        label = ''
        if (i==iMin): label='Isotherm'
        plot(d[ii,1]/rho0,d[ii,3],'sb-', linewidth=1, markersize=0,color='b',mec='b',mew=1, mfc='yellow',label=label);

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

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

a.xaxis.set_minor_locator(MultipleLocator(0.2));
a.set_xlim(       floor(         min(hug[:,cx]/2))  ,     ceil(     1.2*max(hug[:,cx]))  )
a.set_ylim( 10.0**floor(np.log10(min(hug[:,cy]))) , 10.0**ceil(np.log10(max(hug[:,cy]))) )

a.set_yscale('log')
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')

legend(loc='upper left',frameon=True,handlelength=2.0,ncol=1,handletextpad=0.4,numpoints=1)

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