** Equation of state for fully ionized electron-ion plasmas (EOS EIP)
* A.Y.Potekhin & G.Chabrier, Contrib. Plasma Phys., 50 (2010) 82, 
*       and references therein
* Please communicate comments/suggestions to Alexander Potekhin:
*                                            palex@astro.ioffe.ru
* Previously distributed versions (obsolete):
*   eos2000, eos2002, eos2004, eos2006, eos2007, eos2009, eos10,
*     and eos11.
* Last updated 22.12.12
**   L I S T   O F   S U B R O U T I N E S :
*  MAIN (normally commented-out) - example driving routine.
*  MELANGE9 - for arbitrary ionic mixture, renders total (ion+electron)
*          pressure, internal energy, entropy, heat capacity (all
*          normalized to the ionic ideal-gas values), logarithmic
*          derivatives of pressure over temperature and density.
*  EOSFI8 - nonideal (ion-ion + ion-electron + electron-electron)
*          contributions to the free and internal energies, pressure,
*          entropy, heat capacity, derivatives of pressure over
*          logarithm of temperature and over logarithm of density (all
*          normalized to the ionic ideal-gas values) for one ionic
*          component in a mixture.
*  FITION9 - ion-ion interaction contributions to the free and internal
*          energies, pressure, entropy, heat capacity, derivatives of
*          pressure over logarithms of temperature and density.
*  FSCRliq8 - ion-electron (screening) contributions to the free and
*          internal energies, pressure, entropy, heat capacity,
*          derivatives of pressure over logarithms of temperature and
*          density in the liquid phase for one ionic component in a
*          mixture.
*  FSCRsol8 - ion-electron (screening) contributions to the free and
*          internal energies, pressure, entropy, heat capacity,
*          derivatives of pressure over logarithms of temperature and
*          density for monoionic solid.
*  FHARM12 - harmonic (including static-lattice and zero-point)
*          contributions to the free and internal energies, pressure,
*          entropy, heat capacity, derivatives of pressure over
*          logarithms of temperature and density for solid OCP.
*  HLfit12 - the same as FHARM12, but only for thermal contributions
*  ANHARM8 - anharmonic contributions to the free and internal energies,
*          pressure, entropy, heat capacity, derivatives of pressure
*          over logarithms of temperature and density for solid OCP.
*  CORMIX - correction to the linear mixing rule for the Coulomb
*          contributions to the thermodynamic functions in the liquid.
*  ELECT11 - for an ideal electron gas of arbitrary degeneracy and
*          relativity at given temperature and electron chemical
*          potential, renders number density (in atomic units), free
*          energy, pressure, internal energy, entropy, heat capacity 
*          (normalized to the electron ideal-gas values), logarithmic
*          derivatives of pressure over temperature and density.
*  EXCOR7 - electron-electron (exchange-correlation) contributions to
*          the free and internal energies, pressure, entropy, heat
*          capacity, derivatives of pressure over logarithm of
*          temperature and over logarithm of density (all normalized
*          to the classical electron ideal-gas values).
*  FERINV7 - inverse non-relativistic Fermi integrals of orders -1/2,
*          1/2, 3/2, 5/2, and their first and second derivatives.
*  BLIN9 - relativistic Fermi-Dirac integrals of orders 1/2, 3/2, 5/2,
*          and their first, second, and some third derivatives.
*  CHEMFIT7 - electron chemical potential at given density and
*          temperature, and its first derivatives over density and
*          temperature and the second derivative over temperature.
**   I M P R O V E M E N T S   O F   2 0 0 8  -  2 0 0 9 :
*  FHARM8 uses a fit HLfit8 to the thermal free energy of the harmonic
*   Coulomb lattice, which is more accurate than its predecessor FHARM7.
*   Resulting corrections amount up to 20% for the ion heat capacity.
*   Accordingly, S/R D3fit and FthCHA7 deleted (not used anymore).
*  BLIN7 upgraded to BLIN8:
*      - cleaned (a never-reached if-else branch deleted);
*      - Sommerfeld (high-\chi) expansion improved;
*      - some third derivatives added.
*  CORMIX added (and MELANGE7 upgraded to MELANGE8 accordingly).
*  ANHARM7 upgraded to ANHARM8, more consistent with Wigner-Kirkwood.
*  Since the T- and rho-dependences of individual Z values in a mixture
*    are not considered, the corresponding inputs (AYLR, AYLT) are
*    excluded from MELANGE8 (and EOSFI7 changed to EOSFI8 accordingly).
*  ELECT7 upgraded to ELECT9 (high-degeneracy behaviour is improved)
**   I M P R O V E M E N T S   O F   2 0 1 0   -   2 0 1 1 :
*  ELECT9 upgraded (smooth match of two fits at chi >> 1)
*  BLIN8 replaced by BLIN9 - smooth fit interfaces at chi=0.6 and 14.
*  MELANGE8 replaced by MELANGE9 - slightly modified input/output
* 08.08.11 - corrected mistake (unlikely to have an effect) in CHEMFIT7
* 16.11.11 - ELECT9 upgraded to ELECT11 (additional output)
**   I M P R O V E M E N T S   O F   2 0 1 2 :
* 20.04.12 - FHARM8 and HLfit8 upgraded to FHARM12 and HLfit12:
*   output of HLfit12 does not include zero-point vibr., but provides U1
* 22.12.12 - MELANGE9 now includes a correction to the linear mixing
*   rule for the Madelung energy in the random bcc multi-ion lattice.
************************************************************************
*                           MAIN program:               Version 02.06.09
* This driving routine allows one to compile and run this code "as is".
* In practice, however, one usually needs to link subroutines from this
* file to another (external) code, therefore the MAIN program is
* normally commented-out.
C%C      implicit double precision (A-H), double precision (O-Z)
C%C      parameter(MAXY=10,UN_T6=.3157746)
C%C      dimension AY(MAXY),AZion(MAXY),ACMI(MAXY)
C%C      write(*,'('' Introduce the chemical composition (up to'',I3,
C%C     *  '' ion species):''/
C%C     *  '' charge number Z_i, atomic weight A_i,'',
C%C     *  '' partial number density x_i, derivatives d x_i / d ln T'',
C%C     *  '' and d x_i / d ln rho''/
C%C     /   '' (non-positive Z, A, or x=1 terminates the input)'')') MAXY
C%C      NMIX=0
C%C    3 continue
C%C      do IX=1,MAXY
C%C         write(*,'(''Z, A ('',I2,''): ''$)') IX
C%C         read*,AZion(IX),ACMI(IX)
C%C        if (AZion(IX).le.0..or.ACMI(IX).le.0.) goto 2
C%C         write(*,'(''x ('',I2,''): ''$)') IX
C%C         read*,AY(IX)
C%C        if (AY(IX).le.0..or.ACMI(IX).le.0.) goto 2
C%C         NMIX=IX
C%C        if (AY(IX).eq.1.d0) goto 2
C%C      enddo
C%C    2 continue
C%C      if (NMIX.eq.0) then
C%C         print*,'There must be at least one set of positive (x,Z,A).'
C%C        goto 3
C%C      endif
C%C      write(*,114)
C%C      do IX=1,NMIX
C%C         write(*,113) IX,AZion(IX),ACMI(IX),AY(IX)
C%C      enddo
C%C   10 continue
C%C      write(*,'('' Input lg(RHO) [g/cc], lg(T) (K) : ''$)')
C%C      read*,RHOlg,Tlg
C%C      T6=10.d0**(Tlg-6.d0)
C%C      RHO=10.d0**RHOlg
C%C      write(*,112)
C%C    1 continue
C%C      TEMP=T6/UN_T6 ! T [au]
C%C      call MELANGE9(NMIX,AY,AZion,ACMI,RHO,TEMP, ! input
C%C     *   PRADnkT, ! additional output - radiative pressure
C%C     *   DENS,Zmean,CMImean,Z2mean,GAMI,CHI,TPT,LIQSOL, ! output param.
C%C     *   PnkT,UNkT,SNk,CV,CHIR,CHIT) ! output dimensionless TD functions
C%C      Tnk=8.31447d13/CMImean*RHO*T6 ! n_i kT [erg/cc]
C%C      P=PnkT*Tnk/1.d12 ! P [Mbar]
C%C      TEGRAD=CHIT/(CHIT**2+CHIR*CV/PnkT) ! from Maxwell relat.
*   --------------------   OUTPUT   --------------------------------   *
* Here in the output we have:
* RHO - mass density in g/cc
* P - total pressure in Mbar (i.e. in 1.e12 dyn/cm^2)
* PnkT=P/nkT, where n is the number density of ions, T temperature
* CV - heat capacity at constant volume, divided by number of ions, /k
* CHIT - logarithmic derivative of pressure \chi_T
* CHIR - logarithmic derivative of pressure \chi_\rho
* UNkT - internal energy divided by NkT, N being the number of ions
* SNk - entropy divided by number of ions, /k
* GAMI - ionic Coulomb coupling parameter
* TPT=T_p/T, where T_p is the ion plasma temperature
* CHI - electron chemical potential, divided by kT
* LIQSOL = 0 in the liquid state, = 1 in the solid state
C%C      write(*,111) RHO,T6,P,PnkT,CV,CHIT,CHIR,UNkT,SNk,GAMI,TPT,CHI,
C%C     *  LIQSOL
C%C      goto 10
C%C  112 format(/
C%C     *  ' rho [g/cc]     T6 [K]      P [Mbar]   P/(n_i kT)  Cv/(N k)',
C%C     *  '     chi_T       chi_r      U/(N k T)    S/(N k)    Gamma_i',
C%C     *  '      T_p/T    chi_e liq/sol')
C%C  111 format(1P,12E12.3,I2)
C%C  113 format(I3,2F8.3,1PE12.4)
C%C  114 format('       Z      CMI        x_j')
C%C      end
      
      subroutine MELANGE9(NMIX,AY,AZion,ACMI,RHO,TEMP,PRADnkT,
     *   DENS,Zmean,CMImean,Z2mean,GAMImean,CHI,TPT,LIQSOL,
     *   PnkT,UNkT,SNk,CV,CHIR,CHIT)
