text: """Classes and functions for accessing and manipulating tabulated EOS data."""
  ### Module for accessing and manipulating tabulated EOS data
  ### STS 09/2019
  ###
  ##
  import numpy as np
  import pandas as pd
  #
  #
  # calculate the structure for one planet
  # make a class to hold the PREM data
  class isentrope_class:
      """Class to isentrope data extracted from EOS table."""  # this is a documentation string for this class
      def __init__(self): # self is the default name of the object for internal referencing of the variables in the class
          """A function to initialize the class object.""" # this is a documentation string for this function
          self.ND = 0 # number of radius points
          self.density     = []   
          self.pressure    = []
          self.temperature = []
          self.soundspeed  = []
          self.energy  = []
          self.partvel = []
          self.region = [] # Tillotson region flag
          # not going to use all the variables in the file
          self.units = '' # I like to keep a text note in a structure about the units
  #
  class EOShugoniot:
      """Class for Hugoniot array from extEOStable."""	
      def __init__(self):
          self.NH = 0
          self.rho = np.zeros(self.NH)   
          self.T = np.zeros(self.NH)   
          self.P = np.zeros(self.NH)   
          self.U = np.zeros(self.NH)   
          self.S = np.zeros(self.NH)   
          self.up = np.zeros(self.NH)   
          self.us = np.zeros(self.NH)
          self.cs = np.zeros(self.NH)
          self.units = ''
  #
  class EOSvaporcurve:
      """Class for vapor curve from ANEOS."""	
      def __init__(self):
          self.NT = 0
          self.NV = 0
          self.T = np.zeros(self.NT)  
          self.rl = np.zeros(self.NT)  
          self.rv = np.zeros(self.NT)  
          self.Pl = np.zeros(self.NT)  
          self.Pv = np.zeros(self.NT)  
          self.Ul = np.zeros(self.NT)  
          self.Uv = np.zeros(self.NT)  
          self.Sl = np.zeros(self.NT)  
          self.Sv = np.zeros(self.NT)
          self.Gl = np.zeros(self.NT)  
          self.Gv = np.zeros(self.NT)
          self.units = ''
  #
  class EOSmeltcurve:
      """Class for melt curve from ANEOS."""	
      def __init__(self):
          self.NT = 0
          self.NV = 0
          self.T  = np.zeros(self.NT)  
          self.rl = np.zeros(self.NT)  
          self.rs = np.zeros(self.NT)  
          self.Pl = np.zeros(self.NT)  
          self.Ps = np.zeros(self.NT)  
          self.Ul = np.zeros(self.NT)  
          self.Us = np.zeros(self.NT)  
          self.Sl = np.zeros(self.NT)  
          self.Ss = np.zeros(self.NT)
          self.units = ''
  #
  class EOS1barcurve:
      """Class for 1bar curve from the EOS."""	
      def __init__(self):
          self.NT    = 0
          self.S     = np.zeros(self.NT)  
          self.T     = np.zeros(self.NT)  
          self.rho   = np.zeros(self.NT)
          self.Tvap  = 0.
          self.Tmelt = 0.
          self.Sim   = 0.
          self.Scm   = 0.
          self.Siv   = 0.
          self.Scv   = 0.
          self.rhoiv   = 0.
          self.rhocv   = 0.
          self.rhocm   = 0.
          self.rhoim   = 0.
          self.units = ''
  #
  class EOScriticalpoint:
      """Class for critical point state from the EOS."""	
      def __init__(self):
          self.P   = 0
          self.S   = 0  
          self.T   = 0 
          self.rho = 0
          self.U   = 0
          self.units = ''
  #
  class EOStriplepoint:
      """Class for triple point state from the EOS."""	
      def __init__(self):
          self.P   = 0
          self.T   = 0 
          self.Sim   = 0.
          self.Scm   = 0.
          self.Siv   = 0.
          self.Scv   = 0.
          self.rhol  = 0.
          self.units = ''
  #
  class EOSaneoshugoniot:
      """Class for Hugoniot calculated in ANEOS."""	
      def __init__(self):
          self.ND  = 0
          self.NV  = 0
          #self.all = np.zeros((self.ND,self.NV))
          self.rho = 0
          self.T   = 0
          self.P   = 0
          self.U   = 0
          self.S   = 0
          self.us  = 0
          self.up  = 0
          self.units = ''
  #
  class extEOStable:
      """Class for accessing EXTENDED SESAME-STYLE EOS tables output from ANEOS"""
      #     ANEOS KPA FLAG
      #                                TABLE          ANEOS
      #     KPAQQ=STATE INDICATOR      =1, 1p    =1, 1p    (eos without melt)
      #                                =2, 2p lv =2, 2p liquid/solid plus vapor
      #                                          =4, 1p solid  (eos with melt)
      #                                          =5, 2p melt   (eos with melt)
      #                                          =6, 1p liquid (eos with melt)
      #                                =-1 bad value of temperature
      #                                =-2 bad value of density
      #                                =-3 bad value of material number
      #
      def __init__(self):
          self.ND  = 0 # integer; number of density points in grid
          self.NT  = 0 # integer; number of temperature points in grid
          self.rho = np.zeros(self.ND)          # g/cm3, density values
          self.T   = np.zeros(self.NT)          # K, temperature values
          self.P   = np.zeros(self.ND*self.NT)  # GPA, pressure(T,rho)
          self.U   = np.zeros(self.ND*self.NT)  # MJ/kg, sp. internal energy(T,rho)
          self.A   = np.zeros(self.ND*self.NT)  # MJ/kg, Helmholtz free energy(T,rho)
          self.S   = np.zeros(self.ND*self.NT)  # MJ/K/kg, sp. entropy(T,rho)
          self.cs  = np.zeros(self.ND*self.NT)  # cm/s, sound speed(T,rho)
          self.cv  = np.zeros(self.ND*self.NT)  # MJ/K/kg, sp. heat capacity(T,rho)
          self.KPA = np.zeros(self.ND*self.NT)  # integer, ANEOS KPA flag(T,rho)
          self.MDQ = np.zeros(self.ND*self.NT)  # integer, Model Development Quality Flag(T,rho)
          self.units = ''
          self.hug = EOShugoniot()
          self.vc  = EOSvaporcurve()
          self.mc  = EOSmeltcurve()
          self.cp  = EOScriticalpoint()
          self.tp  = EOStriplepoint()
          self.onebar = EOS1barcurve()
          self.anhug = EOSaneoshugoniot()
          # these are variables needed for the sesame header
          self.MATID   = 0.
          self.DATE    = 0.
          self.VERSION = 0.
          self.FMN     = 0.
          self.FMW     = 0.
          self.R0REF   = 0.
          self.K0REF   = 0.
          self.T0REF   = 0.
          self.P0REF   = 0.
          # variables needed for the ANEOS gamma function
          self.gamma0 = 0.
          self.theta0 = 0.
          self.C24    = 0.
          self.C60    = 0.
          self.C61    = 0.
          self.beta   = 0.
          # model name/version string
          self.MODELNAME = ''
  
      def loadstdsesame(self, fname, unitstxt=None):
          """Function for loading STD SESAME-STYLE EOS table output from ANEOS"""
          data = ([])
          if unitstxt is None:
              self.units = 'Units: rho g/cm3, T K, P GPa, U MJ/kg, A MJ/kg, S MJ/K/kg, cs cm/s, cv MJ/K/kg, KPA flag. 2D arrays are (NT,ND).'
          else:
              self.units = unitstxt
          sesamefile = open(fname,"r")  
          sesamedata=sesamefile.readlines()
          sesamefile.close()
          nskip = 6 # skip standard header to get to the content of the 301 table
          # num.density, num. temps
          tmp = sesamedata[nskip][0:16]
          dlen = float(tmp)
          tmp = sesamedata[nskip][16:32]
          tlen = float(tmp)
          if (np.mod((dlen*tlen*3.0+dlen+tlen+2.0),5.0) == 0):
              neos = int((dlen*tlen*3.0+dlen+tlen+2.0)/5.0) 
          else:
              neos = int((dlen*tlen*3.0+dlen+tlen+2.0)/5.0) +1
          #print(dlen,tlen,neos,len(sesamedata))
          data = np.zeros((neos,5),dtype=float)
          for j in range(nskip,neos+nskip):
              tmp3 = sesamedata[j]
              tmp4 = list(tmp3.split())
              if len(tmp4) < 5:
                  lentmp4 = len(tmp4)
                  data[j-nskip,0:lentmp4] = np.asarray(tmp4[0:lentmp4])
              else:
                  data[j-nskip,:] = np.asarray(tmp4)
              #print(j,eosarr[j,:])
          #print(data.shape)
          data=np.resize(data,(neos*5))
          #print(data.shape)
          self.ND  = data[0].astype(int)  # now fill the extEOStable class
          self.NT  = data[1].astype(int)
          self.rho = data[2:2+self.ND]
          self.T   = data[2+self.ND : 2+self.ND+self.NT]
          self.P   = data[2+self.ND+self.NT : 2+self.ND+self.NT+self.ND*self.NT
                              ].reshape(self.NT,self.ND)
          self.U   = data[2+self.ND+self.NT+self.ND*self.NT
                              : 2+self.ND+self.NT+2*self.ND*self.NT
                              ].reshape(self.NT,self.ND)
          self.A   = data[2+self.ND+self.NT+2*self.ND*self.NT
                              : 2+self.ND+self.NT+3*self.ND*self.NT
                              ].reshape(self.NT,self.ND)
          #self.S   = data[2+self.ND+self.NT+3*self.ND*self.NT
          #                    : 2+self.ND+self.NT+4*self.ND*self.NT
          #                    ].reshape(self.NT,self.ND)
          #self.cs  = data[2+self.ND+self.NT+4*self.ND*self.NT
          #                    : 2+self.ND+self.NT+5*self.ND*self.NT
          #                    ].reshape(self.NT,self.ND)
          #self.cv  = data[2+self.ND+self.NT+5*self.ND*self.NT
          #                    : 2+self.ND+self.NT+6*self.ND*self.NT
          #                    ].reshape(self.NT,self.ND)
          #self.KPA = data[2+self.ND+self.NT+6*self.ND*self.NT
          #                    : 2+self.ND+self.NT+7*self.ND*self.NT
          #                    ].reshape(self.NT,self.ND)
  #
      def loadextsesame(self, fname, unitstxt=None):
          """Function for loading EXTENDED SESAME-STYLE EOS table output from ANEOS"""
          data = ([])
          if unitstxt is None:
              self.units = 'Units: rho g/cm3, T K, P GPa, U MJ/kg, A MJ/kg, S MJ/K/kg, cs cm/s, cv MJ/K/kg, KPA flag. 2D arrays are (NT,ND).'
          else:
              self.units = unitstxt
          sesamefile = open(fname,"r")  
          sesamedata=sesamefile.readlines()
          sesamefile.close()
          nskip = 6 # skip standard header to get to the content of the 301 table
          # num.density, num. temps
          tmp = sesamedata[nskip][0:16]
          dlen = float(tmp)
          tmp = sesamedata[nskip][16:32]
          tlen = float(tmp)
          if (np.mod((dlen*tlen*4.0+dlen+tlen+2.0),5.0) == 0):
              neos = int((dlen*tlen*4.0+dlen+tlen+2.0)/5.0)
          else:
              neos = int((dlen*tlen*4.0+dlen+tlen+2.0)/5.0) +1
          #print(dlen,tlen,neos,len(sesamedata))
          data = np.zeros((neos,5),dtype=float)
          for j in range(nskip,neos+nskip):
              tmp3 = sesamedata[j]
              tmp4 = list(tmp3.split())
              if len(tmp4) < 5:
                  lentmp4 = len(tmp4)
                  data[j-nskip,0:lentmp4] = np.asarray(tmp4[0:lentmp4])
              else:
                  data[j-nskip,:] = np.asarray(tmp4)        
              #print(j,eosarr[j,:])
          #print(data.shape)
          data=np.resize(data,(neos*5))
          #print(data.shape)
          self.ND  = data[0].astype(int)  # now fill the extEOStable class
          self.NT  = data[1].astype(int)
          self.rho = data[2:2+self.ND]
          self.T   = data[2+self.ND : 2+self.ND+self.NT]
          #self.P   = data[2+self.ND+self.NT : 2+self.ND+self.NT+self.ND*self.NT
          #                    ].reshape(self.NT,self.ND)
          #self.U   = data[2+self.ND+self.NT+self.ND*self.NT
          #                    : 2+self.ND+self.NT+2*self.ND*self.NT
          #                    ].reshape(self.NT,self.ND)
          #self.A   = data[2+self.ND+self.NT+2*self.ND*self.NT
          #                    : 2+self.ND+self.NT+3*self.ND*self.NT
          #                    ].reshape(self.NT,self.ND)
          self.S   = data[2+self.ND+self.NT+0*self.ND*self.NT
                              : 2+self.ND+self.NT+1*self.ND*self.NT
                              ].reshape(self.NT,self.ND)
          self.cs  = data[2+self.ND+self.NT+1*self.ND*self.NT
                              : 2+self.ND+self.NT+2*self.ND*self.NT
                              ].reshape(self.NT,self.ND)
          self.cv  = data[2+self.ND+self.NT+2*self.ND*self.NT
                              : 2+self.ND+self.NT+3*self.ND*self.NT
                              ].reshape(self.NT,self.ND)
          self.KPA = data[2+self.ND+self.NT+3*self.ND*self.NT
                              : 2+self.ND+self.NT+4*self.ND*self.NT
                              ].reshape(self.NT,self.ND)
  #
      def view(self, q='P', Tlow=None, Thigh=None, rholow=None, rhohigh=None):
          """Function for printing values from EXTENDED SESAME-STYLE EOS table."""
          if Tlow is None:
              Tlow = self.T.min()
          if Thigh is None:
              Thigh = self.T.max()
          if rholow is None:
              rholow = self.rho.min()
          if rhohigh is None:
              rhohigh = self.rho.max()
          print(self.units)
          if q == 'P':
              print('T:', self.T[np.logical_and(self.T >= Tlow,self.T <= Thigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho <= rhohigh)
                                    ])
              print('P:', (self.P[np.logical_and(self.T >= Tlow,
                                                  self.T<=Thigh)
                                 ])[:, np.logical_and(self.rho >= rholow,
                                                       self.rho <= rhohigh)
                                   ])
          if q == 'U':
              print('T:', self.T[np.logical_and(self.T >= Tlow,self.T <= Thigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho <= rhohigh)
                                    ])
              print('U:', (self.U[np.logical_and(self.T >= Tlow,
                                                  self.T <= Thigh)
                                 ])[:, np.logical_and(self.rho >= rholow,
                                                       self.rho <= rhohigh)
                                   ])
          if q == 'A':
              print('T:', self.T[np.logical_and(self.T >= Tlow,self.T <= Thigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho <= rhohigh)
                                    ])
              print('A:', (self.A[np.logical_and(self.T >= Tlow,
                                                  self.T <= Thigh)
                                 ])[:, np.logical_and(self.rho >= rholow,
                                                       self.rho <= rhohigh)
                                   ])
          if q == 'S':
              print('T:', self.T[np.logical_and(self.T >= Tlow,self.T <= Thigh)])
              print('rho:', self.rho[np.logical_and(self.rho
                                     >= rholow,self.rho<=rhohigh)
                                    ])
              print('S:', (self.S[np.logical_and(self.T >= Tlow,self.T <= Thigh)
                                 ])[:, np.logical_and(self.rho >= rholow,
                                                       self.rho <= rhohigh)
                                   ])      
          if q == 'cs':
              print('T:', self.T[np.logical_and(self.T >= Tlow,self.T <= Thigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho<=rhohigh)])
              print('cs:', (self.cs[np.logical_and(self.T >= Tlow,
                                                    self.T <= Thigh)
                                   ])[:, np.logical_and(self.rho >= rholow,
                                                         self.rho <= rhohigh)
                                     ])
          if q == 'cv':
              print('T:', self.T[np.logical_and(self.T >= Tlow,self.T <= Thigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho<=rhohigh)])
              print('cv:', (self.cv[np.logical_and(self.T >= Tlow,
                                                    self.T <= Thigh)
                                   ])[:, np.logical_and(self.rho >= rholow,
                                                         self.rho <= rhohigh)
                                     ])
          if q == 'KPA':
              print('T:', self.T[np.logical_and(self.T >= Tlow,self.T <= Thigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho<=rhohigh)])
              print('KPA:', (self.KPA[np.logical_and(self.T >= Tlow,
                                                    self.T <= Thigh)
                                   ])[:, np.logical_and(self.rho >= rholow,
                                                         self.rho <= rhohigh)
                                     ])
  
      def calchugoniot(self, r0=None, t0=None, pmax=None, writefilename=None):
          """Function for calculating a Hugoniot from EXTENDED SESAME-STYLE EOS table."""
          if r0 is None:
              return 'Must provide r0 and t0.'
          if t0 is None:
              return 'Must provide r0 and t0.'
          if pmax is None:
              pmax=1.E4 # GPa
          self.hug.rho = []
          self.hug.P = []
          self.hug.T = []
          self.hug.U = []
          self.hug.S = []
          self.hug.up = []
          self.hug.us = []
          self.hug.cs = []
    
          it0 = int(np.round(np.interp(t0,self.T,np.arange(self.NT)))) # uses nearest value if t0 not in array
          ir0 = int(np.round(np.interp(r0,self.rho,np.arange(self.ND)))) # uses nearest value if r0 not in the array
          p0  = self.P[it0,ir0] # GPa
          #print(self.P[it0,ir0])
          e0  = self.U[it0,ir0]#np.interp(p0,self.P[it0,:],self.U[it0,:])
          s0  = self.S[it0,ir0]#np.interp(p0,self.P[it0,:],self.S[it0,:])
          up0 = 0. # no initial particle velocity
          us0 = self.cs[it0,ir0]/1.e5 # cm/s->km/s use sound velocity for initial
          cs0 = self.cs[it0,ir0]/1.e5 # cm/s->km/s use sound velocity for initial
          #print(ir0,it0,r0,t0,p0,e0,up0,us0)
          self.hug.rho = np.append(self.hug.rho, self.rho[ir0])
          self.hug.P = np.append(self.hug.P, p0)
          self.hug.T = np.append(self.hug.T, self.T[it0])
          self.hug.U = np.append(self.hug.U, e0)
          self.hug.S = np.append(self.hug.S, s0)
          self.hug.up = np.append(self.hug.up, up0)
          self.hug.us = np.append(self.hug.us, us0)
          self.hug.cs = np.append(self.hug.cs, cs0)
          
          #for iir in range(ir0+1,self.ND):
          iir=ir0+1
          pnew=p0
          while pnew<pmax:
              ediff =0.5*(self.P[it0::,iir]+p0)*(1./r0-1./self.rho[iir])+e0 -(self.U[it0::,iir])  # MJ/kg
              #if not(pd.Series(ediff).is_monotonic_decreasing and pd.Series(ediff).is_unique):
              #    print('bad ediff:', self.T[it0], self.rho[iir], pnew, tnew)
              # np.interp wants x values increasing
              pnew = np.interp(0.,np.flip(ediff),np.flip(self.P[it0::,iir]))
              tnew = np.interp(0.,np.flip(ediff),np.flip(self.T[it0::]))
              enew = np.interp(0.,np.flip(ediff),np.flip(self.U[it0::,iir]))
              snew = np.interp(0.,np.flip(ediff),np.flip(self.S[it0::,iir]))
              upnew = np.sqrt((pnew-p0)*(1./r0-1./self.rho[iir]))
              usnew = (1./r0)*np.sqrt((pnew-p0)/(1./r0-1./self.rho[iir]))
              csnew = np.interp(0.,np.flip(ediff),np.flip(self.cs[it0::,iir]))/1.E5 # km/s
              #print(self.rho[iir],tnew,pnew,enew,upnew,usnew)
              self.hug.rho = np.append(self.hug.rho, self.rho[iir])
              self.hug.P = np.append(self.hug.P, pnew)
              self.hug.T = np.append(self.hug.T, tnew)
              self.hug.U = np.append(self.hug.U, enew)
              self.hug.S = np.append(self.hug.S, snew)
              self.hug.up = np.append(self.hug.up, upnew)
              self.hug.us = np.append(self.hug.us, usnew)
              self.hug.cs = np.append(self.hug.cs, csnew) # km/s
              iir += 1
          self.hug.NH=len(self.hug.P)
          self.hug.units='units: T K, rho g/cm3, P GPa, U MJ/kg, S MJ/K/kg, Up km/s, Us km/s, cs km/s'
  
      def calcporoushugoniot(self, r0=None, t0=None, pmax=None, writefilename=None, r00=None):
          """Function for calculating a Hugoniot from EXTENDED SESAME-STYLE EOS table."""
          if r0 is None:
              return 'Must provide r0 and t0.'
          if t0 is None:
              return 'Must provide r0 and t0.'
          if pmax is None:
              pmax=1.E4 # GPa
          self.hug.rho = []
          self.hug.P = []
          self.hug.T = []
          self.hug.U = []
          self.hug.S = []
          self.hug.up = []
          self.hug.us = []
          self.hug.cs = []
    
          it0 = int(np.round(np.interp(t0,self.T,np.arange(self.NT)))) # uses nearest value if t0 not in array
          ir0 = int(np.round(np.interp(r0,self.rho,np.arange(self.ND)))) # uses nearest value if r0 not in the array
          ir00 = int(np.round(np.interp(r00,self.rho,np.arange(self.ND)))) # uses nearest value if r0 not in the array
          p0  = self.P[it0,ir0] # GPa
          #print(self.P[it0,ir0])
          e0  = self.U[it0,ir0]#np.interp(p0,self.P[it0,:],self.U[it0,:])
          s0  = self.S[it0,ir0]#np.interp(p0,self.P[it0,:],self.S[it0,:])
          up0 = 0. # no initial particle velocity
          us0 = self.cs[it0,ir0]/1.e5 # cm/s->km/s use sound velocity for initial
          cs0 = self.cs[it0,ir0]/1.e5 # cm/s->km/s use sound velocity for initial
          #print(ir0,it0,r0,t0,p0,e0,up0,us0)
          self.hug.rho = np.append(self.hug.rho, self.rho[ir00])
          self.hug.P = np.append(self.hug.P, p0)
          self.hug.T = np.append(self.hug.T, self.T[it0])
          self.hug.U = np.append(self.hug.U, e0)
          self.hug.S = np.append(self.hug.S, s0)
          self.hug.up = np.append(self.hug.up, up0)
          self.hug.us = np.append(self.hug.us, us0)
          self.hug.cs = np.append(self.hug.cs, cs0)
          
          #for iir in range(ir0+1,self.ND):
          iir=ir00+1
          pnew=p0
          while pnew<pmax:
  #            ediff =0.5*(self.P[it0::,iir]+p0)*(1./r0-1./self.rho[iir])+e0 -(self.U[it0::,iir])  # MJ/kg
              ediff =0.5*(self.P[it0::,iir]+p0)*(1./r00-1./self.rho[iir])+e0 -(self.U[it0::,iir])  # MJ/kg
              # np.interp wants x values increasing
              pnew = np.interp(0.,np.flip(ediff),np.flip(self.P[it0::,iir]))
              tnew = np.interp(0.,np.flip(ediff),np.flip(self.T[it0::]))
              enew = np.interp(0.,np.flip(ediff),np.flip(self.U[it0::,iir]))
              snew = np.interp(0.,np.flip(ediff),np.flip(self.S[it0::,iir]))
              upnew = np.sqrt((pnew-p0)*(1./r00-1./self.rho[iir]))
              usnew = (1./r00)*np.sqrt((pnew-p0)/(1./r00-1./self.rho[iir]))
              csnew = np.interp(0.,np.flip(ediff),np.flip(self.cs[it0::,iir]))/1.E5 # km/s
              #print(self.rho[iir],tnew,pnew,enew,upnew,usnew)
              self.hug.rho = np.append(self.hug.rho, self.rho[iir])
              self.hug.P = np.append(self.hug.P, pnew)
              self.hug.T = np.append(self.hug.T, tnew)
              self.hug.U = np.append(self.hug.U, enew)
              self.hug.S = np.append(self.hug.S, snew)
              self.hug.up = np.append(self.hug.up, upnew)
              self.hug.us = np.append(self.hug.us, usnew)
              self.hug.cs = np.append(self.hug.cs, csnew) # km/s
              iir += 1
          self.hug.NH=len(self.hug.P)
          self.hug.units='units: T K, rho g/cm3, P GPa, U MJ/kg, S MJ/K/kg, Up km/s, Us km/s, cs km/s'
  
          if writefilename:
              print('Writing Hugoniot to file: ',writefilename)
              hugoniotfile = open(writefilename,"w")  
              hugoniotfile.writelines('  Hugoniot \n') 
              hugoniotfile.writelines('  Temperature,    Density,        Pressure,       Int. Energy,    Sp. Entropy,    Part. Vel.,     Shock Vel. \n') 
              hugoniotfile.writelines('  K,              g/cm3,          GPa,            MJ/kg,          MJ/K/kg,        km/s,           km/s\n') 
              for iih in range(0,self.hug.NH):
                  hugoniotfile.write("%14.6e, %14.6e, %14.6e, %14.6e, %14.6e, %14.6e, %14.6e\n" % (
                      self.hug.T[iih],self.hug.rho[iih],self.hug.P[iih],self.hug.U[iih],self.hug.S[iih],self.hug.up[iih],self.hug.us[iih]))
              hugoniotfile.close() 
  
      def writestdsesame(self, writestdsesfname=None):
          """Write standard Header-201-301 SESAME EOS TABULAR EOS FILE"""
          # write a standard SESAME ascii file
          #     WRITE STANDARD Header-201-301 SESAME FILE
          #     WRITE SESAME 301 TABLE CONTAINS P, E, HFE
          #sesfile = open("NEW-SESAME-STD-NOTENSION.EOSTXT","w")  
          if writestdsesfname is None:
              print('Please provide a file name.')
              exit(0)
          sesfile = open(writestdsesfname,"w")  
          #     WRITE SESAME HEADER INFORMATION: EOS matid number, number of words in section
          #     could input matid, date, version with the grid
          # these parameters are set in the cell above that sets up the grid for ANEOS
          # THEY SHOULD MATCH.......
          # These variables are needed for the standard table output
          NWDS=9
          SESNTABLES=2.0
          TABLE1 = 201.0
          TABLE2 = 301.0
          #     5 entries in 201 table
          SESNWDS1=5.0
          #     Number of entries in STANDARD 301 table: 3 variables at each T,rho point
          SESNWDS2=2.+self.ND+self.NT+self.ND*self.NT*3.
          #     HEADER SECTION
          #sesfile.write("%14.6e, %14.6e, %14.6e, %14.6e, %14.6e, %14.6e, %14.6e\n" % (antarr[iit],rnew,pnew,enew,snew,upnew,usnew))
          sesfile.write(" INDEX      MATID ={:7d}    NWDS = {:8d}\n".format(int(self.MATID), int(NWDS)))
          sesfile.write("{:16.8e}{:16.8e}{:16.8e}{:16.8e}{:16.8e}\n".format(self.MATID, self.DATE, self.DATE, self.VERSION, SESNTABLES))
          sesfile.write("{:16.8e}{:16.8e}{:16.8e}{:16.8e}\n".format(TABLE1, TABLE2, SESNWDS1, SESNWDS2))
          # 201 SECTION
          sesfile.write(" RECORD     TYPE ={:5d}     NWDS = {:8d}\n".format(int(TABLE1),int(SESNWDS1)))
          sesfile.write("{:16.8e}{:16.8e}{:16.8e}{:16.8e}{:16.8e}\n".format(self.FMN, self.FMW, self.R0REF, self.K0REF, self.T0REF))
          sesfile.write(" RECORD     TYPE ={:5d}     NWDS = {:8d}\n".format(int(TABLE2),int(SESNWDS2)))
          sesfile.write("{:16.8e}{:16.8e}".format(self.ND, self.NT))
          STYLE=2
          #     density array g/cm3
          for k in range(0, int(self.ND)):
              sesfile.write("{:16.8e}".format(self.rho[k]))
              STYLE=STYLE+1
              if (np.mod(STYLE,5) == 0):
                  sesfile.write("\n")
          #     temperature array K
          for j in range(0, int(self.NT)):
              sesfile.write("{:16.8e}".format(self.T[j]))
              STYLE=STYLE+1
              if (np.mod(STYLE,5) == 0):
                  sesfile.write("\n")
          #  pressure array GPa P[tempindex,dindex]
          for j in range(0,int(self.NT)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.P[j,k]))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          #  specific internal energy array MJ/kg U[tempindex,dindex]
          for j in range(0,int(self.NT)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.U[j,k]))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          # Helmholtz free energy array in MJ/kg A[tempindex,dindex]
          for j in range(0,int(self.NT)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.A[j,k]))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          # close the SESAME TABLE FILE
          sesfile.close() 
          print('Done writing the STD SESAME 301 notension table to local directory: ',writestdsesfname)
  
      def writemdqsesame(self, writemdqsesfname=None):
          """Function to write a sesame 301-style ascii file with the MDQ variable"""
          if writemdqsesfname is None:
              print('Please provide a file name.')
              exit(0)
          sesfile = open(writemdqsesfname,"w")  
          sesfile.write("{:16.8e}{:16.8e}".format(self.ND, self.NT))
          STYLE=2
          #     density array g/cm3
          for k in range(0, int(self.ND)):
              sesfile.write("{:16.8e}".format(self.rho[k]))
              STYLE=STYLE+1
              if (np.mod(STYLE,5) == 0):
                  sesfile.write("\n")
          #     temperature array K
          for j in range(0, int(self.NT)):
              sesfile.write("{:16.8e}".format(self.T[j]))
              STYLE=STYLE+1
              if (np.mod(STYLE,5) == 0):
                  sesfile.write("\n")
          #  MDQ Flag[tempindex,dindex]
          for j in range(0,int(self.NT)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.MDQ[j,k]))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          # close the SESAME TABLE FILE
          sesfile.close() 
          print('Done writing the MDQ Flag as a 301-style table to local directory: ',writemdqsesfname)
  
      def loadaneos(self, aneosinfname=None, aneosoutfname=None, silent=False):
          """Function for reading in ANEOS INPUT and OUTPUT FILE DATA into EOS structure."""
          if aneosinfname is None:
              return 'Must provide input file name.'
          if aneosoutfname is None:
              return 'Must provide output file name.'
          # function to gather data from ANEOS input and output files
          # SESAME FILE HEADER INFORMATION MUST BE LOADED INTO THE EOS STRUCTURE BEFORE CALLING THIS FUNCTION
          #
          # READ IN ANEOS INPUT FILE
          aneosinputfile = open(aneosinfname,"r")  
          testin=aneosinputfile.readlines()   # read in the whole ascii file at once because this is fatser
          aneosinputfile.close()
          # gather EOS information from the ANEOS.OUTPUT file
          aneosoutputfile = open(aneosoutfname,"r")  
          testout=aneosoutputfile.readlines() # read everything in at once because this is faster
          aneosoutputfile.close()
          if silent == False:
              print('Done loading ANEOS files.')
  
          # THIS CODE PARSES THE ANEOS.OUTPUT FILE INTO ARRAYS FOR USE IN PLOTTING/USING THE EOS
          if silent == False:
              print('ANEOS WAS CALLED WITH THE FOLLOWING INPUT, LOADED FROM FILE ',aneosinfname)
          # Gather parameters for the gamma function while printing the ANEOS INPUT FILE
          aneoscount=1
          C60=0.
          C61=0.
          betagamma=0.
          for i in np.arange(len(testin)):
              if testin[i].find('ANEOS') == 0:
                  if aneoscount<9:
                      if silent == False:
                          print(' '+testin[i-3],testin[i-2],testin[i-1],testin[i])
                      aneoscount=aneoscount+1
                  else:
                      if silent == False:
                          print(' '+testin[i])
                  if testin[i].find('ANEOS2') == 0:
                      tmp=testin[i]
                      nelem=int(tmp[10:20])
                      #print('nelem=',nelem)
                      rho0=float(tmp[30:40])
                      #print('rho0=',rho0)
                      gamma0=float(tmp[70:80])
                      #print('gamma0=',gamma0)
                      theta0=float(tmp[80:90])
                  if testin[i].find('ANEOS3') == 0:
                      tmp=testin[i]
                      C24=float(tmp[20:30])/3.
                      #print('C24=',C24)
                  if testin[i].find('ANEOS5') == 0:
                      tmp=testin[i]
                      C60=float(tmp[60:70])
                      C61=float(tmp[70:80])
                      #print('C60=',C60)
                  if testin[i].find('ANEOS7') == 0:
                      tmp=testin[i]
                      betagamma=float(tmp[70:80])
  
          # some checks
          if rho0 != self.R0REF:
              print('WARNING: rho0 does not match. STOPPING THIS NOTEBOOK.')
              assert(False) # just a way to stop the notebook
  
          # GUESS A BIG ARRAY SIZE FOR THE PHASE BOUNDARIES AND HUGONIOT IN ANEOS.OUTPUT
          # the melt curve, vapor curve and Hugoniot curves are not fixed length outputs
          nleninit=300
          meltcurve = 0
  
          if silent == False:
              print('READING DATA FROM ANEOS OUTPUT FILE ',aneosoutfname)
  
          # Read in data from the ANEOS.OUTPUT FILE
          imc = -1 # flag for no melt curve in the model
          for i in np.arange(len(testout)):
              if testout[i].find('  Data for ANEOS number') == 0:
                  tmp = testout[i+2][0:50]
                  eosname = tmp.strip()
              if testout[i] == '  TWO-PHASE BOUNDARIES\n':
                  nvc = nleninit
                  ivc = i
                  vcarrtmp = np.zeros((nvc,12),dtype=float)
                  flag=0
                  j=0
                  while flag == 0:
                      if testout[j+i+4].find(' anphas') == 0:
                          print(testout[j+i+4])
                          vcarrtmp[j,:]=vcarrtmp[j-1,:]
                          j=j+1
                      else:
                          tmp=str.replace(testout[j+i+4],'D','E')
                          tmp3 = tmp[0:157]
                          tmp4 = list(tmp3.split())
                          if (tmp4[3].find('E') == -1): # if the number is formatted badly, stop reading in the vapor curve, e.g. 1.00000-321
                              print('Stopped reading in the vapor curve because of badly formatted numbers. Stopped at this line:')
                              print(tmp3)
                              flag=1
                          else:
                              if (len(tmp4) >0) and (float(tmp4[3]) > 0) and (float(tmp4[4]) > 0): # stop if the pressures become negative on the vapor curve
                                  tmp5 = np.asarray(tmp4)
                                  vcarrtmp[j,:] = tmp5[:]
                                  j=j+1
                              else:
                                  flag=1
                  vcarr = np.zeros((j,12),dtype=float)
                  vcarr[:,:] = vcarrtmp[0:j,:]
              if testout[i] == ' LIQUID/SOLID PHASE CURVE\n':
                  nmc = nleninit
                  imc = i
                  meltcurve=1
                  mcarrtmp = np.zeros((nmc,11),dtype=float)
                  flag=0
                  j=0
                  while flag == 0:
                      tmp  = str.replace(testout[j+i+5],'D','E')
                      tmp3 = tmp[0:132]
                      tmp4 = list(tmp3.split())
                      if len(tmp4) > 0:
                          tmp5 = np.asarray(tmp4)
                          mcarrtmp[j,:] = tmp5[:]
                          j=j+1
                      else:
                          flag=1
                  mcarr = np.zeros((j,11),dtype=float)
                  mcarr[:,:] = mcarrtmp[0:j,:]
              if testout[i] == '   HUGONIOT\n':
                  nhc = nleninit
                  ihc = i
                  hcarrtmp = np.zeros((nhc,9),dtype=float)
                  flag=0
                  j=0
                  while flag == 0:
                      tmp=str.replace(testout[j+i+5],'D','E')
                      tmp3 = tmp[0:109]
                      tmp4 = list(tmp3.split())
                      if len(tmp4) > 0:
                          tmp4[3]='0.0' # this column often gives problems with exponential notation so don't read it
                          tmp5 = np.asarray(tmp4)
                          hcarrtmp[j,:] = tmp5[:]
                          j=j+1
                      else:
                          flag=1
                  hcarr = np.zeros((j,9),dtype=float)
                  hcarr[:,:] = hcarrtmp[0:j,:]
  
          # UPDATE THE MAIN EOS STRUCTURE WITH GATHERED INFORMATION
          # Add variables needed to calculate the ANEOS gamma function
          self.gamma0  = gamma0
          self.theta0  = theta0
          self.C24     = C24
          self.C60     = C60
          self.C61     = C61
          self.beta    = betagamma
          #
          # ANEOS.OUTPUT UNITS ARE NOT THE SAME AS THE SESAME TABLE!
          # add the vapor curve to this EOS object extracted from the ANEOS.OUTPUT
          #  TWO-PHASE BOUNDARIES
          #       T         RHOLIQ        RHOVAP        PLIQ         PVAP        ELIQ         EVAP         SLIQ         SVAP        GLIQ         GVAP         PSILIQ      PSIVAP         NTY
          #       K         kg/m**3       kg/m**3       GPa          GPa         J/kg         J/kg        J/kg-K       J/kg-K       J/kg         J/kg
          tmp = vcarr.shape
          #put vapor curve information in nicely named structure
          self.vc.NT  = tmp[0]
          self.vc.T   = vcarr[:,0] # K
          self.vc.rl  = vcarr[:,1]/1.E3 # g/cm3
          self.vc.rv  = vcarr[:,2]/1.E3 # g/cm3
          self.vc.Pl  = vcarr[:,3] # GPa
          self.vc.Pv  = vcarr[:,4] # GPa
          self.vc.Ul  = vcarr[:,5]/1.E6 # MJ/kg
          self.vc.Uv  = vcarr[:,6]/1.E6 # MJ/kg
          self.vc.Sl  = vcarr[:,7]/1.E6 # MJ/K/kg
          self.vc.Sv  = vcarr[:,8]/1.E6 # MJ/K/kg
          self.vc.Gl  = vcarr[:,9]/1.E6 # MJ/kg
          self.vc.Gv  = vcarr[:,10]/1.E6 # MJ/kg
          self.vc.units = 'T K, rho g/cm3, P GPa, U MJ/kg, S MJ/K/kg, G MJ/kg'
          # np.interp wants increasing x values
          self.onebar.Tvap = np.interp(1.E-4,np.flipud(self.vc.Pl),np.flipud(self.vc.T)) # extract boiling point temperature at 1 bar, K
          self.onebar.Siv  = np.interp(1.E-4,np.flipud(self.vc.Pl),np.flipud(self.vc.Sl)) # extract liquid sp. entropy at 1 bar, MJ/K/kg
          self.onebar.Scv  = np.interp(1.E-4,np.flipud(self.vc.Pl),np.flipud(self.vc.Sv)) # extract vapor sp. entropy at 1 bar, MJ/K/kg
          self.onebar.rhoiv  = np.interp(1.E-4,np.flipud(self.vc.Pl),np.flipud(self.vc.rl)) # extract liquid density at 1 bar, g/cm3
          self.onebar.rhocv  = np.interp(1.E-4,np.flipud(self.vc.Pl),np.flipud(self.vc.rv)) # extract vapor density at 1 bar, g/cm3
          #
          # add the ANEOS Hugoniot to this EOS object extracted from the ANEOS.OUTPUT
          #      RHO          T           P          PC           E           S           V           U       RHO/RHOO  #IT  STATE
          #    kg/m**3        K          GPa        GPa          J/kg      J/kg-K       km/sec      km/sec
          #self.anhug.all = hcarr  # 2D array of Hugoniot variables
          tmp = hcarr.shape
          self.anhug.ND = tmp[0] # number of density points on the Hugoniot
          self.anhug.rho = hcarr[:,0]/1.E3 # g/cm3
          self.anhug.T   = hcarr[:,1] # K
          self.anhug.P   = hcarr[:,2] # GPa
          self.anhug.U   = hcarr[:,4]/1.E6 # MJ/kg
          self.anhug.S   = hcarr[:,5]/1.E6 # MJ/K/kg
          self.anhug.us  = hcarr[:,6] # km/s
          self.anhug.up  = hcarr[:,7] # km/s
          self.anhug.units = 'vars: rho g/cm3, T K, P GPa, U MJ/kg, S MJ/K/kg, Us km/s, Up km/s'
          #
          # Add melt curve to EOS objects if available
          # LIQUID/SOLID PHASE CURVE
          #       T         RLIQ       RSOLID      PLIQ       PSOLID      ELIQ        ESOLID       SLIQ       SOLID        GLIQ       GSOLID        #ITER
          #       K        kg/m**3     kg/m**3      GPa         GPa       J/kg         J/kg        J/kg-K     J/kg-K       J/kg        J/kg
          if meltcurve == 1:
              # put the melt curve information in nicely named structure
              tmp=mcarr.shape
              self.mc.NT  = tmp[0]
              self.mc.T   = mcarr[:,0] # K
              self.mc.rl  = mcarr[:,1]/1.E3 # g/cm3
              self.mc.rs  = mcarr[:,2]/1.E3 # g/cm3
              self.mc.Pl  = mcarr[:,3] # GPa
              self.mc.Ps  = mcarr[:,4] # GPa
              self.mc.Ul  = mcarr[:,5]/1.E6 # MJ/kg
              self.mc.Us  = mcarr[:,6]/1.E6 # MJ/kg
              self.mc.Sl  = mcarr[:,7]/1.E6 # MJ/K/kg
              self.mc.Ss  = mcarr[:,8]/1.E6 # MJ/K/kg
              self.mc.units = 'T K, rho g/cm3, P GPa, U MJ/kg, S MJ/K/kg'
              # NOTE THAT TRIPLE POINT AND VAPOR CURVE SOLUTIONS DO NOT ALWAYS MATCH PERFECTLY AT THE TRIPLE POINT
              tmp = np.where(mcarr[:,3] > 0.)[0] # find the triple point first entry with positive pressure
              self.tp.T = mcarr[tmp[0],0] # K
              self.tp.P = np.interp(self.tp.T,np.flipud(self.vc.T),np.flipud(self.vc.Pv))   # this has trouble for forsterite; use the vapor size of the VC mcarr[:,3] # GPa
              self.tp.Sim  = mcarr[tmp[0],8]/1.E6 # extract solid sp. entropy at tp, MJ/K/kg
              self.tp.Scm  = mcarr[tmp[0],7]/1.E6 # extract liquid sp. entropy at tp, MJ/K/kg
              # np.interp wants x values increasing
              self.tp.Siv  = np.interp(self.tp.T,np.flipud(self.vc.T),np.flipud(self.vc.Sl)) # extract liquid sp. entropy at tp, MJ/K/kg
              self.tp.Scv  = np.interp(self.tp.T,np.flipud(self.vc.T),np.flipud(self.vc.Sv)) # extract vapor sp. entropy at tp, MJ/K/kg
              self.tp.rhol = mcarr[tmp[0],1]/1.E3 # extract liquid density at tp, g/cm3
              self.tp.rhos = mcarr[tmp[0],2]/1.E3 # extract solid density at tp, g/cm3
              self.tp.rhov = np.interp(self.tp.T,np.flipud(self.vc.T),np.flipud(self.vc.rv)) # extract vapor density at tp, g/cm3
              self.tp.units = 'T K, P GPa, S MJ/K/kg, rho g/cm3'
              # Extract melting point
              self.onebar.Tmelt = np.interp(1.E-4,self.mc.Pl,self.mc.T) # extract melting point temperature at 1 bar, K
              self.onebar.Sim   = np.interp(1.E-4,self.mc.Pl,self.mc.Ss) # extract liquid sp. entropy at 1 bar BP, MJ/K/kg
              self.onebar.Scm   = np.interp(1.E-4,self.mc.Pl,self.mc.Sl) # extract vapor sp. entropy at 1 bar BP, MJ/K/kg
              self.onebar.rhoim   = np.interp(1.E-4,self.mc.Pl[3::],self.mc.rs[3::]) # extract solid density at 1 bar MP, MJ/K/kg
              self.onebar.rhocm   = np.interp(1.E-4,self.mc.Pl[3::],self.mc.rl[3::]) # extract liquid density at 1 bar MP, MJ/K/kg
          # put the data for the critical point in the EOS structure for easy access
          self.cp.T   = vcarr[0,0] # K
          self.cp.rho = vcarr[0,1]/1.E3 # g/cm3
          self.cp.P   = vcarr[0,3] # GPa
          self.cp.U   = vcarr[0,5]/1.E6 # MJ/kg 
          self.cp.S   = vcarr[0,7]/1.E6 # MJ/K/kg
          self.cp.units = 'T K, rho g/cm3, P GPa, U MJ/kg, S MJ/K/kg'
          #------------------------------------------------------------------
  
  ###########################################################################################
  ########## GADGET STYLE TABLES
  ###########################################################################################
  class extGADtable:
      """Class for v2019 GADGET EOS tables using SESAME UNITS"""
      def __init__(self):
          """Function to initialize GADGET TABULAR EOS structure"""
          self.ND  = 0 # integer; number of density points in grid
          self.NS  = 0 # integer; number of sp. entropy points in grid
          self.rho = np.zeros(self.ND)          # g/cm3, density values
          self.S   = np.zeros(self.NS)          # MJ/K/kg, sp. entropy values
          self.P   = np.zeros(self.ND*self.NS)  # GPA, pressure(S,rho)
          self.T   = np.zeros(self.ND*self.NS)  # K, temperature(S,rho)
          self.U   = np.zeros(self.ND*self.NS)  # MJ/kg, sp. internal energy(S,rho)
          self.A   = np.zeros(self.ND*self.NS)  # MJ/kg, Helmholtz free energy(S,rho)
          self.cs  = np.zeros(self.ND*self.NS)  # cm/s, sound speed(S,rho)
          self.cv  = np.zeros(self.ND*self.NS)  # MJ/K/kg, sp. heat capacity(S,rho)
          self.KPA = np.zeros(self.ND*self.NS)  # integer, ANEOS KPA flag(S,rho)
          self.MDQ = np.zeros(self.ND*self.NS)  # integer, Model Development Quality Flag(S,rho)
          self.units = 'Units: g/cm3, MJ/K/kg, GPa, K, MJ/kg, MJ/kg, cm/s, MJ/K/kg, KPA flag. 2D arrays are (NS,ND).'
          self.MODELNAME = ''
      #
      def view(self, q='P', Slow=None, Shigh=None, rholow=None, rhohigh=None):
          """Function to print variables from GADGET EOS table"""
          if Slow is None:
              Slow = self.S.min()
          if Shigh is None:
              Shigh = self.S.max()
          if rholow is None:
              rholow = self.rho.min()
          if rhohigh is None:
              rhohigh = self.rho.max()
          if q == 'P':
              print('S:', self.S[np.logical_and(self.S >= Slow,self.S <= Shigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho <= rhohigh)
                                    ])
              print('P:', (self.P[np.logical_and(self.S >= Slow,
                                                  self.S<=Shigh)
                                 ])[:, np.logical_and(self.rho >= rholow,
                                                       self.rho <= rhohigh)
                                   ])
          if q == 'U':
              print('S:', self.S[np.logical_and(self.S >= Slow,self.S <= Shigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho <= rhohigh)
                                    ])
              print('U:', (self.U[np.logical_and(self.S >= Slow,
                                                  self.S <= Shigh)
                                 ])[:, np.logical_and(self.rho >= rholow,
                                                       self.rho <= rhohigh)
                                   ])
          if q == 'A':
              print('S:', self.S[np.logical_and(self.S >= Slow,self.S <= Shigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho <= rhohigh)
                                    ])
              print('A:', (self.A[np.logical_and(self.S >= Slow,
                                                  self.S <= Shigh)
                                 ])[:, np.logical_and(self.rho >= rholow,
                                                       self.rho <= rhohigh)
                                   ])
          if q == 'T':
              print('S:', self.S[np.logical_and(self.S >= Slow,self.S <= Shigh)])
              print('rho:', self.rho[np.logical_and(self.rho
                                     >= rholow,self.rho<=rhohigh)
                                    ])
              print('T:', (self.T[np.logical_and(self.S >= Slow,self.S <= Shigh)
                                 ])[:, np.logical_and(self.rho >= rholow,
                                                       self.rho <= rhohigh)
                                   ])      
          if q == 'cs':
              print('S:', self.S[np.logical_and(self.S >= Slow,self.S <= Shigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho<=rhohigh)])
              print('cs:', (self.cs[np.logical_and(self.S >= Slow,
                                                    self.S <= Shigh)
                                   ])[:, np.logical_and(self.rho >= rholow,
                                                         self.rho <= rhohigh)
                                     ])
          if q == 'cv':
              print('S:', self.S[np.logical_and(self.S >= Slow,self.S <= Shigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho<=rhohigh)])
              print('cv:', (self.cv[np.logical_and(self.S >= Slow,
                                                    self.S <= Shigh)
                                   ])[:, np.logical_and(self.rho >= rholow,
                                                         self.rho <= rhohigh)
                                     ])
          if q == 'KPA':
              print('S:', self.S[np.logical_and(self.S >= Slow,self.S <= Shigh)])
              print('rho:', self.rho[np.logical_and(self.rho >= rholow,
                                                     self.rho<=rhohigh)])
              print('KPA:', (self.KPA[np.logical_and(self.S >= Slow,
                                                    self.S <= Shigh)
                                   ])[:, np.logical_and(self.rho >= rholow,
                                                         self.rho <= rhohigh)
                                     ])
      #
      def writestdgadget(self, writestdgadgetfname=None):
          """Function to write an standard gadget2 EOS ascii file in 301 format"""
          #sesfile = open("NEW-GADGET-STD-NOTENSION.TXT","w")  
          if writestdgadgetfname is None:
              print('Please provide a file name.')
              exit(0)
          sesfile = open(writestdgadgetfname,"w")  
          sesfile.write("{:16.8e}{:16.8e}".format(self.ND, self.NS))
          STYLE=2
          #     density array g/cm3
          for k in range(0, int(self.ND)):
              sesfile.write("{:16.8e}".format(self.rho[k]))
              STYLE=STYLE+1
              if (np.mod(STYLE,5) == 0):
                  sesfile.write("\n")
          #     sp. entropy array MJ/K/kg -> erg/K/g
          for j in range(0, int(self.NS)):
              sesfile.write("{:16.8e}".format(self.S[j]*1.E10))
              STYLE=STYLE+1
              if (np.mod(STYLE,5) == 0):
                  sesfile.write("\n")
          #  pressure array GPa -> dynes/cm2 P[sindex,dindex]
          for j in range(0,int(self.NS)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.P[j,k]*1.E10))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          #  temperature array K T[sindex,dindex]
          for j in range(0,int(self.NS)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.T[j,k]))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          #  specific internal energy array MJ/kg -> erg/g U[sindex,dindex]
          for j in range(0,int(self.NS)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.U[j,k]*1.E10))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          #  sound speed array cm/s cs[sindex,dindex]
          for j in range(0,int(self.NS)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.cs[j,k]))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          # close the STD GADGET TABLE FILE
          sesfile.close() 
          print('Done writing the STD GADGET TABLE FILE with notension to local directory: ',writestdgadgetfname)
      #
      def writeextgadget(self, writeextgadgetfname=None):
          """Function to write an extended gadget2 EOS ascii file in 301 format"""
          #sesfile = open("NEW-GADGET-EXT-NOTENSION.TXT","w")  
          if writeextgadgetfname is None:
              print('Please provide a file name.')
              exit(0)
          sesfile = open(writeextgadgetfname,"w")  
          sesfile.write("{:16.8e}{:16.8e}".format(self.ND, self.NS))
          STYLE=2
          #     density array g/cm3
          for k in range(0, int(self.ND)):
              sesfile.write("{:16.8e}".format(self.rho[k]))
              STYLE=STYLE+1
              if (np.mod(STYLE,5) == 0):
                  sesfile.write("\n")
          #     sp. entropy array MJ/K/kg -> erg/K/g
          for j in range(0, int(self.NS)):
              sesfile.write("{:16.8e}".format(self.S[j]*1.E10))
              STYLE=STYLE+1
              if (np.mod(STYLE,5) == 0):
                  sesfile.write("\n")
          # Helmholtz free energy array in MJ/kg -> ergs/g A[tempindex,dindex]
          for j in range(0,int(self.NS)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.A[j,k]*1.E10))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          # sp. heat capacity array MJ/K/kg -> erg/K/g cv[sindex,dindex]
          for j in range(0,int(self.NS)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.cv[j,k]*1.E10))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          # KPA array flag k[sindex,dindex]
          for j in range(0,int(self.NS)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.KPA[j,k]))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          # MDQ array flag MDQ[sindex,dindex]
          for j in range(0,int(self.NS)):
              for k in range(0,int(self.ND)):
                  sesfile.write("{:16.8e}".format(self.MDQ[j,k]))
                  STYLE=STYLE+1
                  if (np.mod(STYLE,5) == 0):
                      sesfile.write("\n")
          # close the EXTENDED GADGET TABLE FILE
          sesfile.close() 
          print('Done writing the EXTENDED GADGET TABLE FILE with notension to local directory: ',writeextgadgetfname)
  #
  ###########################################################################################
  ########## TILLOTSON TABLE CLASS
  ###########################################################################################
  class TillotsonClass:
      """Class for tabulated Tillotson EOS tables using SESAME UNITS"""
      def __init__(self):
          """Function to initialize Tillotson TABULAR EOS structure"""
          # variables stored in SESAME units
          self.rho0  = 0. # reference density g/m3
          self.E0    = 0. # E0 constant MJ/kg
          self.EIV   = 0. # specific internal energy of incipient vaporization MJ/kg
          self.ECV   = 0. # specific internal energy of complete vaporization MJ/kg
          self.AA    = 0. # Bulk modulus K0 GPa
          self.BB    = 0. # B constant GPa
          self.a     = 0. # a constant [-]
          self.b     = 0. # b constant [-]
          self.alpha = 0. # alpha constant [-]
          self.beta  = 0. # beat constant [-]    a+b=Gruneisen gamma at rho0
          self.ND    = 0 # integer; number of density points in grid
          self.NU    = 0 # integer; number of sp. internal energy points in grid
          self.rho   = np.zeros(self.ND)          # g/cm3, density values
          self.U     = np.zeros(self.NU)          # MJ/kg, sp. internal energy values
          self.P     = np.zeros(self.ND*self.NU)  # GPA, pressure(U,rho)
          self.cs    = np.zeros(self.ND*self.NU)  # km/s, sound speed(U,rho)
          self.region= np.zeros(self.ND*self.NU)  # region flag: 1-condensed, 2-interpolated, 3-expanded
          # self.T add later
          self.units = 'Units: g/cm3, MJ/kg, GPa, cm/s. 2D arrays are (NU,ND).'
          self.MODELNAME = ''
          self.hug = EOShugoniot()
  
      def FillTable(self, matparams=None, modelname=None):
          """Function to populate table with Tillotson EOS. Assumes ND,NU,rho,and U are filled."""
          if matparams is None:
              print('Must provide a material parameter list [rho0, E0, EIV, ECV, AA, BB, a, b, alpha, beta].')
              exit(0)
          # matparams: [rho0, E0, EIV, ECV, AA, BB, a, b, alpha, beta]
          # units:     [kg/m3, J/kg, J/kg, J/kg, Pa, Pa, [-]x4]
          # print(self.NU,self.ND)
          self.cs = np.zeros((self.NU,self.ND))
          self.P  = np.zeros((self.NU,self.ND))
          self.region = np.zeros((self.NU,self.ND))
  #        print('Using the iSALE Tillotson implementation.')
          print('Using standard Tillotson implementation (e.g., Hosono et al. 2019, Melosh 1989).')
          print('Tillotson parameters: ',matparams)
          for ie in range(0,self.NU):
              for ir in range(0,self.ND):
                  flag               = 0
                  eng                = self.U[ie]*1.e6 # J/kg
                  dens               = self.rho[ir]*1000. # kg/m3
                  tmp                = Till_P_Hosono(dens,eng,matparams) # returns [P Pa, flag]
                  self.P[ie,ir]      = tmp[0]/1.e9 # Pa to GPa
                  self.region[ie,ir] = tmp[1] # cold, hot, interpolated regions
                  self.cs[ie,ir]     = Till_SoundSpeed(dens,eng,matparams)/1.E3 # km/s ; function should return m/s but it's not *** CHECK
                  #tmp                = Till_P(dens,eng,matparams) # returns [P Pa, flag, cs]
                  #self.P[ie,ir]      = tmp[0]/1.e9 # Pa to GPa
                  #self.region[ie,ir] = tmp[1] # cold, hot, interpolated regions
                  #self.cs[ie,ir]     = tmp[2]/1.E3 # km/s
  
      def calchugoniot(self, r00=None, e0=None, pmax=None, writefilename=None):
          """Function for calculating a Hugoniot from Tillotson EOS table."""
          if r00 is None:
              return 'Must provide r00 in g/cm^3 and e0 in MJ/kg.'
          if e0 is None:
              return 'Must provide r0 in g/cm^3 and e0 in MJ/kg.'
          if pmax is None:
              pmax=1.E4 # GPa
          self.hug.rho = []
          self.hug.P = []
          self.hug.T = []
          self.hug.U = []
          self.hug.up = []
          self.hug.us = []
          self.hug.cs = [] 
          self.hug.region = [] # Tillotson region flag
  
          ie0 = int(np.round(np.interp(e0,self.U,np.arange(self.NU)))) # uses nearest value if e0 not in array
          ir0 = int(np.round(np.interp(r00,self.rho,np.arange(self.ND)))) # uses nearest value if r0 not in the array
          r0  = self.rho[ir0] # g/cm^3
          p0  = self.P[ie0,ir0] # GPa
          e0  = self.U[ie0] # MJ/kg
          up0 = 0. # no initial particle velocity km/s
          us0 = self.cs[ie0,ir0] # km/s use sound velocity for initial
          cs0 = self.cs[ie0,ir0] # km/s use sound velocity for initial
          print('calchugoniot:',ir0,ie0,r0,p0,e0,up0,us0)
          self.hug.rho = np.append(self.hug.rho, self.rho[ir0])
          self.hug.P = np.append(self.hug.P, p0)
          self.hug.U = np.append(self.hug.U, e0)
          self.hug.up = np.append(self.hug.up, up0)
          self.hug.us = np.append(self.hug.us, us0)
          self.hug.cs = np.append(self.hug.cs, cs0)
          
          #for iir in range(ir0+1,self.ND):
          iir=ir0+1
          pnew=p0
          while pnew<pmax:
              ediff =0.5*(self.P[ie0::,iir]+p0)*(1./r0-1./self.rho[iir])+e0 -(self.U[ie0::])  # MJ/kg
              # np.interp wants x values increasing
              pnew = np.interp(0.,np.flip(ediff),np.flip(self.P[ie0::,iir]))
              enew = np.interp(0.,np.flip(ediff),np.flip(self.U[ie0::]))
              csnew = np.interp(0.,np.flip(ediff),np.flip(self.cs[ie0::,iir]))
              upnew = np.sqrt((pnew-p0)*(1./r0-1./self.rho[iir]))
              usnew = (1./r0)*np.sqrt((pnew-p0)/(1./r0-1./self.rho[iir]))
              #print(self.rho[iir],pnew,enew,upnew,usnew)
              self.hug.rho = np.append(self.hug.rho, self.rho[iir])
              self.hug.P = np.append(self.hug.P, pnew)
              self.hug.U = np.append(self.hug.U, enew)
              self.hug.up = np.append(self.hug.up, upnew)
              self.hug.us = np.append(self.hug.us, usnew)
              self.hug.cs = np.append(self.hug.cs, csnew)
              iir += 1
          self.hug.NH=len(self.hug.P)
          self.hug.units='units: rho g/cm3, P GPa, U MJ/kg, Up km/s, Us km/s, cs km/s'
          print('Done calculating Hugoniot with Tillotson EOS table.')
  #
  ###########################################################################################
  ########## TILLOTSON EOS FUNCTIONS
  ###########################################################################################
  # TILLOTSON FUNCTIONS AS IMPLEMENTED BY HOSONO
  # dunite tillotson parameters used by Hosono et al. 2019 
  # these olivine parameters are from Marinova et al. 2011 Icarus 
  # parameters: [rho0, E0, EIV, ECV, AA, BB, a, b, alpha, beta]
  # units:    [kg/m3, J/kg, J/kg, J/kg, Pa, Pa, [-]x4]
  # dunitetill = [3500.0, 550.0e+6, 4.500e+6, 14.50e+6, 131.00e+9,  49.00e+9, 0.5, 1.4, 5.0, 5.0]
  # Basalt parameters from iSALE -- from where? Benz?
  # basalttill = [2650.0, 4.87E8, 4.72E6, 18.2E6, 5.3E10, 5.3E10, 0.6, 0.6, 5., 5.]
  #
  def Till_P_co(dens,eng,tilleos):
      r0 = 0
      E0 = 1
      EIV = 2
      ECV=3
      AA=4
      BB=5
      a=6
      b=7
      alpha=8
      beta=9
      # Hosono implementation - straight up equation
      eta  = dens/tilleos[r0]
      mu   = eta - 1.0
      pco = (tilleos[a] + tilleos[b] / (eng / tilleos[E0] / eta / eta + 1.0)) * dens * eng + \
           tilleos[AA] * mu + tilleos[BB] * mu * mu
  
      return pco
  #
  def Till_P_ex(dens,eng,tilleos):
      import numpy as np
      r0 = 0
      E0 = 1
      EIV = 2
      ECV=3
      AA=4
      BB=5
      a=6
      b=7
      alpha=8
      beta=9
      # hosono equation
      eta  = dens/tilleos[r0]
      mu   = eta - 1.0
      pex = tilleos[a] * dens * eng + \
            (tilleos[b] * dens * eng / (eng / tilleos[E0] / eta / eta + 1.0) + \
             tilleos[AA] * mu * np.exp(- tilleos[beta] * (1.0 / eta - 1.0))) * \
                                         np.exp(- tilleos[alpha] * (1.0 / eta - 1.0) * (1.0 / eta - 1.0))
      return pex
  #
  def Till_P_Hosono(dens,eng,tilleos):
      flag=0 # flag for the region of the EOS
      r0 = 0 # index numbers for variables in the tilleos array
      E0 = 1
      EIV = 2
      ECV=3
      AA=4
      BB=5
      a=6
      b=7
      alpha=8
      beta=9
  
      # Hosono has 2 options
      pmin=1.e-7
      
      # option 1
      if ( (dens>=tilleos[r0]) or (eng < tilleos[EIV]) ):
          # condensed region
          flag=1
          pco = Till_P_co(dens,eng,tilleos)
          if pco < pmin:
              pout = pmin
          else:
              pout = pco
          return [pout,flag]
      else:
          if ( (dens < tilleos[r0]) and (eng > tilleos[ECV]) ):
              # expanded region
              flag=3
              pex = Till_P_ex(dens,eng,tilleos)
              if pex < pmin:
                  pout = pmin
              else:
                  pout = pex
              return [pout,flag]
          else:
              # interpolated region
              flag=2
              Pex = Till_P_ex(dens,eng,tilleos)
              Pco = Till_P_co(dens,eng,tilleos)
              # nobody seem to do this step
              #if Pco < pmin:
              #    Pco = pmin
              Pint = ((eng-tilleos[EIV])*Pex+(tilleos[ECV]-eng)*Pco)/(tilleos[ECV]-tilleos[EIV])
              if Pint < pmin:
                  pout = pmin
              else:
                  pout = Pint
              return [pout,flag]
          
      '''    
      # option 2
      if ( (dens>=tilleos[r0]) or (eng < tilleos[EIV]) ):
          # condensed region
          flag=1
          pout = Till_P_co(dens,eng,tilleos)
          if (dens <= 0.9*tilleos[r0]):
              pout = 1.e-16
              return [pout,flag]
      else:
          if ( (dens < tilleos[r0]) and (eng > tilleos[ECV]) ):
              # expanded region
              flag=3
              pout = Till_P_ex(dens,eng,tilleos)
          else:
              # interpolated region
              flag=2
              Pex = Till_P_ex(dens,eng,tilleos)
              Pco = Till_P_co(dens,eng,tilleos)
              pout = ((eng-tilleos[EIV])*Pex+(tilleos[ECV]-eng)*Pco)/(tilleos[ECV]-tilleos[EIV])
      if (pout < pmin):
          pout = pmin
      return [pout,flag]
      '''
      # END HOSONO TILLOTSON FUNCTION
  
  def Till_dPdrho(dens,eng,tilleos):
      drho=0.0001
      a = Till_P_Hosono(dens+drho,eng,tilleos)
      b = Till_P_Hosono(dens-drho,eng,tilleos)
      return  (a[0] - b[0])/(2.*drho)
  
  def Till_dPdu(dens,eng,tilleos):
      ddu=0.0001
      a=Till_P_Hosono(dens,eng+ddu,tilleos) 
      b=Till_P_Hosono(dens,eng-ddu,tilleos) 
      return (a[0] - b[0])/(2.*ddu)
  
  def Till_SoundSpeed(dens,eng,tilleos):
      if dens<= 100.:
          cs = 1.e-8 # m/s
      else:
          tmp = Till_P_Hosono(dens,eng,tilleos)
          cs = tmp[0]/(dens*dens)*Till_dPdu(dens,eng,tilleos) + Till_dPdrho(dens,eng,tilleos)
          if cs < 0.:
              cs = 0.
          cs = np.sqrt(cs + 1.e-16)
      return cs
  # END HOSONO FUNCTIONS
  # ------------------------------
  # iSALE implementation of Tillotson
  def Till_P(rin,ein,tilleos):
      # returns [pout,flag,csout]
      # output variables
      flag=0 # flag for the region of the EOS
      pout=0.
      csout=0.
      # tout = 0. # add temperature later
      #
      # iSALE implementation of Tillotson
      r0 = 0 # index numbers for variables in the tilleos array
      E0 = 1
      EIV = 2
      ECV=3
      AA=4
      BB=5
      a=6
      b=7
      alpha=8
      beta=9
      #
      # if density is zero or negative return 0 pressure
      if (rin <= 0.):
          pout  = 0.
          csout = 0.
          flag  = -1
          return [pout,flag,csout]
      #
      # define intermediate variables
      eta  = rin/tilleos[r0]
      mu   = eta - 1.0
      imu  = tilleos[r0]/rin - 1.0
      temp = tilleos[E0] * np.power(eta,2.)
      if (temp > 0.):
          erel = ein/temp
      else:
          erel = 0.
      ierel = 1./(erel + 1.)
      grun  = tilleos[a] + tilleos[b]*ierel
      gterm = tilleos[a] * rin * ein
      eterm = tilleos[b] * rin * ein * ierel
      #
      # initialize hot and cold terms to the pressure
      ph = 0.
      pc = 0.
      #
      if ( (eta < 1.) and (ein >= tilleos[EIV]) ):
          # expanded region above EIV
          flag = 3
          tmp  = tilleos[alpha]*np.power(imu,2.)
          if (tmp < 100.):
              exp_a = np.exp(-tmp)
          else:
              exp_a = 0.
          tmp = tilleos[beta]*imu
          if (tmp < 100.):
              exp_b = np.exp(-tmp)
          else:
              exp_b = 0.
          # compute "hot pressure"
          ph  = gterm + (eterm + tilleos[AA]*mu*exp_b)*exp_a  
          #
          # compute sound speed from derivatives
          dpre2   = tilleos[b]*ein*ierel*(1.+2.*(tilleos[alpha]*imu/eta + erel*ierel))*exp_a
          dpre3   = (tilleos[AA]/rin)*(1.+(mu/np.power(eta,2.))*(tilleos[beta]+2.*tilleos[alpha]*imu))*exp_a*exp_b
          dpreh   = tilleos[a]*ein+dpre2+dpre3
          dperh   = rin*(tilleos[a]+tilleos[b]*ierel*(1.-erel*ierel)*exp_a)
          c_dperh = (ph/np.power(rin,2.))*dperh
          cs2_h   = dpreh + c_dperh
          cs2     = cs2_h 
      #
      if ( (eta >= 1.) or (ein < tilleos[ECV]) ):
          # compressed region or expanded below energy of complete vaporization
          flag = 1
          if (eta >= 1.):
              b_temp = tilleos[BB] # in compression
          else:
              b_temp = tilleos[AA] # in expansion -- why does iSALE do this?
          pcold = (tilleos[AA] * mu) + (b_temp * np.power(mu,2.))
          pc    = gterm + eterm + pcold
          #if (pc < 0.):
          #    pc=1.e-10
          #
          # compute sound speed from derivatives
          dprec   = (tilleos[AA]+2.*b_temp*mu)/tilleos[r0] + ein*(grun+2.*tilleos[b]*erel*ierel*ierel)
          dperc   = rin * (grun-tilleos[b]*erel*ierel*ierel)
          c_dperc = (pc/(np.power(rin,2.)))*dperc
          cs2_c   = dprec + c_dperc
          cs2     = cs2_c
      #
      pout = pc + ph
      #
      if (ph*pc != 0.):    
          # interpolated region
          flag=2
          pout = ((ein-tilleos[EIV])*ph+(tilleos[ECV]-ein)*pc)/(tilleos[ECV]-tilleos[EIV])
          # transition sound speed
          dprem = ((ein-tilleos[EIV])*dpreh + (tilleos[ECV]-ein)*dprec) / (tilleos[ECV]-tilleos[EIV])
          dperm = ((ein-tilleos[EIV])*dperh + (tilleos[ECV]-ein)*dperc+ph-pc) / (tilleos[ECV]-tilleos[EIV])
          cs2_m = dprem+(pout/np.power(rin,2.))*dperm
          cs2   = cs2_m
      #
      # limiting after tillotson EOS
      if (pout <= 1.e-10):
          pout = 1.e-10 # pmin
          cs2  = tilleos[AA]/tilleos[r0]
      #
      if cs2>0:
          csout = np.sqrt(cs2)
      else:
          csout = np.sqrt(tilleos[AA]/tilleos[r0])
      return [pout,flag,csout]
      # END OF TILLOTSON EOS ISALE IMPLEMENTATION FUNCTIONS
  ####################################################################################
  ### END eostable.py FILE ###
  ####################################################################################
  