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_Hugoniot.txt',usecols=(23,25)); # up,us

"""
n2 = hug.shape[0]
n1 = n2 - 1
while (hug[n1,1]<100.0 and n1>10):
   n1 = n1 - 1;
"""

n = min(10,hug.shape[0])
while (n<hug.shape[0] and hug[n,1]<50.0):
   n += 1

l = 5
plot(hug[:n,0],hug[:n,1],'-', linewidth=5, markersize=0,color='red',mec='r',mew=2, mfc='white', label='Hugoniot curve');

##############################################################################################################################
hug = numpy.loadtxt('FPEOS_Hugoniot_lower_initial_density.txt',usecols=(23,25)); # up,us

n = min(10,hug.shape[0])
while (n<hug.shape[0] and hug[n,1]<50.0):
   n += 1

l = 2
plot(hug[0:n,0],hug[0:n,1],'-', linewidth=1, markersize=0,color='darkgreen',dashes=(2/l,1/l),mec='r',mew=2, mfc='white', label='Hugoniot with\n'+r'20% lower $\rho_0$');

##############################################################################################################################
hug = numpy.loadtxt('FPEOS_Hugoniot_higher_initial_density.txt',usecols=(23,25)); # up,us

"""
n2 = hug.shape[0]
n1 = n2 - 1
while (hug[n1,1]<100.0 and n1>10):
   n1 = n1 - 1;
"""
n = min(10,hug.shape[0])
while (n<hug.shape[0] and hug[n,1]<50.0):
   n += 1

l = 2
plot(hug[0:n,0],hug[0:n,1],'-', linewidth=1, markersize=0,color='navy',dashes=(16/l,1/l,1/l,1/l,1/l,1/l),mec='r',mew=2, mfc='white', label='Hugoniot with\n'+r'20% higher $\rho_0$');

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

#xlabel(r'Pressure (GPa)')
xlabel(r'Particle velocity u$_p$ (km/s)')
ylabel(r'Shock velocity u$_s$ (km/s)')

# a.set_xscale('log') # not clear whether a log scale for density is better
# 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='lower right',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)