C*********************************************************************************
C*                                                                               *
C*       First principle equation of state of hot, dense helium                  *
C*                                                                               *
C*       Derived by combining two first-principles techiques                     *
C*            a) path integral Monte Carlo and                                   *
C*            b) density functional molecular dynamics                           *
C*                                                                               *
C*       Free energy interpolation code written by B. Militzer                   *
C*                                                                               *
C*       Copy and usage rights by B. Militzer          Berkeley, CA, 02-01-08    *
C*                                                                               *
C*********************************************************************************

      real*8 FUNCTION AUToK(T)
      implicit none
      real*8 T
      AUToK = 315774.660873d0 * T
      return
      end

      real*8 FUNCTION KToAU(TinK)
      implicit none
      real*8 TinK,AUToK
      KToAU = TinK / AUToK(1.0d0)
      return
      end
    
      real*8 FUNCTION RsFromDensity(nn)         
      implicit none
      real*8 nn,rs,pi
      pi = 3.14159265358979323846d0
      rs = (3.0d0/(4.0d0*pi*nn))**(1.0d0/3.0d0)
      RsFromDensity = rs
      return
      end

      real*8 FUNCTION DensityFromRs(rs)         
      implicit none
      real*8 v,rs,pi
      pi = 3.14159265358979323846d0
      v  = 4.0d0/3.0d0*pi*rs*rs*rs
      DensityFromRs = 1.0d0/v;
      return
      end

      real*8 FUNCTION DensityToGCC(nn)         
      implicit none
      real*8 nn
C* the difference between H and He enters here
      DensityToGCC = 22.4263095064477d0 * nn 
      return
      end

      real*8 FUNCTION DensityFromGCC(nn)         
      implicit none
      real*8 DensityToGCC,nn
      DensityFromGCC = nn/DensityToGCC(1.0d0)
      return
      end

      real*8 FUNCTION PressureToGPa(P)         
      implicit none
      real*8 P
      PressureToGPa = 29421.009874d0*P
      return
      end

      real*8 FUNCTION EnergyToeV(E)         
      implicit none
      real *8 E
      EnergyToeV = E * 27.2113834279d0
      return
      end

C****************************************************************************

      real*8 FUNCTION FreeEnergy(nn,T)         
      implicit none
      
      integer i,imin,imax,n,NMAX
      PARAMETER (NMAX=500)
      real*8 nn,T,nns(NMAX),ffs(NMAX)
      real*8 rs,RsFromDensity,TinK
      real*8 AUToK
      real*8 eps,fp,fm
      real*8 nnALL,FALL,FnnALL,splint2
      real*8 fn1,fnn,f

      eps = 0.01d0
      fp  = 1.0d0+eps
      fm  = 1.0d0-eps

      rs = RsFromDensity(nn)
      TinK = AUToK(T)

      if (rs*fp.lt.0.8d0) then
         write(*,*) "Density is too large: rs=",rs
         stop
      endif

      if (rs*fm.gt.2.4d0) then
         write(*,*) "Density is too small: rs=",rs
         stop
      endif

      if (TinK*fp.lt.500d0) then
         write(*,*) "Tempature is too low",TinK
         stop
      endif

      if (TinK*fm.gt.1.28d08) then
         write(*,*) "Tempature is too high",TinK
         stop
      endif

      imin = 1
      imax = 5

      n=0
      do i=imin,imax
         n=n+1
c         nns(n) = nnALL(i)
         nns(n) = log(nnALL(i))
         ffs(n) = FALL(i,T)
c         write(*,*) i,n,nns(n),ffs(n)
      enddo
c      fn1 = FnnALL(imin,T)
c      fnn = FnnALL(imax,T)
      fn1 = FnnALL(imin,T)*nnALL(imin)
      fnn = FnnALL(imax,T)*nnALL(imax)
c      f = splint2(nns,ffs,n,fn1,fnn,nn)
      f = splint2(nns,ffs,n,fn1,fnn,log(nn))

c      if (nn*fp.lt.nns(1).or.nn*fm.gt.nns(n)) then
      if (log(nn*fp).lt.nns(1) .or. log(nn*fm).gt.nns(n)) then
         write(*,*) "Spline problem:",TinK,RsFromDensity(nns(1)),
     *        RsFromDensity(nn),
     *        RsFromDensity(nns(n)),imin,imax
      endif

c        Need to adjust constant "delta T = c*T"
c        This is a shift in the entropy term
c        Adjusted so that T=1.28*10^8 at rs=1.86 agree with Debye+Ideal fermi gas(el and He)
c      f = f - 5604.74529963d0*(TinK/1.28d8)
c      f = f - 5599.79d0*(TinK/1.28d8)
      f = f - 5583.17031902d0*(TinK/1.28d8)

      FreeEnergy = f
      return
      end