*                                                       Version 22.12.12
* Stems from MELANGE8 v.26.12.09.
* Difference: output PRADnkT instead of input KRAD
* + EOS of fully ionized electron-ion plasma mixture.     
* Limitations:
* (a) inapplicable in the regimes of
*      (1) bound-state formation,
*      (2) quantum liquid,
*      (3) presence of positrons;
* (b) for the case of a composition gradually depending on RHO or TEMP,
*  second-order functions (CV,CHIR,CHIT in output) should not be trusted
* Choice of the liquid or solid regime - criterion GAMI [because the
*     choice based on comparison of total (non-OCP) free energies can be
*     sometimes dangerous because of the fit uncertainties ("Local field
*     correction" in solid and quantum effects in liquid are unknown)].
* Input: NMIX - number of different elements;
*        AY - their partial number densities,
*        AZion and ACMI - their charge and mass numbers,
*        RHO - total mass density [g/cc]
*        TEMP - temperature [in a.u.=2Ryd=3.1577e5 K].
* NB: instead of RHO, a true input is CHI, defined below
*     Hence, disagreement between RHO and DENS is the fit error (<0.4%)
* Output:
*         AY - rescaled so that to sum up to 1 and resorted (by AZion)
*         AZion - resorted in ascending order
*         ACMI - resorted in agreement with AZion
*         DENS - electron number density [in a.u.=6.7483346e24 cm^{-3}]
*         Zmean=<Z>, CMImean=<A> - mean ion charge and mass numbers,
*         Z2mean=<Z^2> - mean-square ion charge number
*         GAMImean - effective ion-ion Coulomb coupling constant
*         CHI = mu_e/kT, where mu_e is the electron chem.potential
*         TPT - effective ionic quantum parameter (T_p/T)
*         LIQSOL=0/1 for liquid/solid
*         SNk - dimensionless entropy per 1 ion
*         UNkT - internal energy per kT per ion
*         PnkT - pressure / n_i kT, where n_i is the ion number density
*         PRADnkT - radiative pressure / n_i kT
*         CV - heat capacity per ion, div. by Boltzmann const.
*         CHIR - inverse compressibility -(d ln P / d ln V)_T ("\chi_r")
*         CHIT = (d ln P / d ln T)_V ("\chi_T")
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter(TINY=1.d-7)
      dimension AY(*),AZion(*),ACMI(*)
      parameter (PI=3.141592653d0,C53=5.d0/3.d0,
     *   AUM=1822.888d0, ! a.m.u./m_e
     *   GAMIMELT=175., ! OCP value of Gamma_i for melting
     *   RSIMELT=140., ! ion density parameter of quantum melting
     *   RAD=2.554d-7) ! Radiation constant (=4\sigma/c) (in a.u.)
      Y=0.
      do IX=1,NMIX
         Y=Y+AY(IX)
      enddo
      if (dabs(Y-1.d0).gt.TINY) then
        do IX=1,NMIX
           AY(IX)=AY(IX)/Y
        enddo
         print*,'MELANGE9: partial densities (and derivatives)',
     *    ' are rescaled by factor',1./Y
      endif
* Sort the elements in ascending order in Z_j:
      KSORT=0
      do I=2,NMIX
         J=I
         Z=AZion(J)
         CMI=ACMI(J)
         Y=AY(J)
    1   if (J.le.1.or.AZion(J-1).le.Z) goto 2
           AZion(J)=AZion(J-1)
           ACMI(J)=ACMI(J-1)
           AY(J)=AY(J-1)
           J=J-1
           KSORT=1
        goto 1
    2    AZion(J)=Z
         ACMI(J)=CMI
         AY(J)=Y
      enddo
      if (KSORT.eq.1) write(*,'('' Ions are resorted as follows:''/
     *  '' i    Z_i       A_i       x_i''/(0P,I3,'':'',1P,3E10.3))')
     *  (J,AZion(J),ACMI(J),AY(J),J=1,NMIX)
* Calculation of average values:
      Zmean=0.
      Z2mean=0.
      Z52=0.
      Z53=0.
      Z73=0.
      Z321=0. ! corr.26.12.09
      CMImean=0.
      do IX=1,NMIX
         Zmean=Zmean+AY(IX)*AZion(IX)
         Z2mean=Z2mean+AY(IX)*AZion(IX)**2
         Z13=AZion(IX)**(1./3.)
         Z53=Z53+AY(IX)*Z13**5 
         Z73=Z73+AY(IX)*Z13**7
         Z52=Z52+AY(IX)*dsqrt(AZion(IX))**5
         Z321=Z321+AY(IX)*AZion(IX)*dsqrt(AZion(IX)+1.d0)**3 ! 26.12.09
         CMImean=CMImean+AY(IX)*ACMI(IX)
      enddo
* (0) Photons:
      UINTRAD=RAD*TEMP**4
      PRESSRAD=UINTRAD/3.
C      CVRAD=4.*UINTRAD/TEMP
* (1) ideal electron gas (including relativity and degeneracy)  -----  *
      DENS=RHO/11.20587*Zmean/CMImean ! number density of electrons [au]
      call CHEMFIT(DENS,TEMP,CHI)
* NB: CHI can be used as true input instead of RHO or DENS
      call ELECT11(TEMP,CHI,
     *  DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,
     *  DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
* NB: at this point DENS is redefined (the difference can be ~0.1%)
      DTE=DENS*TEMP
      PRESSE=PEid*DTE ! P_e [a.u.]
      UINTE=UEid*DTE ! U_e / V [a.u.]
* (2) non-ideal Coulomb EIP  ----------------------------------------  *
      RS=(.75/PI/DENS)**.333333 ! r_s - electron density parameter
      RSI=RS*CMImean*Z73*AUM ! R_S - ion density parameter
      GAME=1./RS/TEMP ! electron Coulomb parameter Gamma_e
      GAMImean=Z53*GAME   ! effective Gamma_i - ion Coulomb parameter
      if (GAMImean.lt.GAMIMELT.or.RSI.lt.RSIMELT) then
         LIQSOL=0 ! liquid regime
      else
         LIQSOL=1 ! solid regime
      endif
      TPT2=0.
      do IX=1,NMIX
         TPT2=TPT2+AY(IX)*GAME**2/RS*3./(AUM*ACMI(IX))*AZion(IX)
* In the case of a mixture, this estimate of the quantum ion parameter
* is as crude as the linear mixing model for Fharm employed below.
      enddo
      if (LIQSOL.eq.0.and.TPT2.gt.18.) then
         print*,'MELANGE9: strong quantum effects in liquid!'
         read(*,'(A)')
      endif
      TPT=dsqrt(TPT2) ! effective T_p/T - ion quantum parameter
* Calculate partial thermodynamic quantities and combine them together:
      UINT=UINTE
      PRESS=PRESSE
      CVtot=CVE*DENS
      Stot=SEid*DENS
      PDLT=PRESSE*CHITE ! d P_e[a.u.] / d ln T
      PDLR=PRESSE*CHIRE ! d P_e[a.u.] / d ln\rho
      DENSI=DENS/Zmean ! number density of all ions
      PRESSI=DENSI*TEMP ! ideal-ions total pressure (normalization)
* Add Coulomb+xc nonideal contributions, and ideal free energy:
      do IX=1,NMIX
        if (AY(IX).lt.TINY) goto 10 ! skip this species
         Zion=AZion(IX)
         CMI=ACMI(IX)
         GAMI=Zion**C53/RS/TEMP ! Gamma_i for given ion species
         DNI=DENSI*AY(IX) ! number density of ions of given type
         PRI=DNI*TEMP ! = ideal-ions partial pressure (normalization)
         call EOSFI8(LIQSOL,CMI,Zion,RS,GAMI,
     *     FC1,UC1,PC1,SC1,CV1,PDT1,PDR1,
     *     FC2,UC2,PC2,SC2,CV2,PDT2,PDR2)
* First-order TD functions:
         UINT=UINT+UC2*PRI ! internal energy density (e+i+Coul.)
         Stot=Stot+DNI*(SC2-dlog(AY(IX))) !entropy per unit volume[a.u.]
         PRESS=PRESS+PC2*PRI ! pressure (e+i+Coul.) [a.u.]
* Second-order functions (they take into account compositional changes):
         CVtot=CVtot+DENSI*CV2*AY(IX) ! C_V (e+i+Coul.)/ V
         PDLT=PDLT+PRI*PDT2 ! d P / d ln T
         PDLR=PDLR+PRI*PDR2 ! d P / d ln\rho
   10   continue
      enddo ! next IX
* Corrections to the linear mixing rule:
      if (LIQSOL.eq.0) then ! liquid phase
         call CORMIX(RS,GAME,Zmean,Z2mean,Z52,Z53,Z321,
     *     FMIX,UMIX,PMIX,CVMIX,PDTMIX,PDRMIX)
      else ! solid phase (only Madelung contribution) [22.12.12]
         FMIX=0.
        do I=1,NMIX
          do J=I+1,NMIX
             RZ=AZion(J)/AZion(I)
             X2=AY(J)/(AY(I)+AY(J))
             X1=dim(1.d0,X2)
            if (X2.lt.TINY) goto 11
             X=X2/RZ+(1.d0-1.d0/RZ)*X2**RZ
             GAMI=Zion**C53/RS/TEMP ! Gamma_i for given ion species
             DeltaG=.012*(1.d0-1.d0/RZ**2)*(X1+X2*RZ**C53)
             DeltaG=DeltaG*X/X2*dim(1.d0,X)/X1
             FMIX=FMIX+AY(I)*AY(J)*GAMI*DeltaG
   11       continue
          enddo
   12     continue
        enddo
         UMIX=FMIX
         PMIX=FMIX/3.d0
         CVMIX=0.
         PDTMIX=0.
         PDRMIX=FMIX/2.25d0
      endif
      UINT=UINT+UMIX*PRESSI
      Stot=Stot+DENSI*(UMIX-FMIX)
      PRESS=PRESS+PMIX*PRESSI
      CVtot=CVtot+DENSI*CVMIX
      PDLT=PDLT+PRESSI*PDTMIX
      PDLR=PDLR+PRESSI*PDRMIX
* First-order:
      PRADnkT=PRESSRAD/PRESSI ! radiative pressure / n_i k T
C      CVtot=CVtot+CVRAD
C      Stot=Stot+CVRAD/3.
      PnkT=PRESS/PRESSI ! P / n_i k T
      UNkT=UINT/PRESSI ! U / N_i k T
C      UNkT=UNkT+UINTRAD/PRESSI
      SNk=Stot/DENSI ! S / N_i k
* Second-order:
      CV=CVtot/DENSI ! C_V per ion
      CHIR=PDLR/PRESS ! d ln P / d ln\rho
      CHIT=PDLT/PRESS ! d ln P / d ln T
C      CHIT=CHIT+4.*PRESSRAD/PRESS ! d ln P / d ln T
      return
      end

      subroutine EOSFI8(LIQSOL,CMI,Zion,RS,GAMI,
     *  FC1,UC1,PC1,SC1,CV1,PDT1,PDR1,
     *  FC2,UC2,PC2,SC2,CV2,PDT2,PDR2)
*                                                       Version 16.09.08
*                 call FHARM8 has been replaced by call FHARM12 27.04.12
* Non-ideal parts of thermodynamic functions in the fully ionized plasma
* Stems from EOSFI5 and EOSFI05 v.04.10.05
* Input: LIQSOL=0/1(liquid/solid), 
*        Zion,CMI - ion charge and mass numbers,
*        RS=r_s (electronic density parameter),
*        GAMI=Gamma_i (ion coupling),
* Output: FC1 and UC1 - non-ideal "ii+ie+ee" contribution to the 
*         free and internal energies (per ion per kT),
*         PC1 - analogous contribution to pressure divided by (n_i kT),
*         CV1 - "ii+ie+ee" heat capacity per ion [units of k]
*         PDT1=(1/n_i kT)*(d P_C/d ln T)_V
*         PDR1=(1/n_i kT)*(d P_C/d ln\rho)_T
* FC2,UC2,PC2,SC2,CV2 - analogous to FC1,UC1,PC1,SC1,CV1, but including
* the part corresponding to the ideal ion gas. This is useful for 
* preventing accuracy loss in some cases (e.g., when SC2 << SC1).
* FC2 does not take into account the entropy of mixing S_{mix}: in a
* mixture, S_{mix}/(N_i k) has to be added externally (see MELANGE9).
* FC2 does not take into account the ion spin degeneracy either.
* When needed, the spin term must be added to the entropy externally.
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter(TINY=1.d-20) ! prevents division by zero
      parameter (AUM=1822.888d0) ! a.m.u/m_e
      if (LIQSOL.ne.1.and.LIQSOL.ne.0) stop'EOSFI8: invalid LIQSOL'
      if (CMI.le..1) stop'EOSFI8: too small CMI'
      if (Zion.le..1) stop'EOSFI8: too small Zion'
      if (RS.le..0) stop'EOSFI8: invalid RS'
      if (GAMI.le..0) stop'EOSFI8: invalid GAMI'
      GAME=GAMI/Zion**1.666667
      call EXCOR7(RS,GAME,FXC,UXC,PXC,CVXC,SXC,PDTXC,PDRXC) ! "ee"("xc")
* Calculate "ii" part:
      COTPT=dsqrt(3./AUM/CMI)/Zion**(7./6.) ! auxiliary coefficient
      TPT=GAMI/dsqrt(RS)*COTPT              ! T_p/T
      FidION=1.5*dlog(TPT**2/GAMI)-1.323515
* 1.3235=1+0.5*ln(6/pi); FidION = F_{id.ion gas}/(N_i kT), but without
* the term x_i ln x_i = -S_{mix}/(N_i k).
      if (LIQSOL.eq.0) then                 ! liquid
         call FITION9(GAMI,
     *     FION,UION,PION,CVii,PDTii,PDRii)
         FItot=FION+FidION
         UItot=UION+1.5
         PItot=PION+1.d0
         CVItot=CVii+1.5d0
         SCItot=UItot-FItot
         PDTi=PDTii+1.d0
         PDRi=PDRii+1.d0
         FWK=TPT**2/24.d0 ! Wigner-Kirkwood (quantum diffr.) term
         UWK=2.d0*FWK
         CVWK=-UWK
         PWK=FWK
         PDRWK=2.d0*PWK
         PDTWK=-PWK
      else                                  ! solid
         call FHARM12(GAMI,TPT,
     *     Fharm,Uharm,Pharm,CVharm,Sharm,PDTharm,PDRharm) ! harm."ii"
         call ANHARM8(GAMI,TPT,Fah,Uah,Pah,CVah,PDTah,PDRah) ! anharm.
         FItot=Fharm+Fah
         FION=FItot-FidION
         UItot=Uharm+Uah
         UION=UItot-1.5d0 ! minus 1.5=ideal-gas, in order to get "ii"
         PItot=Pharm+Pah
         PION=PItot-1.d0 ! minus 1=ideal-gas
         PDTi=PDTharm+PDTah
         PDRi=PDRharm+PDRah
         PDTii=PDTi-1.d0 ! minus 1=ideal-gas
         PDRii=PDRi-1.d0 ! minus 1=ideal-gas
         CVItot=CVharm+CVah
         SCItot=Sharm+Uah-Fah
         CVii=CVItot-1.5d0 ! minus 1.5=ideal-gas
         FWK=0.
         UWK=0.
         PWK=0.
         CVWK=0.
         PDRWK=0.
         PDTWK=0.
      endif
* Calculate "ie" part:
      if (LIQSOL.eq.1) then
         call FSCRsol8(RS,GAMI,Zion,TPT,
     *     FSCR,USCR,PSCR,S_SCR,CVSCR,PDTSCR,PDRSCR)
      else
         call FSCRliq8(RS,GAME,Zion,
     *     FSCR,USCR,PSCR,CVSCR,PDTSCR,PDRSCR)
         S_SCR=USCR-FSCR
      endif
* Total excess quantities ("ii"+"ie"+"ee", per ion):
      FC0=FSCR+Zion*FXC+FWK
      UC0=USCR+Zion*UXC+UWK
      PC0=PSCR+Zion*PXC+PWK
      SC0=S_SCR+Zion*SXC+FWK
      CV0=CVSCR+Zion*CVXC+CVWK
      PDT0=PDTSCR+Zion*PDTXC+PDTWK
      PDR0=PDRSCR+Zion*PDRXC+PDRWK
      FC1=FION+FC0
      UC1=UION+UC0
      PC1=PION+PC0
      SC1=(UION-FION)+SC0
      CV1=CVii+CV0
      PDT1=PDTii+PDT0
      PDR1=PDRii+PDR0
* Total excess + ideal-ion quantities
      FC2=FItot+FC0
      UC2=UItot+UC0
      PC2=PItot+PC0
      SC2=SCItot+SC0
      CV2=CVItot+CV0
      PDT2=PDTi+PDT0
      PDR2=PDRi+PDR0
      return
      end

* ==================  ELECTRON-ION COULOMB LIQUID  =================== *
      subroutine FITION9(GAMI,FION,UION,PION,CVii,PDTii,PDRii)
*                                                       Version 11.09.08
* Dummy argument Zion is deleted in 2009.
* Non-ideal contributions to thermodynamic functions of classical OCP.
*   Stems from FITION00 v.24.05.00.
* Input: GAMI - ion coupling parameter
* Output: FION - ii free energy / N_i kT
*         UION - ii internal energy / N_i kT
*         PION - ii pressure / n_i kT
*         CVii - ii heat capacity / N_i k
*         PDTii = PION + d(PION)/d ln T = (1/N_i kT)*(d P_{ii}/d ln T)
*         PDRii = PION + d(PION)/d ln\rho
*   Parameters adjusted to Caillol (1999).
      implicit double precision (A-H),double precision (O-Z)
      save
      parameter (A1=-.907347d0,A2=.62849d0,C1=.004500d0,G1=170.0,
     *  C2=-8.4d-5,G2=.0037,SQ32=.8660254038d0) ! SQ32=sqrt(3)/2
      A3=-SQ32-A1/dsqrt(A2)
      F0=A1*(dsqrt(GAMI*(A2+GAMI))-
     -     A2*dlog(dsqrt(GAMI/A2)+dsqrt(1.+GAMI/A2)))+
     +     2.*A3*(dsqrt(GAMI)-datan(dsqrt(GAMI)))
      U0=dsqrt(GAMI)**3*(A1/dsqrt(A2+GAMI)+A3/(1.d0+GAMI))
*   This is the zeroth approximation. Correction:
      UION=U0+C1*GAMI**2/(G1+GAMI)+C2*GAMI**2/(G2+GAMI**2)
      FION=F0+C1*(GAMI-G1*dlog(1.d0+GAMI/G1))+
     +   C2/2.*dlog(1.d0+GAMI**2/G2)
      CVii=-0.5*dsqrt(GAMI)**3*(A1*A2/dsqrt(A2+GAMI)**3+
     +  A3*(1.d0-GAMI)/(1.d0+GAMI)**2) -
     -  GAMI**2*(C1*G1/(G1+GAMI)**2+C2*(G2-GAMI**2)/(G2+GAMI**2)**2)
      PION=UION/3.
      PDRii=(4.*UION-CVii)/9. ! p_{ii} + d p_{ii} / d ln\rho
      PDTii=CVii/3. ! p_{ii} + d p_{ii} / d ln T
      return
      end

      subroutine FSCRliq8(RS,GAME,Zion,
     *     FSCR,USCR,PSCR,CVSCR,PDTSCR,PDRSCR) ! fit to the el.-ion scr.
*                                                       Version 11.09.08
*                                                       cleaned 16.06.09
* Stems from FSCRliq7 v. 09.06.07. Included a check for RS=0.
*   INPUT: RS - density parameter, GAME - electron Coulomb parameter,
*          Zion - ion charge number,
*   OUTPUT: FSCR - screening (e-i) free energy per kT per 1 ion,
*           USCR - internal energy per kT per 1 ion (screen.contrib.)
*           PSCR - pressure divided by (n_i kT) (screen.contrib.)
*           CVSCR - heat capacity per 1 ion (screen.contrib.)
*           PDTSCR,PDRSCR = PSCR + d PSCR / d ln(T,\rho)
      implicit double precision(A-H),double precision(O-Z)
      save
      parameter(XRS=.0140047,TINY=1.d-19)
      if (RS.lt.0.) stop'FSCRliq8: RS < 0'
      if (RS.lt.TINY) then
         FSCR=0.
         USCR=0.
         PSCR=0.
         CVSCR=0.
         PDTSCR=0.
         PDRSCR=0.
         return
      endif
      SQG=sqrt(GAME)
      SQR=sqrt(RS)
      SQZ1=dsqrt(1.+Zion)
      SQZ=dsqrt(Zion)
      CDH0=Zion/1.73205 ! 1.73205=sqrt(3.)
      CDH=CDH0*(SQZ1**3-SQZ**3-1.)
      SQG=sqrt(GAME)
      ZLN=dlog(Zion)
      Z13=exp(ZLN/3.) ! Zion**(1./3.)
      X=XRS/RS ! relativity parameter
      CTF=Zion**2*.2513*(Z13-1.+.2/sqrt(Z13))
* Thomas-Fermi constant; .2513=(18/175)(12/\pi)^{2/3}
      P01=1.11*exp(.475*ZLN)
      P03=0.2+0.078*ZLN**2
      PTX=1.16+.08*ZLN
      TX=GAME**PTX
      TXDG=PTX*TX/GAME
      TXDGG=(PTX-1.)*TXDG/GAME
      TY1=1./(1.d-3*Zion**2+2.*GAME)
      TY1DG=-2.*TY1**2
      TY1DGG=-4.*TY1*TY1DG
      TY2=1.+6.*RS**2
      TY2DX=-12.*RS**2/X
      TY2DXX=-3.*TY2DX/X
      TY=RS**3/TY2*(1.+TY1)
      TYX=3./X+TY2DX/TY2
      TYDX=-TY*TYX
      TYDG=RS**3*TY1DG/TY2
      P1=(Zion-1.)/9.
      COR1=1.+P1*TY
      COR1DX=P1*TYDX
      COR1DG=P1*TYDG
      COR1DXX=P1*(TY*(3./X**2+(TY2DX/TY2)**2-TY2DXX/TY2)-TYDX*TYX)
      COR1DGG=P1*RS**3*TY1DGG/TY2
      COR1DXG=-P1*TYDG*TYX
      U0=.78*sqrt(GAME/Zion)*RS**3
      U0DX=-3.*U0/X
      U0DG=.5*U0/GAME
      U0DXX=-4.*U0DX/X
      U0DGG=-.5*U0DG/GAME
      U0DXG=-3.*U0DG/X
      D0DG=Zion**3
      D0=GAME*D0DG+21.*RS**3
      D0DX=-63.*RS**3/X
      D0DXX=252.*RS**3/X**2
      COR0=1.+U0/D0
      COR0DX=(U0DX-U0*D0DX/D0)/D0
      COR0DG=(U0DG-U0*D0DG/D0)/D0
      COR0DXX=(U0DXX-(2.*U0DX*D0DX+U0*D0DXX)/D0+2.*(D0DX/D0)**2)/D0
      COR0DGG=(U0DGG-2.*U0DG*D0DG/D0+2.*U0*(D0DG/D0)**2)/D0
      COR0DXG=(U0DXG-(U0DX*D0DG+U0DG*D0DX)/D0+2.*U0*D0DX*D0DG/D0**2)/D0
* Relativism:
      RELE=dsqrt(1.d0+X**2)
      Q1=.18/dsqrt(dsqrt(Zion))
      Q2=.2+.37/dsqrt(Zion)
      H1U=1.+X**2/5.
      H1D=1.+Q1*X+Q2*X**2
      H1=H1U/H1D
      H1X=.4*X/H1U-(Q1+2.*Q2*X)/H1D
      H1DX=H1*H1X
      H1DXX=H1DX*H1X+
     +  H1*(.4/H1U-(.4*X/H1U)**2-2.*Q2/H1D+((Q1+2.*Q2*X)/H1D)**2)
      UP=CDH*SQG+P01*CTF*TX*COR0*H1
      UPDX=P01*CTF*TX*(COR0DX*H1+COR0*H1DX)
      UPDG=.5*CDH/SQG+P01*CTF*(TXDG*COR0+TX*COR0DG)*H1
      UPDXX=P01*CTF*TX*(COR0DXX*H1+2.*COR0DX*H1DX+COR0*H1DXX)
      UPDGG=-.25*CDH/(SQG*GAME)+
     +  P01*CTF*(TXDGG*COR0+2.*TXDG*COR0DG+TX*COR0DGG)*H1
      UPDXG=P01*CTF*(TXDG*(COR0DX*H1+COR0*H1DX)+
     +  TX*(COR0DXG*H1+COR0DG*H1DX))
      DN1=P03*SQG+P01/RS*TX*COR1
      DN1DX=P01*TX*(COR1/XRS+COR1DX/RS)
      DN1DG=.5*P03/SQG+P01/RS*(TXDG*COR1+TX*COR1DG)
      DN1DXX=P01*TX/XRS*(2.*COR1DX+X*COR1DXX)
      DN1DGG=-.25*P03/(GAME*SQG)+
     +  P01/RS*(TXDGG*COR1+2.*TXDG*COR1DG+TX*COR1DGG)
      DN1DXG=P01*(TXDG*(COR1/XRS+COR1DX/RS)+TX*(COR1DG/XRS+COR1DXG/RS))
      DN=1.+DN1/RELE
      DNDX=DN1DX/RELE-X*DN1/RELE**3
      DNDXX=(DN1DXX-((2.*X*DN1DX+DN1)-3.*X**2*DN1/RELE**2)/RELE**2)/RELE
      DNDG=DN1DG/RELE
      DNDGG=DN1DGG/RELE
      DNDXG=DN1DXG/RELE-X*DN1DG/RELE**3
      FSCR=-UP/DN*GAME
      FX=(UP*DNDX/DN-UPDX)/DN
      FXDG=((UPDG*DNDX+UPDX*DNDG+UP*DNDXG-2.*UP*DNDX*DNDG/DN)/DN-
     -  UPDXG)/DN
      FDX=FX*GAME
      FG=(UP*DNDG/DN-UPDG)/DN
      FDG=FG*GAME-UP/DN
      FDGDH=SQG*DNDG/DN**2 ! d FDG / d CDH
      FDXX=((UP*DNDXX+2.*(UPDX*DNDX-UP*DNDX**2/DN))/DN-UPDXX)/DN*GAME
      FDGG=2.*FG+GAME*((2.*DNDG*(UPDG-UP*DNDG/DN)+UP*DNDGG)/DN-UPDGG)/DN
      FDXG=FX+GAME*FXDG
      USCR=GAME*FDG
      CVSCR=-GAME**2*FDGG
      PSCR=(X*FDX+GAME*FDG)/3.
      PDTSCR=-GAME**2*(X*FXDG+FDGG)/3.
      PDRSCR=(12.*PSCR+X**2*FDXX+2.*X*GAME*FDXG+GAME**2*FDGG)/9.
      return
      end

* ==============   SUBROUTINES FOR THE SOLID STATE   ================= *
      subroutine FSCRsol8(RS,GAMI,Zion,TPT,
     *     FSCR,USCR,PSCR,S_SCR,CVSCR,PDTSCR,PDRSCR)
*                                                       Version 28.05.08
*                    undefined zero variable Q1DXG is wiped out 21.06.10
* Fit to the el.-ion screening in bcc or fcc Coulomb solid
* Stems from FSCRsol8 v.09.06.07. Included a check for RS=0.
*   INPUT: RS - el. density parameter, GAMI - ion coupling parameter,
*          Zion - ion charge, TPT=T_p/T - ion quantum parameter
*   OUTPUT: FSCR - screening (e-i) free energy per kT per 1 ion,
*           USCR - internal energy per kT per 1 ion (screen.contrib.)
*           PSCR - pressure divided by (n_i kT) (screen.contrib.)
*           S_SCR - screening entropy contribution / (N_i k)
*           CVSCR - heat capacity per 1 ion (screen.contrib.)
*           PDTSCR,PDRSCR = PSCR + d PSCR / d ln(T,\rho)
      implicit double precision(A-H),double precision(O-Z)
      save
      dimension AP(4) ! parameters of the fit
      parameter (ENAT=2.7182818285d0,TINY=1.d-19)
      data AP/1.1866,.684,17.9,41.5/,PX/.205/ ! for bcc lattice
cc      data AP/1.1857,.663,17.1,40./,PX/.212/ ! for fcc lattice
      if (RS.lt.0.) stop'FSCRliq8: RS < 0'
      if (RS.lt.TINY) then
         FSCR=0.
         USCR=0.
         PSCR=0.
         S_SCR=0.
         CVSCR=0.
         PDTSCR=0.
         PDRSCR=0.
         return
      endif
      XSR=.0140047/RS ! relativity parameter
      Z13=Zion**.3333333
      P1=.00352*(1.-AP(1)/Zion**.267+.27/Zion)
      P2=1.+2.25/Z13*
     *(1.+AP(2)*Zion**5+.222*Zion**6)/(1.+.222*Zion**6)
      ZLN=dlog(Zion)
      Finf=sqrt(P2/XSR**2+1.)*Z13**2*P1 ! The TF limit
      FinfX=-P2/((P2+XSR**2)*XSR)
      FinfDX=Finf*FinfX
      FinfDXX=FinfDX*FinfX-FinfDX*(P2+3.*XSR**2)/((P2+XSR**2)*XSR)
      R1=AP(4)/(1.+ZLN)
      R2=.395*ZLN+.347/Zion/sqrt(Zion)
      R3=1./(1.+ZLN*sqrt(ZLN)*.01+.097/Zion**2)
      Q1U=R1+AP(3)*XSR**2
      Q1D=1.+R2*XSR**2
      Q1=Q1U/Q1D
      Q1X=2.*XSR*(AP(3)/Q1U-R2/Q1D)
      Q1XDX=Q1X/XSR+4.*XSR**2*((R2/Q1D)**2-(AP(3)/Q1U)**2)
      Q1DX=Q1*Q1X
      Q1DXX=Q1DX*Q1X+Q1*Q1XDX
* New quantum factor, in order to suppress CVSCR at TPT >> 1
      if (TPT.lt.6./PX) then
         Y0=(PX*TPT)**2
         Y0DX=Y0/XSR
         Y0DG=2.*Y0/GAMI
         Y0DXX=0.
         Y0DGG=Y0DG/GAMI
         Y0DXG=Y0DG/XSR
         Y1=dexp(Y0)
         Y1DX=Y1*Y0DX
         Y1DG=Y1*Y0DG
         Y1DXX=Y1*(Y0DX**2+Y0DXX)
         Y1DGG=Y1*(Y0DG**2+Y0DGG)
         Y1DXG=Y1*(Y0DX*Y0DG+Y0DXG)
         SA=1.d0+Y1
         SUPA=dlog(SA)
         SUPADX=Y1DX/SA
         SUPADG=Y1DG/SA
         SUPADXX=(Y1DXX-Y1DX**2/SA)/SA
         SUPADGG=(Y1DGG-Y1DG**2/SA)/SA
         SUPADXG=(Y1DXG-Y1DX*Y1DG/SA)/SA
         EM2=ENAT-2.d0
         SB=ENAT-EM2/Y1
         SUPB=dlog(SB)
         EM2Y1=EM2/(Y1**2*SB)
         SUPBDX=EM2Y1*Y1DX
         SUPBDG=EM2Y1*Y1DG
         SUPBDXX=EM2Y1*(Y1DXX-2.*Y1DX**2/Y1-Y1DX*SUPBDX)
         SUPBDGG=EM2Y1*(Y1DGG-2.*Y1DG**2/Y1-Y1DG*SUPBDG)
         SUPBDXG=EM2Y1*(Y1DXG-2.*Y1DX*Y1DG/Y1-Y1DG*SUPBDX)
         SUP=dsqrt(SUPA/SUPB)
         SUPX=.5*(SUPADX/SUPA-SUPBDX/SUPB)
         SUPDX=SUP*SUPX
         SUPG=.5*(SUPADG/SUPA-SUPBDG/SUPB)
         SUPDG=SUP*SUPG
         SUPDXX=SUPDX*SUPX+
     +     SUP*.5*(SUPADXX/SUPA-(SUPADX/SUPA)**2-
     -             SUPBDXX/SUPB+(SUPBDX/SUPB)**2)
         SUPDGG=SUPDG*SUPG+
     +     SUP*.5*(SUPADGG/SUPA-(SUPADG/SUPA)**2-
     -             SUPBDGG/SUPB+(SUPBDG/SUPB)**2)
         SUPDXG=SUPDX*SUPG+
     +     SUP*.5*((SUPADXG-SUPADX*SUPADG/SUPA)/SUPA-
     -             (SUPBDXG-SUPBDX*SUPBDG/SUPB)/SUPB)
      else
         SUP=PX*TPT
         SUPDX=.5*PX*TPT/XSR
         SUPDG=PX*TPT/GAMI
         SUPDXX=-.5*SUPDX/XSR
         SUPDGG=0.
         SUPDXG=SUPDX/GAMI
      endif
      GR3=(GAMI/SUP)**R3
      GR3X=-R3*SUPDX/SUP
      GR3DX=GR3*GR3X
      GR3DXX=GR3DX*GR3X-R3*GR3*(SUPDXX/SUP-(SUPDX/SUP)**2)
      GR3G=R3*(1./GAMI-SUPDG/SUP)
      GR3DG=GR3*GR3G
      GR3DGG=GR3DG*GR3G+GR3*R3*((SUPDG/SUP)**2-SUPDGG/SUP-1./GAMI**2)
      GR3DXG=GR3DG*GR3X+GR3*R3*(SUPDX*SUPDG/SUP**2-SUPDXG/SUP)
      W=1.+Q1/GR3
      WDX=Q1DX/GR3-Q1*GR3DX/GR3**2
      WDG=-Q1*GR3DG/GR3**2
      WDXX=Q1DXX/GR3-(2.*Q1DX*GR3DX+Q1*(GR3DXX-2.*GR3DX**2/GR3))/GR3**2
      WDGG=Q1*(2.*GR3DG**2/GR3-GR3DGG)/GR3**2
      WDXG=-(Q1DX*GR3DG+Q1*(GR3DXG-2.*GR3DX*GR3DG/GR3))/GR3**2
      FSCR=-GAMI*Finf*W
      FDX=-GAMI*(FinfDX*W+Finf*WDX)
      FDXX=-GAMI*(FinfDXX*W+2.*FinfDX*WDX+Finf*WDXX)
      FDG=-Finf*W-GAMI*Finf*WDG
      FDGG=-2.*Finf*WDG-GAMI*Finf*WDGG
      FDXG=-FinfDX*W-Finf*WDX-GAMI*(FinfDX*WDG+Finf*WDXG)
      S_SCR=-GAMI**2*Finf*WDG
      USCR=S_SCR+FSCR
      CVSCR=-GAMI**2*FDGG
      PSCR=(XSR*FDX+GAMI*FDG)/3.
      PDTSCR=GAMI**2*(XSR*Finf*(FinfX*WDG+WDXG)-FDGG)/3.
      PDRSCR=(12.*PSCR+XSR**2*FDXX+2.*XSR*GAMI*FDXG+
     +  GAMI**2*FDGG)/9.
      return
      end

      subroutine ANHARM8(GAMI,TPT,Fah,Uah,Pah,CVah,PDTah,PDRah)
* ANHARMONIC free energy                                Version 27.07.07
*                                                       cleaned 16.06.09
* Stems from ANHARM8b. Difference: AC=0., B1=.12 (.1217 - over accuracy)
* Input: GAMI - ionic Gamma, TPT=Tp/T - ionic quantum parameter
* Output: anharm.free en. Fah=F_{AH}/(N_i kT), internal energy Uah,
*   pressure Pah=P_{AH}/(n_i kT), specific heat CVah = C_{V,AH}/(N_i k),
*   PDTah = Pah + d Pah / d ln T, PDRah = Pah + d Pah / d ln\rho
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter(NM=3)
      dimension AA(NM)
      data AA/10.9,247.,1.765d5/ ! Farouki & Hamaguchi'93
      data B1/.12/ ! coeff.at \eta^2/\Gamma at T=0
      CK=B1/AA(1) ! fit coefficient
      TPT2=TPT**2
      TPT4=TPT2**2
      TQ=B1*TPT2/GAMI ! quantum dependence
      TK2=CK*TPT2
      SUP=dexp(-TK2) ! suppress.factor of class.anharmonicity
      Fah=0.
      Uah=0.
      Pah=0.
      CVah=0.
      PDTah=0.
      PDRah=0.
      SUPGN=SUP
      do N=1,NM
         CN=N
         SUPGN=SUPGN/GAMI ! SUP/Gamma^n
         ACN=AA(N)
         Fah=Fah-ACN/CN*SUPGN
         Uah=Uah+(ACN*(1.+2.*TK2/CN))*SUPGN
         PN=AA(N)/3.+TK2*AA(N)/CN
         Pah=Pah+PN*SUPGN
         CVah=CVah+((CN+1.)*AA(N)+(4.-2./CN)*AA(N)*TK2+
     +     4.*AA(N)*CK**2/CN*TPT4)*SUPGN
         PDTah=PDTah+(PN*(1.+CN+2.*TK2)-2./CN*AA(N)*TK2)*SUPGN
         PDRah=PDRah+(PN*(1.-CN/3.-TK2)+AA(N)/CN*TK2)*SUPGN
      enddo
      Fah=Fah-TQ
      Uah=Uah-TQ
      Pah=Pah-TQ/1.5
      PDRah=PDRah-TQ/4.5
      return
      end

      subroutine FHARM12(GAMI,TPT,
     *   Fharm,Uharm,Pharm,CVth,Sth,PDTharm,PDRharm)
* Thermodynamic functions of a harmonic crystal, incl.stat.Coul.lattice
* 
*                                                       Version 27.04.12
* Stems from FHARM8 v.15.02.08
* Replaced HLfit8 with HLfit12: rearranged output.
* Input: GAMI - ionic Gamma, TPT=T_{p,i}/T
* Output: Fharm=F/(N_i T), Uharm=U/(N_i T), Pharm=P/(n_i T),
* CVth=C_V/N_i, Sharm=S/N_i
* PDTharm = Pharm + d Pharm / d ln T, PDRharm = Pharm + d Pharm/d ln\rho
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter(CM=.895929256d0) ! Madelung
      call HLfit12(TPT,F,U,CVth,Sth,U1,CW,1)
      U0=-CM*GAMI ! perfect lattice
      E0=1.5d0*U1*TPT ! zero-point energy
      Uth=U+E0
      Fth=F+E0
      Uharm=U0+Uth
      Fharm=U0+Fth
      Pharm=U0/3.d0+Uth/2.d0
      PDTharm=.5d0*CVth
      PDRharm=U0/2.25d0+.75d0*Uth-.25d0*CVth
      return
      end

      subroutine HLfit12(eta,F,U,CV,S,U1,CW,LATTICE)
*                                                       Version 24.04.12
* Stems from HLfit8 v.03.12.08;
*   differences: E0 excluded from  U and F;
*   U1 and d(CV)/d\ln(T) are added on the output.
* Fit to thermal part of the thermodynamic functions.
* Baiko, Potekhin, & Yakovlev (2001).
* Zero-point lattice quantum energy 1.5u_1\eta EXCLUDED (unlike HLfit8).
* Input: eta=Tp/T, LATTICE=1 for bcc, 2 for fcc
* Output: F and U (normalized to NkT) - due to phonon excitations,
*   CV and S (normalized to Nk) in the HL model,
*   U1 - the 1st phonon moment,
*   CW=d(CV)/d\ln(T)
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter(EPS=1.d-5,TINY=1.d-99)
      if (LATTICE.eq.1) then ! bcc lattice
         CLM=-2.49389d0 ! 3*ln<\omega/\omega_p>
         U1=.5113875d0
         ALPHA=.265764d0
         BETA=.334547d0
         GAMMA=.932446d0
         A1=.1839d0
         A2=.593586d0
         A3=.0054814d0
         A4=5.01813d-4
         A6=3.9247d-7
         A8=5.8356d-11
         B0=261.66d0
         B2=7.07997d0
         B4=.0409484d0
         B5=.000397355d0
         B6=5.11148d-5
         B7=2.19749d-6
         C9=.004757014d0
         C11=.0047770935d0
      elseif (LATTICE.eq.2) then ! fcc lattice
         CLM=-2.45373d0
         U1=.513194d0
         ALPHA=.257591d0
         BETA=.365284d0
         GAMMA=.9167070d0
         A1=.0
         A2=.532535d0
         A3=.0
         A4=3.76545d-4
         A6=2.63013d-7
         A8=6.6318d-11
         B0=303.20d0
         B2=7.7255d0
         B4=.0439597d0
         B5=.000114295d0
         B6=5.63434d-5
         B7=1.36488d-6
         C9=.00492387d0
         C11=.00437506d0
      else
         stop'HLfit: unknown lattice type'
      endif
      if (eta.gt.1./EPS) then ! asymptote of Eq.(13) of BPY'01
         U=3./(C11*eta**3)
         F=-U/3.
         CV=4.*U
        goto 50
      elseif (eta.lt.EPS) then ! Eq.(17) of BPY'01
        if (eta.lt.TINY) stop'HLfit: eta is too small'
         F=3.*dlog(eta)+CLM-1.5*U1*eta+eta**2/24. 
         U=3.-1.5*U1*eta+eta**2/12.
         CV=3.-eta**2/12.
         goto 50
      endif
      eta2=eta**2
      eta3=eta2*eta
      eta4=eta3*eta
      eta5=eta4*eta
      eta6=eta5*eta
      eta7=eta6*eta
      eta8=eta7*eta
      B9=A6*C9
      B11=A8*C11
      UP=1.+A1*eta+A2*eta2+A3*eta3+A4*eta4+A6*eta6+A8*eta8
      DN=B0+B2*eta2+B4*eta4+B5*eta5+B6*eta6+
     +  B7*eta7+eta8*(B9*eta+B11*eta3)
      EA=dexp(-ALPHA*eta)
      EB=dexp(-BETA*eta)
      EG=dexp(-GAMMA*eta)
      F=dlog(1.d0-EA)+dlog(1.d0-EB)+dlog(1.-EG)-UP/DN ! F_{thermal}/NT
      UP1=A1+
     + 2.*A2*eta+3.*A3*eta2+4.*A4*eta3+6.*A6*eta5+8.*A8*eta7
      UP2=2.*A2+6.*A3*eta+12.*A4*eta2+30.*A6*eta4+56.*A8*eta6
      UP3=6.*A3+24.*A4*eta+120.*A6*eta3+336*A8*eta5
      DN1=2.*B2*eta+4.*B4*eta3+5.*B5*eta4+6.*B6*eta5+
     +  7.*B7*eta6+eta8*(9.*B9+11.*B11*eta2)
      DN2=2.*B2+12.*B4*eta2+20.*B5*eta3+30.*B6*eta4+
     +  42.*B7*eta5+72.*B9*eta7+110.*B11*eta8*eta
      DN3=24.*B4*eta+60.*B5*eta2+120.*B6*eta3+
     +  210.*B7*eta4+504.*B9*eta6+990.*B11*eta8
      DF1=ALPHA*EA/(1.d0-EA)+BETA*EB/(1.d0-EB)+GAMMA*EG/(1.d0-EG)-
     -  (UP1*DN-DN1*UP)/DN**2 ! int.en./NT/eta = df/d\eta
      DF2=ALPHA**2*EA/(1.d0-EA)**2+BETA**2*EB/(1.d0-EB)**2+
     +  GAMMA**2*EG/(1.d0-EG)**2+
     +  ((UP2*DN-DN2*UP)*DN-2.*(UP1*DN-DN1*UP)*DN1)/DN**3 ! -d2f/d\eta^2
      U=DF1*eta
      CV=DF2*eta2
      DF3=-ALPHA**3*EA/(1.d0-EA)**3*(1.+EA)-
     -  BETA**3*EB/(1.d0-EB)**3*(1.+EB)-
     -  GAMMA**3*EG/(1.d0-EG)**3*(1.+EG)+
     +  UP3/DN-(3.*UP2*DN1+3.*UP1*DN2+UP*DN3)/DN**2+
     +  6.*DN1*(UP1*DN1+UP*DN2)/DN**3-6.*UP*DN1**3/DN**4 ! -d3f/d\eta^3
      CW=-2.*CV-eta3*DF3
   50 continue
      S=U-F
      return
      end
	
      subroutine CORMIX(RS,GAME,Zmean,Z2mean,Z52,Z53,Z321,
     *  FMIX,UMIX,PMIX,CVMIX,PDTMIX,PDRMIX)
*                                                       Version 02.07.09
* Correction to the linear mixing rule for moderate to small Gamma
* Input: RS=r_s (if RS=0, then OCP, otherwise EIP)
*        GAME=\Gamma_e
*        Zmean=<Z> (average Z of all ions, without electrons)
*        Z2mean=<Z^2>, Z52=<Z^2.5>, Z53=<Z^{5/3}>, Z321=<Z(Z+1)^1.5>
* Output: FMIX=\Delta f - corr.to the reduced free energy f=F/N_{ion}kT
*         UMIX=\Delta u - corr.to the reduced internal energy u
*         PMIX=\Delta u - corr.to the reduced pressure P=P/n_{ion}kT
*         CVMIX=\Delta c - corr.to the reduced heat capacity c_V
*         PDTMIX=(1/n_{ion}kT)d\Delta P / d ln T
*               = \Delta p +  d \Delta p / d ln T
*         PDRMIX=(1/n_{ion}kT)d\Delta P / d ln n_e
* (composition is assumed fixed: Zmean,Z2mean,Z52,Z53=constant)
      implicit double precision (A-H), double precision (O-Z)
      parameter (TINY=1.d-9)
      GAMImean=GAME*Z53
      if (RS.lt.TINY) then ! OCP
         Dif0=Z52-dsqrt(Z2mean**3/Zmean)
      else
         Dif0=Z321-dsqrt((Z2mean+Zmean)**3/Zmean)
      endif
      DifR=Dif0/Z52
      DifFDH=Dif0*GAME*sqrt(GAME/3.) ! F_DH - F_LM(DH)
      D=Z2mean/Zmean**2
      if (dabs(D-1.d0).lt.TINY) then ! no correction
         FMIX=0.
         UMIX=0.
         PMIX=0.
         CVMIX=0.
         PDTMIX=0.
         PDRMIX=0.
         return
      endif
      P3=D**(-0.2)
      D0=(2.6*DifR+14.*DifR**3)/(1.d0-P3)
      GP=D0*GAMImean**P3
      FMIX0=DifFDH/(1.+GP)
      Q=D**2*.0117
      R=1.5/P3-1.
      GQ=Q*GP
      FMIX=FMIX0/(1.+GQ)**R
      G=1.5-P3*GP/(1.+GP)-R*P3*GQ/(1.+GQ)
      UMIX=FMIX*G
      PMIX=UMIX/3.d0
      GDG=-P3**2*(GP/(1.d0+GP)**2+R*GQ/(1.d0+GQ)**2) ! d G /d ln Gamma
      UDG=UMIX*G+FMIX*GDG ! d u_mix /d ln Gamma
      CVMIX=UMIX-UDG
      PDTMIX=PMIX-UDG/3.
      PDRMIX=PMIX+UDG/9.
      return
      end

* ===================  IDEAL ELECTRON GAS  =========================== *
      subroutine ELECT11(TEMP,CHI,
     *  DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,
     *  DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
*                                                       Version 17.11.11
* ELECT9 v.04.03.09 + smooth match of two fits at chi >> 1 + add.outputs
* Compared to ELECTRON v.06.07.00, this S/R is completely rewritten: 
*        numerical differentiation is avoided now.
* Compared to ELECT7 v.06.06.07,
*    - call BLIN7 is changed to call BLIN9,
*    - Sommerfeld expansion is used at chi >~ 28 i.o. 1.e4
*    - Sommerfeld expansion is corrected: introduced DeltaEF, D1 and D2.
* Ideal electron-gas EOS.
* Input: TEMP - T [a.u.], CHI=\mu/T
* Output: DENS - electron number density n_e [a.u.],
*         FEid - free energy / N_e kT, UEid - internal energy / N_e kT,
*         PEid - pressure (P_e) / n_e kT, SEid - entropy / N_e k,
*         CVE - heat capacity / N_e k,
*         CHITE=(d ln P_e/d ln T)_V, CHIRE=(d ln P_e/d ln n_e)_T
*         DlnDH=(d ln n_e/d CHI)_T = (T/n_e) (d n_e/d\mu)_T
*         DlnDT=(d ln n_e/d ln T)_CHI
*         DlnDHH=(d^2 ln n_e/d CHI^2)_T
*         DlnDTT=(d^2 ln n_e/d (ln T)^2)_CHI
*         DlnDHT=d^2 ln n_e/d (ln T) d CHI
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter (CHI2=28.d0,XMAX=20.d0)
      parameter (DCHI2=CHI2-1.d0)
      parameter (XSCAL2=XMAX/DCHI2)
      X2=(CHI-CHI2)*XSCAL2
      if (X2.lt.-XMAX) then
         call ELECT11a(TEMP,CHI,
     *     DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,
     *     DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
      elseif (X2.gt.XMAX) then
         call ELECT11b(TEMP,CHI,
     *     DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,
     *     DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
      else
         call FERMI10(X2,XMAX,FP,FM)
         call ELECT11a(TEMP,CHI,
     *     DENSa,FEida,PEida,UEida,SEida,CVEa,CHITEa,CHIREa,
     *     DlnDHa,DlnDTa,DlnDHHa,DlnDTTa,DlnDHTa)
         call ELECT11b(TEMP,CHI,
     *     DENSb,FEidb,PEidb,UEidb,SEidb,CVEb,CHITEb,CHIREb,
     *     DlnDHb,DlnDTb,DlnDHHb,DlnDTTb,DlnDHTb)
         DENS=DENSa*FP+DENSb*FM
         FEid=FEida*FP+FEidb*FM
         PEid=PEida*FP+PEidb*FM
         UEid=UEida*FP+UEidb*FM
         SEid=SEida*FP+SEidb*FM
         CVE=CVEa*FP+CVEb*FM
         CHITE=CHITEa*FP+CHITEb*FM
         CHIRE=CHIREa*FP+CHIREb*FM
         DlnDH=DlnDHa*FP+DlnDHb*FM
         DlnDT=DlnDTa*FP+DlnDTb*FM
         DlnDHH=DlnDHHa*FP+DlnDHHb*FM
         DlnDHT=DlnDHTa*FP+DlnDHTb*FM
         DlnDTT=DlnDTTa*FP+DlnDTTb*FM
      endif
      return
      end

      subroutine ELECT11a(TEMP,CHI,
     *  DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,
     *  DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
*                                                       Version 16.11.11
* This is THE FIRST PART of ELECT9 v.04.03.09.
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter (BOHR=137.036,PI=3.141592653d0)
      parameter (PI2=PI**2,BOHR2=BOHR**2,BOHR3=BOHR2*BOHR) !cleaned 15/6
      TEMR=TEMP/BOHR2 ! T in rel.units (=T/mc^2)
      call BLIN9(TEMR,CHI,
     *  W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT,
     *  W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,
     *  W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,
     *  W0XXX,W0XTT,W0XXT)
      TPI=TEMR*dsqrt(2.d0*TEMR)/PI2 ! common pre-factor
      DENR=TPI*(W1*TEMR+W0)
      PR=TEMR*TPI/3.*(W2*TEMR+2.*W1)
      U=TEMR*TPI*(W2*TEMR+W1)
* (these are density, pressure, and internal energy in the rel.units)
      PEid=PR/(DENR*TEMR)
      UEid=U/(DENR*TEMR)
      FEid=CHI-PEid
      DENS=DENR*BOHR3 ! converts from rel.units to a.u.
      SEid=UEid-FEid
* derivatives over T at constant chi:
      dndT=TPI*(1.5*W0/TEMR+2.5*W1+W0DT+TEMR*W1DT) ! (d n_e/dT)_\chi
      dPdT=TPI/3.*(5.*W1+2.*TEMR*W1DT+3.5*TEMR*W2+TEMR**2*W2DT)!dP/dT
      dUdT=TPI*(2.5*W1+TEMR*W1DT+3.5*TEMR*W2+TEMR**2*W2DT)!dU/dT_\chi
* derivatives over chi at constant T and second derivatives:
      dndH=TPI*(W0DX+TEMR*W1DX) ! (d n_e/d\chi)_T
      dndHH=TPI*(W0DXX+TEMR*W1DXX) ! (d^2 n_e/d\chi)_T
      dndTT=TPI*(.75*W0/TEMR**2+3.*W0DT/TEMR+W0DTT+
     +  3.75*W1/TEMR+5.*W1DT+TEMR*W1DTT)
      dndHT=TPI*(1.5*W0DX/TEMR+W0DXT+2.5*W1DX+TEMR*W1DXT)
      DlnDH=dndH/DENR ! (d ln n_e/d\chi)_T
      DlnDT=dndT*TEMR/DENR ! (d ln n_e/d ln T)_\chi
      DlnDHH=dndHH/DENR-DlnDH**2 ! (d^2 ln n_e/d\chi^2)_T
      DlnDTT=TEMR**2/DENR*dndTT+DlnDT-DlnDT**2 ! d^2 ln n_e/d ln T^2
      DlnDHT=TEMR/DENR*(dndHT-dndT*DlnDH) ! d^2 ln n_e/d\chi d ln T
      dPdH=TPI/3.*TEMR*(2.*W1DX+TEMR*W2DX) ! (d P_e/d\chi)_T
      dUdH=TPI*TEMR*(W1DX+TEMR*W2DX) ! (d U_e/d\chi)_T
      CVE=(dUdT-dUdH*dndT/dndH)/DENR
      CHITE=TEMR/PR*(dPdT-dPdH*dndT/dndH)
      CHIRE=DENR/PR*dPdH/dndH ! (dndH*TEMR*PEid) ! DENS/PRE*dPdH/dndH
      return
      end

      subroutine ELECT11b(TEMP,CHI,
     *  DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,
     *  DlnDH,DlnDT,DlnDHH,DlnDTT,DlnDHT)
*                                                       Version 17.11.11
* Stems from ELECT9b v.19.01.10, Diff. - additional output.
* Sommerfeld expansion at very large CHI.
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter (BOHR=137.036,PI=3.141592653d0)
      parameter (PI2=PI**2,BOHR2=BOHR**2,BOHR3=BOHR2*BOHR) !cleaned 15/6
      TEMR=TEMP/BOHR2 ! T in rel.units (=T/mc^2)
      EF=CHI*TEMR ! Fermi energy in mc^2 - zeroth aprox. = CMU1
      DeltaEF=PI2*TEMR**2/6.d0*(1.d0+2.d0*EF*(2.d0+EF))/
     /  (EF*(1.d0+EF)*(2.d0+EF)) ! corr. [p.125, equiv.Eq.(6) of PC'10]
      EF=EF+DeltaEF ! corrected Fermi energy (14.02.09)
      G=1.d0+EF ! electron Lorentz-factor
      if (EF.gt.1.d-5) then ! relativistic expansion (Yak.&Shal.'89)
        PF=dsqrt(G**2-1.d0) ! Fermi momentum [rel.un.=mc]
        F=(PF*(1.+2.d0*PF**2)*G-PF**3/.375d0-dlog(PF+G))/8.d0/PI2!F/V
        DF=-TEMR**2*PF*G/6.d0 ! thermal correction to F/V
        P=(PF*G*(PF**2/1.5d0-1.d0)+dlog(PF+G))/8.d0/PI2 ! P(T=0)
        DP=TEMR**2*PF*(PF**2+2.d0)/G/18.d0 ! thermal correction to P
        CVE=PI2*TEMR*G/PF**2
      else ! nonrelativistic limit
        PF=dsqrt(2.d0*EF)
        F=PF**5*0.1d0/PI2
        DF=-TEMR**2*PF/6.d0
        P=F/1.5d0
        DP=TEMR**2*PF/9.d0
        CVE=PI2*TEMR/EF/2.d0
      endif
      F=F+DF
      P=P+DP
      S=-2.d0*DF ! entropy per unit volume [rel.un.]
      U=F+S
      CHIRE=PF**5/(9.d0*PI2*P*G)
      CHITE=2.d0*DP/P
      DENR=PF**3/3.d0/PI2 ! n_e [rel.un.=\Compton^{-3}]
      DENS=DENR*BOHR3 ! conversion to a.u.(=\Bohr_radius^{-3})
* derivatives over chi at constant T and T at constant chi:
      TPI=TEMR*dsqrt(2.d0*TEMR)/PI2 ! common pre-factor
      call SOMMERF(TEMR,CHI,
     *  W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT,
     *  W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,
     *  W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,
     *  W0XXX,W0XTT,W0XXT)
      dndH=TPI*(W0DX+TEMR*W1DX) ! (d n_e/d\chi)_T
      dndT=TPI*(1.5*W0/TEMR+2.5*W1+W0DT+TEMR*W1DT) ! (d n_e/dT)_\chi
      dndHH=TPI*(W0DXX+TEMR*W1DXX) ! (d^2 n_e/d\chi)_T
      dndTT=TPI*(.75*W0/TEMR**2+3.*W0DT/TEMR+W0DTT+
     +  3.75*W1/TEMR+5.*W1DT+TEMR*W1DTT)
      dndHT=TPI*(1.5*W0DX/TEMR+W0DXT+2.5*W1DX+TEMR*W1DXT)
      DlnDH=dndH/DENR ! (d ln n_e/d\chi)_T
      DlnDT=dndT*TEMR/DENR ! (d ln n_e/d ln T)_\chi
      DlnDHH=dndHH/DENR-DlnDH**2 ! (d^2 ln n_e/d\chi^2)_T
      DlnDTT=TEMR**2/DENR*dndTT+DlnDT-DlnDT**2 ! d^2 ln n_e/d ln T^2
      DlnDHT=TEMR/DENR*(dndHT-dndT*DlnDH) ! d^2 ln n_e/d\chi d ln T
      DT=DENR*TEMR
      PEid=P/DT
      UEid=U/DT
      FEid=F/DT
      SEid=S/DT
* Empirical corrections of 16.02.09:
      D1=DeltaEF/EF
      D2=D1*(4.d0-2.d0*(PF/G))
      CVE=CVE/(1.d0+D2)
      SEid=SEid/(1.d0+D1)
      CHITE=CHITE/(1.d0+D2)
      return
      end

      subroutine SOMMERF(TEMR,CHI,
     *  W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT,
     *  W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,
     *  W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,
     *  W0XXX,W0XTT,W0XXT)
*                                                       Version 17.11.11
* Sommerfeld expansion for the Fermi-Dirac integrals
* Input: TEMR=T/mc^2; CHI=(\mu-mc^2)/T
* Output: Wk - Fermi-Dirac integral of the order k+1/2
*         WkDX=dWk/dCHI, WkDT = dWk/dT, WkDXX=d^2 Wk / d CHI^2,
*         WkDTT=d^2 Wk / d T^2, WkDXT=d^2 Wk /dCHIdT,
*         W0XXX=d^3 W0 / d CHI^3, W0XTT=d^3 W0 /(d CHI d^2 T),
*         W0XXT=d^3 W0 /dCHI^2 dT
* [Draft source: yellow book pages 124-127]
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter(PI=3.141592653d0,EPS=1.d-4)
      parameter(PI2=PI**2)
      if (CHI.lt..5d0) stop'SOMMERF: non-degenerate (small CHI)'
      if (TEMR.le.0.d0) stop'SOMMERF: T < 0'
      CMU1=CHI*TEMR ! chemical potential in rel.units
      CMU=1.d0+CMU1
      call SUBFERMJ(CMU1,
     *  CJ00,CJ10,CJ20,
     *  CJ01,CJ11,CJ21,
     *  CJ02,CJ12,CJ22,
     *  CJ03,CJ13,CJ23,
     *  CJ04,CJ14,CJ24,CJ05)
      PIT26=(PI*TEMR)**2/6.d0
CCC      PITAU4=PIT26**2*0.7d0
      CN0=dsqrt(.5d0/TEMR)/TEMR
      CN1=CN0/TEMR
      CN2=CN1/TEMR
      W0=CN0*(CJ00+PIT26*CJ02) ! +CN0*PITAU4*CJ04
      W1=CN1*(CJ10+PIT26*CJ12) ! +CN1*PITAU4*CJ14
      W2=CN2*(CJ20+PIT26*CJ22) ! +CN2*PITAU4*CJ24
      W0DX=CN0*TEMR*(CJ01+PIT26*CJ03) ! +CN0*PITAU4*CJ05
      W1DX=CN0*(CJ11+PIT26*CJ13)
      W2DX=CN1*(CJ21+PIT26*CJ23)
      W0DT=CN1*(CMU1*CJ01-1.5d0*CJ00+PIT26*(CMU1*CJ03+.5d0*CJ02))
CCC     +  CN1*PITAU4*(CMU1*CJ05+2.5d0*CJ04)
      W1DT=CN2*(CMU1*CJ11-2.5d0*CJ10+PIT26*(CMU1*CJ13-.5d0*CJ12))
      W2DT=CN2/TEMR*(CMU1*CJ21-3.5d0*CJ20+PIT26*(CMU1*CJ23-1.5d0*CJ22))
      W0DXX=CN0*TEMR**2*(CJ02+PIT26*CJ04)
      W1DXX=CN0*TEMR*(CJ12+PIT26*CJ14)
      W2DXX=CN0*(CJ22+PIT26*CJ24)
      W0DXT=CN0*(CMU1*CJ02-.5d0*CJ01+PIT26*(CMU1*CJ04+1.5d0*CJ03))
      W1DXT=CN1*(CMU1*CJ12-1.5d0*CJ11+PIT26*(CMU1*CJ14+.5d0*CJ13))
      W2DXT=CN2*(CMU1*CJ22-2.5d0*CJ21+PIT26*(CMU1*CJ24-.5d0*CJ23))
      W0DTT=CN2*(3.75d0*CJ00-3.d0*CMU1*CJ01+CMU1**2*CJ02+
     +  PIT26*(-.25d0*CJ02+CMU1*CJ03+CMU1**2*CJ04))
      W1DTT=CN2/TEMR*(8.75d0*CJ10-5.d0*CMU1*CJ11+CMU1**2*CJ12+
     +  PIT26*(.75d0*CJ12-CMU1*CJ13+CMU1**2*CJ14))
      W2DTT=CN2/TEMR**2*(15.75d0*CJ20-7.d0*CMU1*CJ21+CMU1**2*CJ22+
     +  PIT26*(3.75d0*CJ22-3.d0*CMU1*CJ23+CMU1**2*CJ24))
      W0XXX=CN0*TEMR**3*(CJ03+PIT26*CJ05)
      W0XXT=CN0*TEMR*(CMU1*CJ03+.5d0*CJ02+PIT26*(CMU1*CJ05+2.5d0*CJ04))
      W0XTT=CN1*(.75d0*CJ01-CMU1*CJ02+CMU1**2*CJ03+
     +  PIT26*(.75d0*CJ03+3.d0*CMU1*CJ04+CMU1**2*CJ05))
      return
      end

      subroutine SUBFERMJ(CMU1,
     *  CJ00,CJ10,CJ20,
     *  CJ01,CJ11,CJ21,
     *  CJ02,CJ12,CJ22,
     *  CJ03,CJ13,CJ23,
     *  CJ04,CJ14,CJ24,CJ05)
*                                                       Version 17.11.11
* Supplement to SOMMERF
      implicit double precision (A-H), double precision (O-Z)
      save
      if (CMU1.le.0.d0) stop'SUBFERMJ: small CHI'
      CMU=1.d0+CMU1
      X0=dsqrt(CMU1*(2.d0+CMU1))
      X3=X0**3
      X5=X0**5
      if (X0.lt.EPS) then
         CJ00=X3/3.d0
         CJ10=.1d0*X5
         CJ20=X0**7/28.d0
      else
         CL=dlog(X0+CMU)
         CJ00=.5d0*(X0*CMU-CL) ! J_{1/2}^0
         CJ10=X3/3.d0-CJ00 ! J_{3/2}^0
         CJ20=(.75d0*CMU-2.d0)/3.d0*X3+1.25d0*CJ00 ! J_{5/2}^0
      endif
      CJ01=X0 ! J_{1/2}^1
      CJ11=CJ01*CMU1 ! J_{3/2}^1
      CJ21=CJ11*CMU1 ! J_{5/2}^1
      CJ02=CMU/X0 ! J_{1/2}^2
      CJ12=CMU1/X0*(3.d0+2.d0*CMU1) ! J_{3/2}^2
      CJ22=CMU1**2/X0*(5.d0+3.d0*CMU1) ! J_{5/2}^2
      CJ03=-1.d0/X3 ! J_{1/2}^3
      CJ13=CMU1/X3*(2.d0*CMU1**2+6.d0*CMU1+3.d0)
      CJ23=CMU1**2/X3*(6.d0*CMU1**2+2.d1*CMU1+1.5d1)
      CJ04=3.d0*CMU/X5
      CJ14=-3.d0*CMU1/X5
      CJ24=CMU1**2/X5*(6.d0*CMU1**3+3.d1*CMU1**2+45.d0*CMU1+15.d0)
      CJ05=(-12.d0*CMU1**2-24.d0*CMU1-15.d0)/(X5*X0**2)
      return
      end

      subroutine FERMI10(X,XMAX,FP,FM)
*                                                       Version 20.01.10
* Fermi distribution function and its 3 derivatives
* Input: X - argument f(x)
*        XMAX - max|X| where it is assumed that 0 < f(x) < 1.
* Output: FP = f(x)
*         FM = 1-f(x)
      implicit double precision (A-H), double precision (O-Z)
      save
      if (XMAX.lt.3.d0) stop'FERMI10: XMAX'
      if (X.gt.XMAX) then
         FP=0.d0
         FM=1.d0
      elseif (X.lt.-XMAX) then
         FP=1.d0
         FM=0.d0
      else
         FP=1.d0/(dexp(X)+1.d0)
         FM=1.d0-FP
      endif
      return
      end

* ==============  ELECTRON EXCHANGE AND CORRELATION   ================ *
      subroutine EXCOR7(RS,GAME,FXC,UXC,PXC,CVXC,SXC,PDTXC,PDRXC)
*                                                       Version 09.06.07
* Exchange-correlation contribution for the electron gas
* Stems from TANAKA1 v.03.03.96. Added derivatives.
* Input: RS - electron density parameter =electron-sphere radius in a.u.
*        GAME - electron Coulomb coupling parameter
* Output: FXC - excess free energy of e-liquid per kT per one electron
*             according to Tanaka & Ichimaru 85-87 and Ichimaru 93 FXC,
*         UXC - free and internal energy contr.[per 1 electron, kT]
*         PXC - pressure contribution divided by (n_e kT)
*         CVXC - heat capacity divided by N_e k
*         SXC - entropy divided by N_e k
*         PDTXC,PDRXC = PXC + d ln PXC / d ln(T,\rho)
      implicit double precision(A-H),double precision(O-Z)
      save
      THETA=.543*RS/GAME ! non-relativistic degeneracy parameter
      SQTH=dsqrt(THETA)
      THETA2=THETA**2
      THETA3=THETA2*THETA
      THETA4=THETA3*THETA
      if (THETA.gt..005) then
         CHT1=dcosh(1.d0/THETA)
         SHT1=dsinh(1.d0/THETA)
         CHT2=dcosh(1.d0/SQTH)
         SHT2=dsinh(1.d0/SQTH)
         T1=SHT1/CHT1 ! dtanh(1.d0/THETA)
         T2=SHT2/CHT2 ! dtanh(1./dsqrt(THETA))
         T1DH=-1./(THETA*CHT1)**2 ! d T1 / d\theta
         T1DHH=2./(THETA*CHT1)**3*(CHT1-SHT1/THETA)
         T2DH=-.5*SQTH/(THETA*CHT2)**2
         T2DHH=(.75*SQTH*CHT2-.5*SHT2)/(THETA*CHT2)**3
      else
         T1=1.
         T2=1.
         T1DH=0.
         T2DH=0.
         T1DHH=0.
         T2DHH=0.
      endif
      A0=.75+3.04363*THETA2-.09227*THETA3+1.7035*THETA4
      A0DH=6.08726*THETA-.27681*THETA2+6.814*THETA3
      A0DHH=6.08726-.55362*THETA+20.442*THETA2
      A1=1.+8.31051*THETA2+5.1105*THETA4
      A1DH=16.62102*THETA+20.442*THETA3
      A1DHH=16.62102+61.326*THETA2
      A=.610887*A0/A1*T1 ! HF fit of Perrot and Dharma-wardana
      AH=A0DH/A0-A1DH/A1+T1DH/T1
      ADH=A*AH
      ADHH=ADH*AH+A*(A0DHH/A0-(A0DH/A0)**2-A1DHH/A1+(A1DH/A1)**2+
     +  T1DHH/T1-(T1DH/T1)**2)
      B0=.341308+12.070873d0*THETA2+1.148889d0*THETA4
      B0DH=24.141746d0*THETA+4.595556d0*THETA3
      B0DHH=24.141746d0+13.786668d0*THETA2
      B1=1.+10.495346d0*THETA2+1.326623*THETA4
      B1DH=20.990692d0*THETA+5.306492*THETA3
      B1DHH=20.990692d0+15.919476d0*THETA2
      B=SQTH*T2*B0/B1
      BH=.5/THETA+T2DH/T2+B0DH/B0-B1DH/B1
      BDH=B*BH
      BDHH=BDH*BH+B*(-.5/THETA2+T2DHH/T2-(T2DH/T2)**2+
     +  B0DHH/B0-(B0DH/B0)**2-B1DHH/B1+(B1DH/B1)**2)
      D0=.614925+16.996055d0*THETA2+1.489056*THETA4
      D0DH=33.99211d0*THETA+5.956224d0*THETA3
      D0DHH=33.99211d0+17.868672d0*THETA2
      D1=1.+10.10935*THETA2+1.22184*THETA4
      D1DH=20.2187*THETA+4.88736*THETA3
      D1DHH=20.2187+14.66208*THETA2
      D=SQTH*T2*D0/D1
      DH=.5/THETA+T2DH/T2+D0DH/D0-D1DH/D1
      DDH=D*DH
      DDHH=DDH*DH+D*(-.5/THETA2+T2DHH/T2-(T2DH/T2)**2+
     +  D0DHH/D0-(D0DH/D0)**2-D1DHH/D1+(D1DH/D1)**2)
      E0=.539409+2.522206*THETA2+.178484*THETA4
      E0DH=5.044412*THETA+.713936*THETA3
      E0DHH=5.044412+2.141808*THETA2
      E1=1.+2.555501*THETA2+.146319*THETA4
      E1DH=5.111002*THETA+.585276*THETA3
      E1DHH=5.111002+1.755828*THETA2
      E=THETA*T1*E0/E1
      EH=1./THETA+T1DH/T1+E0DH/E0-E1DH/E1
      EDH=E*EH
      EDHH=EDH*EH+E*(T1DHH/T1-(T1DH/T1)**2+E0DHH/E0-(E0DH/E0)**2-
     -  E1DHH/E1+(E1DH/E1)**2-1./THETA2)
      EXP1TH=dexp(-1./THETA)
      C=(.872496+.025248*EXP1TH)*E
      CDH=.025248*EXP1TH/THETA2*E+C*EDH/E
      CDHH=.025248*EXP1TH/THETA2*(EDH+(1.-2.*THETA)/THETA2*E)+
     +  CDH*EDH/E+C*EDHH/E-C*(EDH/E)**2
      DISCR=dsqrt(4.*E-D**2)
      DIDH=.5/DISCR*(4.*EDH-2.*D*DDH)
      DIDHH=(-((2.*EDH-D*DDH)/DISCR)**2+2.*EDHH-DDH**2-D*DDHH)/DISCR
      S1=-C/E*GAME
      S1H=CDH/C-EDH/E
      S1DH=S1*S1H
      S1DHH=S1DH*S1H+S1*(CDHH/C-(CDH/C)**2-EDHH/E+(EDH/E)**2)
      S1DG=-C/E ! => S1DGG=0
      S1DHG=S1DG*(CDH/C-EDH/E)
      B2=B-C*D/E
      B2DH=BDH-(CDH*D+C*DDH)/E+C*D*EDH/E**2
      B2DHH=BDHH-(CDHH*D+2.*CDH*DDH+C*DDHH)/E+
     +  (2.*(CDH*D+C*DDH-C*D*EDH/E)*EDH+C*D*EDHH)/E**2
      SQGE=dsqrt(GAME)
      S2=-2./E*B2*SQGE
      S2H=B2DH/B2-EDH/E
      S2DH=S2*S2H
      S2DHH=S2DH*S2H+S2*(B2DHH/B2-(B2DH/B2)**2-EDHH/E+(EDH/E)**2)
      S2DG=.5*S2/GAME
      S2DGG=-.5*S2DG/GAME
      S2DHG=.5*S2DH/GAME
      R3=E*GAME+D*SQGE+1.
      R3DH=EDH*GAME+DDH*SQGE
      R3DHH=EDHH*GAME+DDHH*SQGE
      R3DG=E+.5*D/SQGE
      R3DGG=-.25*D/(GAME*SQGE)
      R3DHG=EDH+.5*DDH/SQGE
      B3=A-C/E
      B3DH=ADH-CDH/E+C*EDH/E**2
      B3DHH=ADHH-CDHH/E+(2.*CDH*EDH+C*EDHH)/E**2-2.*C*EDH**2/E**3
      C3=(D/E*B2-B3)/E ! =D*B2/E**2-B3/E
      C3DH=(DDH*B2+D*B2DH+B3*EDH)/E**2-2.*D*B2*EDH/E**3-B3DH/E
      C3DHH=(-B3DHH+
     +  (DDHH*B2+2.*DDH*B2DH+D*B2DHH+B3DH*EDH+B3*EDHH+B3DH*EDH)/E-
     -  2.*((DDH*B2+D*B2DH+B3*EDH+DDH*B2+D*B2DH)*EDH+D*B2*EDHH)/E**2+
     +  6.*D*B2*EDH**2/E**3)/E
      S3=C3*dlog(R3)
      S3DH=S3*C3DH/C3+C3*R3DH/R3
      S3DHH=(S3DH*C3DH+S3*C3DHH)/C3-S3*(C3DH/C3)**2+
     +  (C3DH*R3DH+C3*R3DHH)/R3-C3*(R3DH/R3)**2
      S3DG=C3*R3DG/R3
      S3DGG=C3*(R3DGG/R3-(R3DG/R3)**2)
      S3DHG=(C3DH*R3DG+C3*R3DHG)/R3-C3*R3DG*R3DH/R3**2
      B4=2.-D**2/E
      B4DH=EDH*(D/E)**2-2.*D*DDH/E
      B4DHH=EDHH*(D/E)**2+2.*EDH*(D/E)**2*(DDH/D-EDH/E)-
     -  2.*(DDH**2+D*DDHH)/E+2.*D*DDH*EDH/E**2
      C4=2.*E*SQGE+D
      C4DH=2.*EDH*SQGE+DDH
      C4DHH=2.*EDHH*SQGE+DDHH
      C4DG=E/SQGE
      C4DGG=-.5*E/(GAME*SQGE)
      C4DHG=EDH/SQGE
      S4A=2./E/DISCR
      S4AH=EDH/E+DIDH/DISCR
      S4ADH=-S4A*S4AH
      S4ADHH=-S4ADH*S4AH-
     -  S4A*(EDHH/E-(EDH/E)**2+DIDHH/DISCR-(DIDH/DISCR)**2)
      S4B=D*B3+B4*B2
      S4BDH=DDH*B3+D*B3DH+B4DH*B2+B4*B2DH
      S4BDHH=DDHH*B3+2.*DDH*B3DH+D*B3DHH+B4DHH*B2+2.*B4DH*B2DH+B4*B2DHH
      S4C=datan(C4/DISCR)-datan(D/DISCR)
      UP1=C4DH*DISCR-C4*DIDH
      DN1=DISCR**2+C4**2
      UP2=DDH*DISCR-D*DIDH
      DN2=DISCR**2+D**2
      S4CDH=UP1/DN1-UP2/DN2
      S4CDHH=(C4DHH*DISCR-C4*DIDHH)/DN1-
     -  UP1*2.*(DISCR*DIDH+C4*C4DH)/DN1**2-
     -  (DDHH*DISCR-D*DIDHH)/DN2+UP2*2.*(DISCR*DIDH+D*DDH)/DN2**2
      S4CDG=C4DG*DISCR/DN1
      S4CDGG=C4DGG*DISCR/DN1-2.*C4*DISCR*(C4DG/DN1)**2
      S4CDHG=(C4DHG*DISCR+C4DG*DIDH-
     -  C4DG*DISCR/DN1*2.*(DISCR*DIDH+C4*C4DH))/DN1
      S4=S4A*S4B*S4C
      S4DH=S4ADH*S4B*S4C+S4A*S4BDH*S4C+S4A*S4B*S4CDH
      S4DHH=S4ADHH*S4B*S4C+S4A*S4BDHH*S4C+S4A*S4B*S4CDHH+
     +  2.*(S4ADH*S4BDH*S4C+S4ADH*S4B*S4CDH+S4A*S4BDH*S4CDH)
      S4DG=S4A*S4B*S4CDG
      S4DGG=S4A*S4B*S4CDGG
      S4DHG=S4A*S4B*S4CDHG+S4CDG*(S4ADH*S4B+S4A*S4BDH)
      FXC=S1+S2+S3+S4
      FXCDH=S1DH+S2DH+S3DH+S4DH
      FXCDG=S1DG+S2DG+S3DG+S4DG
      FXCDHH=S1DHH+S2DHH+S3DHH+S4DHH
      FXCDGG=S2DGG+S3DGG+S4DGG
      FXCDHG=S1DHG+S2DHG+S3DHG+S4DHG
      PXC=(GAME*FXCDG-2.*THETA*FXCDH)/3.
      UXC=GAME*FXCDG-THETA*FXCDH
      SXC=(GAME*S2DG-S2+GAME*S3DG-S3+S4A*S4B*(GAME*S4CDG-S4C))-
     -  THETA*FXCDH
      if (dabs(SXC).lt.1.d-9*dabs(THETA*FXCDH)) SXC=0. ! accuracy loss
      CVXC=2.*THETA*(GAME*FXCDHG-FXCDH)-THETA**2*FXCDHH-GAME**2*FXCDGG
      if (dabs(CVXC).lt.1.d-9*dabs(GAME**2*FXCDGG)) CVXC=0. ! accuracy
      PDLH=THETA*(GAME*FXCDHG-2.*FXCDH-2.*THETA*FXCDHH)/3.
      PDLG=GAME*(FXCDG+GAME*FXCDGG-2.*THETA*FXCDHG)/3.
      PDRXC=PXC+(PDLG-2.*PDLH)/3.
      PDTXC=GAME*(THETA*FXCDHG-GAME*FXCDGG/3.)-
     -  THETA*(FXCDH/.75+THETA*FXCDHH/1.5)
      return
      end

* ======================  AUXILIARY SUBROUTINES   ==================== *
      subroutine FERINV7(F,N,X,XDF,XDFF) ! Inverse Fermi intergals
*                                                       Version 24.05.07
* X_q(f)=F^{-1}_q(f) : H.M.Antia 93 ApJS 84, 101
* q=N-1/2=-1/2,1/2,3/2,5/2 (N=0,1,2,3)
* Input: F - argument, N=q+1/2
* Output: X=X_q, XDF=dX/df, XDFF=d^2 X / df^2
* Relative error: N = 0     1      2      3
*        for X:    3.e-9, 4.2e-9, 2.3e-9, 6.2e-9
* jump at f=4:
*         for XDF: 6.e-7, 5.4e-7, 9.6e-8, 3.1e-7
*       for XDFF: 4.7e-5, 4.8e-5, 2.3e-6, 1.5e-6
      implicit double precision (A-H), double precision (O-Z)
      save
      dimension A(0:5,0:3),B(0:6,0:3),C(0:6,0:3),D(0:6,0:3),
     *  LA(0:3),LB(0:3),LD(0:3)
      data A/-1.570044577033d4,1.001958278442d4,-2.805343454951d3,
     *            4.121170498099d2,-3.174780572961d1,1.d0, ! X_{-1/2}
     *    1.999266880833d4,5.702479099336d3,6.610132843877d2,
     *            3.818838129486d1,1.d0,0., ! X_{1/2}
     *    1.715627994191d2,1.125926232897d2,2.056296753055d1,1.d0,0.,0.,
     *    2.138969250409d2,3.539903493971d1,1.d0,0.,0.,0./, ! X_{5/2}
     *  B/-2.782831558471d4,2.886114034012d4,-1.274243093149d4,
     *            3.063252215963d3,-4.225615045074d2,3.168918168284d1,
     *            -1.008561571363d0, ! X_{-1/2}
     *    1.771804140488d4,-2.014785161019d3,9.130355392717d1,
     *            -1.670718177489d0,0.,0.,0., ! X_{1/2}
     *    2.280653583157d2,1.193456203021d2,1.16774311354d1,
     *            -3.226808804038d-1,3.519268762788d-3,0.,0., ! X_{3/2}
     *    7.10854551271d2,9.873746988121d1,1.067755522895d0,
     *            -1.182798726503d-2,0.,0.,0./, ! X_{5/2}
     *  C/2.206779160034d-8,-1.437701234283d-6,6.103116850636d-5,
     *     -1.169411057416d-3,1.814141021608d-2,-9.588603457639d-2,1.d0,
     *  -1.277060388085d-2,7.187946804945d-2,-4.262314235106d-1,
     *      4.997559426872d-1,-1.285579118012d0,-3.930805454272d-1,1.d0,
     *  -6.321828169799d-3,-2.183147266896d-2,-1.05756279932d-1,
     *       -4.657944387545d-1,-5.951932864088d-1,3.6844711771d-1,1.d0,
     *  -3.312041011227d-2,1.315763372315d-1,-4.820942898296d-1,
     *       5.099038074944d-1,5.49561349863d-1,-1.498867562255d0,1.d0/,
     *  D/8.827116613576d-8,-5.750804196059d-6,2.429627688357d-4,
     *       -4.601959491394d-3,6.932122275919d-2,-3.217372489776d-1,
     *       3.124344749296d0, ! X_{-1/2}
     *    -9.745794806288d-3,5.485432756838d-2,-3.29946624326d-1,
     *       4.077841975923d-1,-1.145531476975d0,-6.067091689181d-2,0.,
     *    -4.381942605018d-3,-1.5132365041d-2,-7.850001283886d-2,
     *     -3.407561772612d-1,-5.074812565486d-1,-1.387107009074d-1,0.,
     *    -2.315515517515d-2,9.198776585252d-2,-3.835879295548d-1,
     *       5.415026856351d-1,-3.847241692193d-1,3.739781456585d-2,
     *       -3.008504449098d-2/, ! X_{5/2}
     *  LA/5,4,3,2/,LB/6,3,4,3/,LD/6,5,5,6/
      if (N.lt.0.or.N.gt.3) stop'FERINV7: Invalid subscript'
      if (F.le.0.) stop'FERINV7: Non-positive argument'
      if (F.lt.4.) then
         T=F
         UP=0.
         UP1=0.
         UP2=0.
         DOWN=0.
         DOWN1=0.
         DOWN2=0.
         do I=LA(N),0,-1
            UP=UP*T+A(I,N)
           if (I.ge.1) UP1=UP1*T+A(I,N)*I
           if (I.ge.2) UP2=UP2*T+A(I,N)*I*(I-1)
         enddo
         do I=LB(N),0,-1
            DOWN=DOWN*T+B(I,N)
           if (I.ge.1) DOWN1=DOWN1*T+B(I,N)*I
           if (I.ge.2) DOWN2=DOWN2*T+B(I,N)*I*(I-1)
         enddo
         X=dlog(T*UP/DOWN)
         XDF=1.d0/T+UP1/UP-DOWN1/DOWN
         XDFF=-1.d0/T**2+UP2/UP-(UP1/UP)**2-DOWN2/DOWN+(DOWN1/DOWN)**2
      else
         P=-1./(.5+N) ! = -1/(1+\nu) = power index
         T=F**P ! t - argument of the rational fraction
         T1=P*T/F ! dt/df
         T2=P*(P-1.)*T/F**2 ! d^2 t / df^2
         UP=0.
         UP1=0.
         UP2=0.
         DOWN=0.
         DOWN1=0.
         DOWN2=0.
         do I=6,0,-1
            UP=UP*T+C(I,N)
           if (I.ge.1) UP1=UP1*T+C(I,N)*I
           if (I.ge.2) UP2=UP2*T+C(I,N)*I*(I-1)
         enddo
         do I=LD(N),0,-1
            DOWN=DOWN*T+D(I,N)
           if (I.ge.1) DOWN1=DOWN1*T+D(I,N)*I
           if (I.ge.2) DOWN2=DOWN2*T+D(I,N)*I*(I-1)
         enddo
         R=UP/DOWN
         R1=(UP1-UP*DOWN1/DOWN)/DOWN ! dR/dt
         R2=(UP2-(2.*UP1*DOWN1+UP*DOWN2)/DOWN+2.*UP*(DOWN1/DOWN)**2)/
     /     DOWN
         X=R/T
         RT=(R1-R/T)/T
         XDF=T1*RT
         XDFF=T2*RT+T1**2*(R2-2.*RT)/T
      endif
      return
      end

      subroutine BLIN9(TEMP,CHI,
     *  W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT,
     *  W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,
     *  W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,
     *  W0XXX,W0XTT,W0XXT)
*                                                       Version 21.01.10
* Stems from BLIN8 v.24.12.08
* Difference - smooth matching of different CHI ranges
* Input: TEMP=T/mc^2; CHI=(\mu-mc^2)/T
* Output: Wk - Fermi-Dirac integral of the order k+1/2
*         WkDX=dWk/dCHI, WkDT = dWk/dT, WkDXX=d^2 Wk / d CHI^2,
*         WkDTT=d^2 Wk / d T^2, WkDXT=d^2 Wk /dCHIdT,
*         W0XXX=d^3 W0 / d CHI^3, W0XTT=d^3 W0 /(d CHI d^2 T),
*         W0XXT=d^3 W0 /dCHI^2 dT
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter (CHI1=0.6d0,CHI2=14.d0,XMAX=30.d0)
      parameter (DCHI1=.1d0,DCHI2=CHI2-CHI1-DCHI1)
      parameter (XSCAL1=XMAX/DCHI1,XSCAL2=XMAX/DCHI2)
      X1=(CHI-CHI1)*XSCAL1
      X2=(CHI-CHI2)*XSCAL2
      if (X1.lt.-XMAX) then
         call BLIN9a(TEMP,CHI,
     *     W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT,
     *     W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,
     *     W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,
     *     W0XXX,W0XTT,W0XXT)
      elseif (X2.lt.XMAX) then ! match two fits
        if (X1.lt.XMAX) then ! match fits "a" and "b"
           call FERMI10(X1,XMAX,FP,FM)
           call BLIN9a(TEMP,CHI,
     *       W0a,W0DXa,W0DTa,W0DXXa,W0DTTa,W0DXTa,
     *       W1a,W1DXa,W1DTa,W1DXXa,W1DTTa,W1DXTa,
     *       W2a,W2DXa,W2DTa,W2DXXa,W2DTTa,W2DXTa,
     *       W0XXXa,W0XTTa,W0XXTa)
           call BLIN9b(TEMP,CHI,
     *       W0b,W0DXb,W0DTb,W0DXXb,W0DTTb,W0DXTb,
     *       W1b,W1DXb,W1DTb,W1DXXb,W1DTTb,W1DXTb,
     *       W2b,W2DXb,W2DTb,W2DXXb,W2DTTb,W2DXTb,
     *       W0XXXb,W0XTTb,W0XXTb)
        else ! match fits "b" and "c"
           call FERMI10(X2,XMAX,FP,FM)
           call BLIN9b(TEMP,CHI,
     *       W0a,W0DXa,W0DTa,W0DXXa,W0DTTa,W0DXTa,
     *       W1a,W1DXa,W1DTa,W1DXXa,W1DTTa,W1DXTa,
     *       W2a,W2DXa,W2DTa,W2DXXa,W2DTTa,W2DXTa,
     *       W0XXXa,W0XTTa,W0XXTa)
           call BLIN9c(TEMP,CHI,
     *       W0b,W0DXb,W0DTb,W0DXXb,W0DTTb,W0DXTb,
     *       W1b,W1DXb,W1DTb,W1DXXb,W1DTTb,W1DXTb,
     *       W2b,W2DXb,W2DTb,W2DXXb,W2DTTb,W2DXTb,
     *       W0XXXb,W0XTTb,W0XXTb)
        endif
         W0=W0a*FP+W0b*FM
         W0DX=W0DXa*FP+W0DXb*FM !! +(W0a-W0b)*F1
         W0DT=W0DTa*FP+W0DTb*FM
         W0DXX=W0DXXa*FP+W0DXXb*FM !! +2.d0*(W0DXa-W0DXb)*F1+(W0a-W0b)*F2
         W0DTT=W0DTTa*FP+W0DTTb*FM
         W0DXT=W0DXTa*FP+W0DXTb*FM !! +(W0DTa-W0DTb)*F1
         W0XXX=W0XXXa*FP+W0XXXb*FM !! +3.d0*(W0DXXa-W0DXXb)*F1+3.d0*(W0DXa-W0DXb)*F2+(W0a-W0b)*F3
         W0XTT=W0XTTa*FP+W0XTTb*FM !! +(W0DTTa-W0DTTb)*F1
         W0XXT=W0XXTa*FP+W0XXTb*FM !! +2.d0*(W0DXTa-W0DXTb)*F1+(W0DTa-W0DTb)*F2
         W1=W1a*FP+W1b*FM
         W1DX=W1DXa*FP+W1DXb*FM !! +(W1a-W1b)*F1
         W1DT=W1DTa*FP+W1DTb*FM
         W1DXX=W1DXXa*FP+W1DXXb*FM !! +2.d0*(W1DXa-W1DXb)*F1+(W1a-W1b)*F2
         W1DTT=W1DTTa*FP+W1DTTb*FM
         W1DXT=W1DXTa*FP+W1DXTb*FM !! +(W1DTa-W1DTb)*F1
         W2=W2a*FP+W2b*FM
         W2DX=W2DXa*FP+W2DXb*FM !! +(W2a-W2b)*F1
         W2DT=W2DTa*FP+W2DTb*FM
         W2DXX=W2DXXa*FP+W2DXXb*FM !! +2.d0*(W2DXa-W2DXb)*F1+(W2a-W2b)*F2
         W2DTT=W2DTTa*FP+W2DTTb*FM
         W2DXT=W2DXTa*FP+W2DXTb*FM !! 
      else
         call BLIN9c(TEMP,CHI,
     *     W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT,
     *     W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,
     *     W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,
     *     W0XXX,W0XTT,W0XXT)
      endif
      return
      end

      subroutine BLIN9a(TEMP,CHI,
     *  W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT,
     *  W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,
     *  W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,
     *  W0XXX,W0XTT,W0XXT)
*                                                       Version 19.01.10
* First part of BILN9: small CHI. Stems from BLIN9 v.24.12.08
      implicit double precision (A-H), double precision (O-Z)
      save
      dimension AC(5,0:2),AU(5,0:2),AA(5,0:2)
      data AC/.37045057 d0, .41258437 d0,
     &        9.777982 d-2, 5.3734153 d-3, 3.8746281 d-5, ! c_i^0
     &        .39603109 d0, .69468795 d0, 
     &        .22322760 d0, 1.5262934 d-2, 1.3081939 d-4, ! c_i^1
     &        .76934619 d0, 1.7891437 d0, 
     &        .70754974 d0, 5.6755672 d-2, 5.5571480 d-4/ ! c_i^2
      data AU/.43139881 d0, 1.7597537 d0, 
     &        4.1044654 d0, 7.7467038 d0, 13.457678 d0, ! \chi_i^0
     &        .81763176 d0, 2.4723339 d0, 
     &        5.1160061 d0, 9.0441465 d0, 15.049882 d0, ! \chi_i^1
     &        1.2558461 d0, 3.2070406 d0, 
     &        6.1239082 d0, 10.316126 d0, 16.597079 d0/ ! \chi_i^2
      data KRUN/0/
      KRUN=KRUN+1
      if (KRUN.eq.1) then ! initialize
        do J=0,2
        do I=1,5
           AA(I,J)=dexp(-AU(I,J))
        enddo
        enddo
      endif
        do K=0,2
           W=0.
           WDX=0.
           WDT=0.
           WDXX=0.
           WDTT=0.
           WDXT=0.
           WDXXX=0.
           WDXTT=0.
           WDXXT=0.
             ECHI=dexp(-CHI)
            do I=1,5
               SQ=dsqrt(1.d0+AU(I,K)*TEMP/2.)
               DN=AA(I,K)+ECHI ! e^{-\chi_i}+e^{-\chi})
               W=W+AC(I,K)*SQ/DN
               WDX=WDX+AC(I,K)*SQ/DN**2
               WDT=WDT+AC(I,K)*AU(I,K)/(SQ*DN)
               WDXX=WDXX+AC(I,K)*SQ*(ECHI-AA(I,K))/DN**3
               WDTT=WDTT-AC(I,K)*AU(I,K)**2/(DN*SQ**3)
               WDXT=WDXT+AC(I,K)*AU(I,K)/(SQ*DN**2)
               WDXXX=WDXXX+AC(I,K)*SQ*
     *           (ECHI**2-4.*ECHI*AA(I,K)+AA(I,K)**2)/DN**4
               WDXTT=WDXTT-AC(I,K)*AU(I,K)**2/(DN**2*SQ**3)
               WDXXT=WDXXT+AC(I,K)*AU(I,K)*(ECHI-AA(I,K))/(SQ*DN**3)
            enddo
             WDX=WDX*ECHI
             WDT=WDT/4.
             WDXX=WDXX*ECHI
             WDTT=WDTT/16.
             WDXT=WDXT/4.*ECHI
             WDXXX=WDXXX*ECHI
             WDXTT=WDXTT*ECHI/16.
             WDXXT=WDXXT/4.*ECHI
          if (K.eq.0) then
             W0=W
             W0DX=WDX
             W0DT=WDT
             W0DXX=WDXX
             W0DTT=WDTT
             W0DXT=WDXT
             W0XXX=WDXXX
             W0XTT=WDXTT
             W0XXT=WDXXT
          elseif (K.eq.1) then
             W1=W
             W1DX=WDX
             W1DT=WDT
             W1DXX=WDXX
             W1DTT=WDTT
             W1DXT=WDXT
          else
             W2=W
             W2DX=WDX
             W2DT=WDT
             W2DXX=WDXX
             W2DTT=WDTT
             W2DXT=WDXT
          endif
        enddo ! next K
      return
      end

      subroutine BLIN9b(TEMP,CHI,
     *  W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT,
     *  W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,
     *  W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,
     *  W0XXX,W0XTT,W0XXT)
*                                                       Version 19.01.10
* Second part of BILN9: intermediate CHI. Stems from BLIN8 v.24.12.08
      implicit double precision (A-H), double precision (O-Z)
      save
      dimension AX(5),AXI(5),AH(5),AV(5)
      data AX/7.265351 d-2, .2694608 d0, 
     &        .533122 d0, .7868801 d0, .9569313 d0/ ! x_i
      data AXI/.26356032 d0, 1.4134031 d0, 
     &         3.5964258 d0, 7.0858100 d0, 12.640801 d0/ ! \xi_i
      data AH/3.818735 d-2, .1256732 d0, 
     &        .1986308 d0, .1976334 d0, .1065420 d0/ ! H_i
      data AV/.29505869 d0, .32064856 d0, 7.3915570 d-2, 
     &        3.6087389 d-3, 2.3369894 d-5/ ! \bar{V}_i
      parameter (EPS=1.d-3)
      if (CHI.lt.EPS) stop'BLIN9b: CHI is too small'
        do K=0,2
           W=0.
           WDX=0.
           WDT=0.
           WDXX=0.
           WDTT=0.
           WDXT=0.
           WDXXX=0.
           WDXTT=0.
           WDXXT=0.
             SQCHI=dsqrt(CHI)
            do I=1,5
               CE=AX(I)-1.d0
               ECHI=dexp(CE*CHI)
               DE=1.d0+ECHI
               D=1.d0+AX(I)*CHI*TEMP/2.
               H=CHI**(K+1)*SQCHI*dsqrt(D)/DE
               HX=(K+1.5)/CHI+.25*AX(I)*TEMP/D-ECHI*CE/DE
               HDX=H*HX
               HXX=(K+1.5)/CHI**2+.125*(AX(I)*TEMP/D)**2+ECHI*(CE/DE)**2
               HDXX=HDX*HX-H*HXX
               HT=.25*AX(I)*CHI/D
               HDT=H*HT
               HDTT=-H*HT**2
               HTX=1./CHI-.5*AX(I)*TEMP/D
               HDXT=HDX*HT+HDT*HTX
               HDXXT=HDXX*HT+HDX*HT*HTX+HDXT*HTX+
     +           HDT*(.25*(AX(I)*TEMP/D)**2-1./CHI**2)
               HDXTT=HDXT*HT-HDX*.125*(AX(I)*CHI/D)**2+HDTT*HTX+
     +           HDT*.5*AX(I)*(TEMP*.5*AX(I)*CHI/D**2-1./D)
               HXXX=(2*K+3)/CHI**3+.125*(AX(I)*TEMP/D)**3-
     -           ECHI*(1.d0-ECHI)*(CE/DE)**3
               HDXXX=HDXX*HX-2.*HDX*HXX+H*HXXX
               XICHI=AXI(I)+CHI
               DXI=1.d0+XICHI*TEMP/2.
               V=XICHI**K*dsqrt(XICHI*DXI)
               VX=(K+.5)/XICHI+.25*TEMP/DXI
               VDX=V*VX
               VT=.25*XICHI/DXI
               VDT=V*VT
               VXX=(K+.5)/XICHI**2+.125*(TEMP/DXI)**2
               VDXX=VDX*VX-V*VXX
               VDXXX=VDXX*VX-2.*VDX*VXX+
     +           V*((2*K+1)/XICHI**3+.125*(TEMP/DXI)**3)
               VXXT=(1.-.5*TEMP*XICHI/DXI)/DXI
               VDTT=-V*VT**2
               VXT=1./XICHI-.5*TEMP/DXI
               VDXT=VDT*VXT+VDX*VT
               VDXXT=VDXT*VX+VDX*.25*VXXT-VDT*VXX-V*.25*TEMP/DXI*VXXT
               VDXTT=VDTT*VXT-VDT*.5*VXXT+VDXT*VT-
     -           VDX*.125*(XICHI/DXI)**2
               W=W+AH(I)*AX(I)**K*H+AV(I)*V
               WDX=WDX+AH(I)*AX(I)**K*HDX+AV(I)*VDX
               WDT=WDT+AH(I)*AX(I)**K*HDT+AV(I)*VDT
               WDXX=WDXX+AH(I)*AX(I)**K*HDXX+AV(I)*VDXX
               WDTT=WDTT+AH(I)*AX(I)**K*HDTT+AV(I)*VDTT
               WDXT=WDXT+AH(I)*AX(I)**K*HDXT+AV(I)*VDXT
               WDXXX=WDXXX+AH(I)*AX(I)**K*HDXXX+AV(I)*VDXXX
               WDXTT=WDXTT+AH(I)*AX(I)**K*HDXTT+AV(I)*VDXTT
               WDXXT=WDXXT+AH(I)*AX(I)**K*HDXXT+AV(I)*VDXXT
            enddo
          if (K.eq.0) then
             W0=W
             W0DX=WDX
             W0DT=WDT
             W0DXX=WDXX
             W0DTT=WDTT
             W0DXT=WDXT
             W0XXX=WDXXX
             W0XTT=WDXTT
             W0XXT=WDXXT
          elseif (K.eq.1) then
             W1=W
             W1DX=WDX
             W1DT=WDT
             W1DXX=WDXX
             W1DTT=WDTT
             W1DXT=WDXT
          else
             W2=W
             W2DX=WDX
             W2DT=WDT
             W2DXX=WDXX
             W2DTT=WDTT
             W2DXT=WDXT
          endif
        enddo ! next K
      return
      end

      subroutine BLIN9c(TEMP,CHI,
     *  W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT,
     *  W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,
     *  W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,
     *  W0XXX,W0XTT,W0XXT)
*                                                       Version 19.01.10
* Third part of BILN9: large CHI. Stems from BLIN8 v.24.12.08
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter (PI=3.141592653d0,PI26=PI*PI/6.)
      dimension AM(0:2),AMDX(0:2),AMDT(0:2),
     *  AMDXX(0:2),AMDTT(0:2),AMDXT(0:2)
      if (CHI*TEMP.lt..1) then
        do K=0,2
           W=0.
           WDX=0.
           WDT=0.
           WDXX=0.
           WDTT=0.
           WDXT=0.
           WDXXX=0.
           WDXTT=0.
           WDXXT=0.
          do J=0,4 ! for nonrel.Fermi integrals from k+1/2 to k+4.5
             CNU=K+J+.5 ! nonrelativistic Fermi integral index \nu
             CHINU=CHI**(K+J)*dsqrt(CHI) ! \chi^\nu
             F=CHINU*(CHI/(CNU+1.)+PI26*CNU/CHI+ ! nonrel.Fermi
     +         .7*PI26**2*CNU*(CNU-1.)*(CNU-2.)/CHI**3)
             FDX=CHINU*(1.+PI26*CNU*(CNU-1.)/CHI**2+
     +         .7*PI26**2*CNU*(CNU-1.)*(CNU-2.)*(CNU-3.)/CHI**4)
             FDXX=CHINU/CHI*CNU*(1.+PI26*(CNU-1.)*(CNU-2.)/CHI**2+
     +         .7*PI26**2*(CNU-1.)*(CNU-2.)*(CNU-3.)*(CNU-4.)/CHI**4)
             FDXXX=CHINU/CHI**2*CNU*(CNU-1.)*
     *         (1.+PI26*(CNU-2.)*(CNU-3.)/CHI**2+
     +         .7*PI26**2*(CNU-2.)*(CNU-3.)*(CNU-4.)*(CNU-5.)/CHI**4)
            if (J.eq.0) then
               W=F
               WDX=FDX
               WDXX=FDXX
               WDXXX=FDXXX
            elseif (J.eq.1) then
               C=.25*TEMP
               W=W+C*F ! Fermi-Dirac, expressed through Fermi
               WDX=WDX+C*FDX
               WDXX=WDXX+C*FDXX
               WDT=F/4.
               WDXT=FDX/4.
               WDTT=0.
               WDXXX=WDXXX+C*FDXXX
               WDXXT=FDXX/4.
               WDXTT=0.
            else
               C=-C/J*(2*J-3)/4.*TEMP
               W=W+C*F
               WDX=WDX+C*FDX
               WDT=WDT+C*J/TEMP*F
               WDXX=WDXX+C*FDXX
               WDTT=WDTT+C*J*(J-1)/TEMP**2*F
               WDXT=WDXT+C*J/TEMP*FDX
               WDXXX=WDXXX+C*FDXXX
               WDXTT=WDXTT+C*J*(J-1)/TEMP**2*FDX
               WDXXT=WDXXT+C*J/TEMP*FDXX
            endif
          enddo ! next J
          if (K.eq.0) then
             W0=W
             W0DX=WDX
             W0DT=WDT
             W0DXX=WDXX
             W0DTT=WDTT
             W0DXT=WDXT
             W0XXX=WDXXX
             W0XTT=WDXTT
             W0XXT=WDXXT
          elseif (K.eq.1) then
             W1=W
             W1DX=WDX
             W1DT=WDT
             W1DXX=WDXX
             W1DTT=WDTT
             W1DXT=WDXT
          else
             W2=W
             W2DX=WDX
             W2DT=WDT
             W2DXX=WDXX
             W2DTT=WDTT
             W2DXT=WDXT
          endif
        enddo ! next K
*   ----------------------------------------------------------------   *
      else ! CHI > 14, CHI*TEMP > 0.1: general high-\chi expansion
         D=1.d0+CHI*TEMP/2.d0
         R=dsqrt(CHI*D)
         RX=.5d0/CHI+.25d0*TEMP/D
         RDX=R*RX
         RDT=.25d0*CHI**2/R
         RXX=-.5d0/CHI**2-.125d0*(TEMP/D)**2
         RDXX=RDX*RX+R*RXX
         RDTT=-.25d0*RDT*CHI/D
         RXT=.25d0/D-.125d0*CHI*TEMP/D**2
         RDXT=RDT*RX+R*RXT
         RXXX=1.d0/CHI**3+.125d0*(TEMP/D)**3
         RDXXX=RDXX*RX+2.d0*RDX*RXX+R*RXXX
         RXTT=-.25d0/D**2*CHI+.125d0*CHI**2*TEMP/D**3
         RDXTT=RDTT*RX+2.d0*RDT*RXT+R*RXTT
         RXXT=-RXT*TEMP/D
         RDXXT=RDXT*RX+RDX*RXT+RDT*RXX+R*RXXT
        do K=0,2
           DM=K+.5d0+(K+1.d0)*CHI*TEMP/2.d0
           AM(K)=CHI**K*DM/R
           FMX1=.5d0*(K+1.)*TEMP/DM
           FMX2=.25d0*TEMP/D
           FMX=(K-.5d0)/CHI+FMX1-FMX2
           AMDX(K)=AM(K)*FMX
           CKM=.5d0*(K+1.d0)/DM
           FMT1=CKM*CHI
           FMT2=.25d0*CHI/D
           FMT=FMT1-FMT2
           AMDT(K)=AM(K)*FMT
           FMXX=-(K-.5d0)/CHI**2-FMX1**2+2.d0*FMX2**2
           AMDXX(K)=AMDX(K)*FMX+AM(K)*FMXX
           FMTT=2.d0*FMT2**2-FMT1**2
           AMDTT(K)=AMDT(K)*FMT+AM(K)*FMTT
           AMDXT(K)=AMDX(K)*FMT+AM(K)*(CKM*(1.d0-CKM*CHI*TEMP)-
     -       .25d0/D+.125d0*CHI*TEMP/D**2)
          if (K.eq.0) then
             FMXXX=(2*K-1)/CHI**3+2.d0*FMX1**3-8.d0*FMX2**3
             AMDXXX=AMDXX(K)*FMX+2.d0*AMDX(K)*FMXX+AM(K)*FMXXX
             FMT1DX=CKM-TEMP*CHI*CKM**2
             FMT2DX=(.25d0-CHI*TEMP*.125d0/D)/D
             FMXT=FMT1DX-FMT2DX
             FMTTX=4.d0*FMT2*FMT2DX-2.d0*FMT1*FMT1DX
             AMDXTT=AMDXT(K)*FMT+AMDT(K)*FMXT+AMDX(K)*FMTT+AM(K)*FMTTX
             FMX1DT=CKM-CHI*TEMP*CKM**2
             FMX2DT=.25d0/D*(1.d0-.5d0*CHI*TEMP/D)
             FMXXT=4.d0*FMX2*FMX2DT-2.d0*FMX1*FMX1DT
             AMDXXT=AMDXT(K)*FMX+AMDX(K)*FMXT+AMDT(K)*FMXX+AM(K)*FMXXT
          endif
        enddo
         SQ2T=dsqrt(2.d0*TEMP)
           A=1.d0+CHI*TEMP+SQ2T*R
           ADX=TEMP+SQ2T*RDX
           ADT=CHI+R/SQ2T+SQ2T*RDT
           ADXX=SQ2T*RDXX
           ADTT=-R/SQ2T**3+2.d0/SQ2T*RDT+SQ2T*RDTT
           ADXT=1.d0+RDX/SQ2T+SQ2T*RDXT
           ADXTT=-RDX/SQ2T**3+2.d0/SQ2T*RDXT+SQ2T*RDXTT
           ADXXT=RDXX/SQ2T+SQ2T*RDXXT
           XT1=CHI+1.d0/TEMP
           Aln=dlog(A)
           FJ0=.5d0*XT1*R-Aln/SQ2T**3
           ASQ3=A*SQ2T**3
           ASQ3DX=ADX*SQ2T**3
           FJ0DX=.5d0*(R+XT1*RDX)-ADX/ASQ3
           FJ0DT=.5d0*(XT1*RDT-R/TEMP**2)-ADT/ASQ3+
     +       .75d0/(TEMP**2*SQ2T)*Aln
           FJ0DXX=RDX+.5d0*XT1*RDXX+(ADX/A)**2/SQ2T**3-ADXX/ASQ3
           FJ0DTT=R/TEMP**3-RDT/TEMP**2+.5d0*XT1*RDTT+
     +       3.d0/(ASQ3*TEMP)*ADT+
     +     (ADT/A)**2/SQ2T**3-ADTT/ASQ3-1.875d0/(TEMP**3*SQ2T)*Aln
           BXT=1.5d0/TEMP*ADX+ADX*ADT/A-ADXT
           BXXT=1.5d0/TEMP*ADXX+(ADXX*ADT+ADX*ADXT)/A-
     -       (ADX/A)**2*ADT-ADXXT
           FJ0DXT=.5d0*(RDT-RDX/TEMP**2+XT1*RDXT)+BXT/ASQ3
           FJ0XXX=RDXX*1.5d0+.5d0*XT1*RDXXX+
     +      (2.d0*ADX*(ADXX/A-(ADX/A)**2)-
     -      SQ2T*RDXXX+ADXX/ASQ3*ASQ3DX)/ASQ3
           FJ0XTT=RDX/TEMP**3-RDXT/TEMP**2+.5d0*(RDTT+XT1*RDXTT)+
     +      3.d0/TEMP*(ADXT-ADT/ASQ3*ASQ3DX)/ASQ3+
     +      (2.d0*ADT*(ADXT/A-ADT*ADX/A**2)-
     -      ADXTT+ADTT*ASQ3DX/ASQ3)/ASQ3-1.875d0/(TEMP**3*SQ2T)*ADX/A
           FJ0XXT=.5d0*(RDXT-RDXX/TEMP**2+RDXT+XT1*RDXXT)+
     +      (BXXT-BXT*ASQ3DX/ASQ3)/ASQ3
         W0=FJ0+PI26*AM(0)
         W0DX=FJ0DX+PI26*AMDX(0)
         W0DT=FJ0DT+PI26*AMDT(0)
         W0DXX=FJ0DXX+PI26*AMDXX(0)
         W0DTT=FJ0DTT+PI26*AMDTT(0)
         W0DXT=FJ0DXT+PI26*AMDXT(0)
         W0XXX=FJ0XXX+PI26*AMDXXX
         W0XTT=FJ0XTT+PI26*AMDXTT
         W0XXT=FJ0XXT+PI26*AMDXXT
           FJ1=(R**3/1.5d0-FJ0)/TEMP
           FJ1DX=(2.d0*R**2*RDX-FJ0DX)/TEMP
           FJ1DT=(2.d0*R**2*RDT-FJ0DT-FJ1)/TEMP
           FJ1DXX=(4.d0*R*RDX**2+2.d0*R**2*RDXX-FJ0DXX)/TEMP
           FJ1DTT=(4.d0*R*RDT**2+2.d0*R**2*RDTT-FJ0DTT-2.d0*FJ1DT)/TEMP
           FJ1DXT=(4.d0*R*RDX*RDT+2.d0*R**2*RDXT-FJ0DXT-FJ1DX)/TEMP
         W1=FJ1+PI26*AM(1)
         W1DX=FJ1DX+PI26*AMDX(1)
         W1DT=FJ1DT+PI26*AMDT(1)
         W1DXX=FJ1DXX+PI26*AMDXX(1)
         W1DTT=FJ1DTT+PI26*AMDTT(1)
         W1DXT=FJ1DXT+PI26*AMDXT(1)
           FJ2=(.5d0*CHI*R**3-1.25d0*FJ1)/TEMP
           FJ2DX=(.5d0*R**3+1.5d0*CHI*R**2*RDX-1.25d0*FJ1DX)/TEMP
           FJ2DT=(1.5d0*CHI*R**2*RDT-1.25d0*FJ1DT-FJ2)/TEMP
           FJ2DXX=(3.d0*R*RDX*(R+CHI*RDX)+1.5d0*CHI*R**2*RDXX-
     -       1.25d0*FJ1DXX)/TEMP
          FJ2DTT=(3.d0*CHI*R*(RDT**2+.5d0*R*RDTT)-
     -      1.25d0*FJ1DTT-2.d0*FJ2DT)/TEMP
           FJ2DXT=(1.5d0*R*RDT*(R+2.d0*CHI*RDX)+1.5d0*CHI*R**2*RDXT-
     -       1.25d0*FJ1DXT-FJ2DX)/TEMP
         W2=FJ2+PI26*AM(2)
         W2DX=FJ2DX+PI26*AMDX(2)
         W2DT=FJ2DT+PI26*AMDT(2)
         W2DXX=FJ2DXX+PI26*AMDXX(2)
         W2DTT=FJ2DTT+PI26*AMDTT(2)
         W2DXT=FJ2DXT+PI26*AMDXT(2)
      endif
      return
      end

      subroutine CHEMFIT(DENS,TEMP,CHI)
*                                                       Version 07.06.07
* This is merely an interface to CHEMFIT7 for compatibility purposes.
* Input:  DENS - electron density [a.u.=6.7483346e24 cm^{-3}],
* TEMP - temperature [a.u.=2Ryd=3.1577e5 K]
* Output: CHI=\mu/TEMP, where \mu - electron chem.pot.w/o rest-energy
      implicit double precision (A-H), double precision (O-Z)
      save
      DENR=DENS/2.5733806d6 ! n_e in rel.un.=\lambda_{Compton}^{-3}
      TEMR=TEMP/1.8778865d4 ! T in rel.un.=(mc^2/k)=5.93e9 K
      call CHEMFIT7(DENR,TEMR,CHI,CMU1,0,CMUDENR,CMUDT,CMUDTT)
      return
      end

      subroutine CHEMFIT7(DENR,TEMR,CHI,CMU1,KDERIV,
     *  CMUDENR,CMUDT,CMUDTT)
*                                                       Version 28.05.07
*                                                     corrected 08.08.11
* Fit to the chemical potential of free electron gas described in:
*     G.Chabrier & A.Y.Potekhin, Phys.Rev.E 58, 4941 (1998)
* Stems from CHEMFIT v.10.10.96. The main difference - derivatives.
*  All quantities are by default in relativistic units
* Input:  DENR - electron density, TEMR - temperature
*         KDERIV=0 if the derivatives are not required
* Output: CHI=CMU1/TEMR, where CMU1 = \mu-1 - chem.pot.w/o rest-energy
*         CMUDENR = (d\mu/d n_e)_T
*         CMUDT = (d\mu/dT)_V
*         CMUDTT = (d^2\mu/dT^2)_V
* CMUDENR,CMUDT, and CMUDTT =0 on output, if KREDIV=0
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter (PARA=1.612d0,PARB=6.192d0,PARC=.0944d0,
     *  PARF=5.535d0,PARG=.698d0)
      PF0=(29.6088132d0*DENR)**.33333333d0 ! Classical Fermi momentum
      if (PF0.gt.1.d-4) then
         TF=dsqrt(1.d0+PF0**2)-1.d0 ! Fermi temperature
      else
         TF=.5d0*PF0**2
      endif
      THETA=TEMR/TF
      THETA32=THETA*dsqrt(THETA)
      Q2=12.d0+8.d0/THETA32
      T1=dexp(-THETA) ! former ('96) 1/T
      U3=T1**2+PARA
      THETAC=THETA**PARC
      THETAG=THETA**PARG
      D3=PARB*THETAC*T1**2+PARF*THETAG
      Q3=1.365568127d0-U3/D3 ! 1.365...=2/\pi^{1/3}
      if (THETA.gt.1.d-5) then 
         Q1=1.5d0*T1/(1.d0-T1)
      else
         Q1=1.5d0/THETA
      endif
      SQT=dsqrt(TEMR)
      G=(1.d0+Q2*TEMR*Q3+Q1*SQT)*TEMR
      H=(1.d0+.5d0*TEMR/THETA)*(1.d0+Q2*TEMR)
      CT=1.d0+G/H
      F=.666666667d0/THETA32
      call FERINV7(F,1,X,XDF,XDFF)
      CHI=X ! non-relativistic result
     -  -    1.5d0*dlog(CT) ! Relativistic fit
      CMU1=TEMR*CHI ! Fit to chemical potential w/o mc^2
      if (KDERIV.eq.0) then ! DISMISS DERIVATIVES
         CMUDENR=0.
         CMUDT=0.
         CMUDTT=0.
         return
      endif
* CALCULATE DERIVATIVES:
* 1: derivatives of CHI over THETA and T
* (a): Non-relativistic result:
      THETA52=THETA32*THETA
      CHIDY=-XDF/THETA52 ! d\chi/d\theta
      CHIDYY=(XDFF/THETA**4-2.5d0*CHIDY)/THETA ! d^2\chi/d\theta^2
* (b): Relativistic corrections:
      if (THETA.gt.1.d-5) then 
         Q1D=-Q1/(1.d0-T1)
         Q1DD=-Q1D*(1.d0+T1)/(1.d0-T1)
      else
         Q1D=-1.5d0/THETA**2
         Q1DD=-2.d0*Q1D/THETA ! sign corrected 08.08.11
      endif
      Q2D=-12.d0/THETA52 ! d q_2 / d \theta
      Q2DD=30.d0/(THETA52*THETA) ! d^2 q_2 / d \theta^2
      U3D=-2.d0*T1**2
      D3D=PARF*PARG*THETAG/THETA+PARB*T1**2*THETAC*(PARC/THETA-2.d0)
      D3DD=PARF*PARG*(PARG-1.d0)*THETAG/THETA**2+
     +PARB*T1**2*THETAC*(PARC*(PARC-1.d0)/THETA**2-4.d0*PARC/THETA+4.d0)
      Q3D=(D3D*U3/D3-U3D)/D3
      Q3DD=(2.d0*U3D+(2.d0*U3D*D3D+U3*D3DD)/D3-2.d0*U3*(D3D/D3)**2)/D3
      GDY=TEMR*(Q1D*SQT+(Q2D*Q3+Q2*Q3D)*TEMR) ! dG/d\theta
      GDT=1.d0+1.5d0*Q1*SQT+2.d0*Q2*Q3*TEMR
      GDYY=TEMR*(Q1DD*SQT+(Q2DD*Q3+2.d0*Q2D*Q3D+Q2*Q3DD)*TEMR)
      GDTT=.75d0*Q1/SQT+2.d0*Q2*Q3
      GDYT=1.5d0*Q1D*SQT+2.d0*(Q2D*Q3+Q2*Q3D)*TEMR
      HDY=(-.5d0/THETA**2+Q2D+.5d0*(Q2D-Q2/THETA)/THETA*TEMR)*TEMR
      HDT=(.5d0+Q2*TEMR)/THETA+Q2
      HDYY=TEMR/THETA**3+Q2DD*TEMR+
     +  TEMR**2*(.5d0*Q2DD-Q2D/THETA+Q2/THETA**2)/THETA
      HDTT=Q2/THETA
      HDYT=Q2D*(1.d0+TEMR/THETA)-(.5d0+Q2*TEMR)/THETA**2
      CTY=GDY/G-HDY/H
      CTT=GDT/G-HDT/H
      GH=G/H
      CTDY=GH*CTY
      CTDT=GH*CTT
      CTDYY=CTDY*CTY+GH*(GDYY/G-(GDY/G)**2-HDYY/H+(HDY/H)**2)
      CTDTT=CTDT*CTT+GH*(GDTT/G-(GDT/G)**2-HDTT/H+(HDT/H)**2)
      CTDYT=CTDT*CTY+GH*(GDYT/G-GDY*GDT/G**2-HDYT/H+HDY*HDT/H**2)
      CHIDY=CHIDY-1.5d0*CTDY/CT
      CHIDT=-1.5d0*CTDT/CT
      CHIDYY=CHIDYY+1.5d0*((CTDY/CT)**2-CTDYY/CT)
      CHIDTT=1.5d0*((CTDT/CT)**2-CTDTT/CT)
      CHIDYT=1.5d0*(CTDY*CTDT/CT**2-CTDYT/CT)
      CMUDENR=-(THETA*PF0)**2/(3.d0*DENR*(1.d0+TF))*CHIDY
      CMUDT=CHI+THETA*CHIDY+TEMR*CHIDT
      CMUDTT=2.d0*(CHIDY/TF+CHIDT+THETA*CHIDYT)+
     +  THETA/TF*CHIDYY+TEMR*CHIDTT
      return
      end
