c f77 msmfunctions.f -o msmfunctions
c
c see appendices A, B  of Martin, Schwarz, and Mandy 1996, ApJ, 461, 265.

c input file: contable.f
c
c interactive input of temperature and density and selection of 
c desired function (by an integer 1--5 corresponding to the order of
c the gammas in Table 1)
c
c sample driver program to read contable.dat and evaluate the functions
c
      character*5 namprog(5),namsubs(5),namgrid(5)
      real*8 ff
      dimension rm(5),fx(5),fn(5)
      dimension ym(21,5),y(21)
      n=21
c
      write(6,*) 
     $     ' Evaluate the functions in Appendices A and B of'
      write(6,*) 
     $     ' Martin, Schwarz, and Mandy 1996, ApJ, 461, 265'
      write(6,*) 
     $     ' using the constants in Table 1'
      write(6,*) ' '
c
c read the data from the common table
c
 2    format(1x, 5(1x,2x,a5,7x))
 4    format(1x, 5(1x,1pd14.6))
      open(unit=16,file='contable.dat',status='old')
      write(6,*) ' reading contable.dat'
      write(6,*) ' '
c
c skim off the top titles
c
      read(16,2) (namsubs(j),j=1,5)
      read(16,2) (namprog(j),j=1,5)
      read(16,2) (namgrid(j),j=1,5)
c
c skim off the rms values
c
      read(16,4) (rm(j),j=1,5)
      read(16,4) (fx(j),j=1,5)
      read(16,4) (fn(j),j=1,5)
c
c get the array of 21 coefficients
c
      do 1 i=1,n
         read(16,4) (ym(i,j),j=1,5)
 1    continue
c
      write(6,*) ' begin interactive input'
      write(6,*) ' '
c
c loop back to here
c
 20   continue
c
c
c select which function
c
      write(6,*) ' choose one of the 5 functions (negative to quit)'
      write(6,*) ' po=1, cid=2, dt=3, ie=4, r=5'
c     write(6,*) ' '
      read(5,*) j
      if(j.lt.1) go to 30
      if(j.gt.5) go to 20
c
c choose appropriate constants from array
c
      do 10 i=1,n
         y(i)=ym(i,j)
 10   continue
c
c loop back to here
c
 50   continue
c
c
c temperature and density
c
      write(6,*) ' input temperature and density (cm-3)'
      write(6,*) ' (choose negative temperature to get new function)'
      read(5,*) te,de
      if(te.lt.0) go to 20
c
c check that temperature and density are in range over which fit is valid
c for interpolation (avoid extrapolation)
c
      if(j.eq.1) then
         if(te.lt.450. .or. te.gt.45000.) 
     $        write(6,*) 
     $ ' extrapolation warning: temperature out of range 450 - 45000'
      else
         if(te.lt.600. .or. te.gt.45000.) 
     $        write(6,*) 
     $ ' extrapolation warning: temperature out of range 600 - 45000'
      endif
      if(de.lt.1.d-3 .or. de.gt.1.e9) 
     $     write(6,*) 
     $ ' extrapolation warning: density out of range 10-3 - 10+9'
c
c evaluate the function
c
      if(j.eq.5) then
         f=ftdr(19,y,te,de)
      else
         f=ftd(n,y,te,de)
      endif
      ff=10.d0**f
      write(6,40) namsubs(j),f,ff
 40   format(1x,' log and value of function ',a5,1x,f10.4,1pd13.4)
      write(6,*) ' '
c
      go to 50
c
 30   continue
      write(6,*) ' goodbye'
      stop
      end
c
      function ftd(n,y,te,de)
c evaluate log_f at a single te, de
c po cid dt ie only
c
      dimension y(n)
c
      tl=alog10(te)
c high
      a =   y(1) 
     $     +y(2)*tl
     $     +y(3)*tl*tl
     $     +y(4)*tl*tl*tl
     $     +y(5)*alog10(1.0
     $     +y(6)/te)
      a1=   y(7)/te
c low
      b =   y(8)
     $     +y(9)*tl
     $     +y(10)*tl*tl
     $     +y(11)*alog10(1.0
     $     +y(12)/te)       
      b1=   y(13)/te
c critical
      c =   y(14) 
     $     +y(15)*tl
     $     +y(16)*tl*tl
     $     +y(17)/te
      c1 =  y(18)+c
c exponent
      d =   y(19)
     $     +y(20)*exp(-te/1850.0)
     $     +y(21)*exp(-te/440.0)    
      ftd = a- (a-b)  /(1.0+(de/10.**c )**d)
     $     +a1-(a1-b1)/(1.0+(de/10.**c1)**d)
c
      return
      end
c
      function ftdr(n,y,te,de)
c evaluate log_f at a single te, de
c ra only
      dimension y(n)
c
      tl=alog10(te)
      dl=alog10(de)
c     
c high (lambda form)
      a =   y(1) 
     $     +y(2)*tl
     $     +y(3)*tl*tl
     $     +y(4)*tl*tl*tl
     $     +y(5)*alog10(1.0
     $     +y(6)/te)
     $     +y(7)/te
c low
      b =   y(8) 
     $     +y(9)*tl
     $     +y(10)*tl*tl
     $     +y(11)*tl*tl*tl
     $     +y(12)*alog10(1.0
     $     +y(13)/te)
     $     +y(14)/te
c log(n_H / n_crit) == log(n_h*low/high)
c log density relative to vertex
      f = dl - b 
c general purpose constant
      w=    y(15)
     $     +y(16)*tl
     $     +y(17)*tl*tl
     $     +y(18)*tl*tl*tl
     $     +y(19)/te
c combine to make hyperbola
      ftdr = a + 0.5*(f - sqrt(f*f+2.*w*w))
c change to gamma
      ftdr = ftdr -dl
      return
      end
c