C****************************************************************************

      real*8 FUNCTION Pressure(nn,T)
      implicit none
      real*8 nn,T,dnn,nn1,nn2,F1,F2,dFdnn,P,eps
      real*8 FreeEnergy
      eps = 1.0d-6
      dnn = nn*eps
      nn1 = nn-dnn
      nn2 = nn+dnn
      F1 = FreeEnergy(nn1,T)
      F2 = FreeEnergy(nn2,T)
      dFdnn = (F2-F1)/(2.0d0*dnn)
      P = dFdnn*nn*nn
      Pressure = P
      return
      end

      real*8 FUNCTION Entropy(nn,T)
      implicit none
      real*8 nn,T,dT,T1,T2,F1,F2,dFdT,FreeEnergy,eps
      eps = 1.0d-6
      dT = T*eps
      T1 = T-dT
      T2 = T+dT
      F1 = FreeEnergy(nn,T1)
      F2 = FreeEnergy(nn,T2)
      dFdT = (F2-F1)/(2.0d0*dT)
      Entropy = -dFdT
      return
      end

      real*8 FUNCTION InternalEnergy(nn,T)
      implicit none
      real*8 nn,T,E,FreeEnergy,Entropy
      E = FreeEnergy(nn,T) + T * Entropy(nn,T)
      InternalEnergy = E
      return
      end
      
C****************************************************************************
      real*8 FUNCTION MinimumTemperatureForDensity(nn)
      implicit none
      real*8 rs,RsFromDensity,nn
      real*8 T,TinK,AUToK,KToAU
C      rs = RsFromDensity(nn)
      
      TinK = 500.0d0

      T = KToAU(TinK)
      MinimumTemperatureForDensity = T
      return
      end


      real*8 FUNCTION MaximumTemperatureForDensity(nn)
      implicit none
      real*8 rs,RsFromDensity,nn
      real*8 T,TinK,AUToK,KToAU
C      rs = RsFromDensity(nn)
      
      TinK = 1.28d08
      T = KToAU(TinK)
      MaximumTemperatureForDensity = T
      return
      end

      PROGRAM PrintEOS
      implicit none
      integer i,j,np
      real*8 nni,rs,rho
      real*8 T1,T2,T,TinK,AUToK,KToAU
      real*8 F,E,S,P
      real*8 PinGPa,EineV
      real*8 nnALL,RsFromDensity,DensityFromRs
      real*8 DensityToGCC,DensityFromGCC
      real*8 MinimumTemperatureForDensity,MaximumTemperatureForDensity
      real*8 FreeEnergy,Pressure,InternalEnergy,Entropy
      real*8 PressureToGPa,EnergyToeV

      write(*,*)                                                                     
      write(*,*) '+---------------------------------------------------+'
      write(*,*) '|                                                   |'
      write(*,*) '|     Equation of state of hot, dense HELIUM        |'
      write(*,*) '|                                                   |'
      write(*,*) '| Free energy interpolation code                    |'
      write(*,*) '| written by B. Militzer         Berkeley, 02-01-08 |'
      write(*,*) '|                                                   |'
      write(*,*) '+---------------------------------------------------+'

      do i=1,5
         nni = nnALL(i)
c         nni = DensityFromGCC(0.83201349337d0)
c         nni  = DensityFromRs(1.0d0)
c         nni  = DensityFromRs(1.86d0)
         rs  = RsFromDensity(nni)
         rho = DensityToGCC(nni)
c         write(*,*) nni,rs,rho
         T1  = MinimumTemperatureForDensity(nni)
         T2  = MaximumTemperatureForDensity(nni)
         np = 50
c         np = (T2-T1)/KToAU(500.0d0)+0.5d0
c         write(*,*) np
         write(*,*) 
         do j=0,np
            T = T1*exp(log(T2/T1)*real(j)/real(np))
c            T = T1+(T2-T1)*real(j)/real(np)
            TinK = AUToK(T)
C            write(*,*) "Trying",rs,TinK
            F = FreeEnergy(nni,T)
            P = Pressure(nni,T)
            E = InternalEnergy(nni,T)
            S = Entropy(nni,T)
            PinGPa = PressureToGPa(P)
            EineV  = EnergyToeV(E)
c            write(*,*) rs,TinK,PinGPa,E,EineV
 10         format(a,f8.6,a,f9.6,a,f13.3,a,f10.6,a,f14.6,
     *           a,f12.6,a,f10.6)
            write(*,10) 
     *           'rs= ',rs,
     *           '  rhoHe(g/cc)= ',rho,
     *           '  T(K)= ',TinK,
     *           '  E(Ha/e)= ',E,
     *           '  P(GPa)= ',PinGPa,
c     *           '  [F+cT](Ha/e)= ',F,
c     *           '  [S-c](Ha/kb/e)= ',S
     *           '  [F](Ha/e)= ',F,
     *           '  [S](Ha/kb/e)= ',S

         enddo
      enddo
      stop
      end
            
