*** SUPPLEMENT TO eos12.f FOR THE CASE OF ARBITRARY MAGNETIC FIELD
*
*             EQUATION OF STATE FOR ELECTRON-ION PLASMAS
*
*                 USES SUBROUTINES FROM eos12.f
*               (!!! eos12.f needs to be linked !!!)
*
* A.Y.Potekhin & G.Chabrier, in preparation (2012)
* Please send comments/suggestions to Alexander Potekhin:
*                                                   palex@astro.ioffe.ru
* Last updated 02.10.12
**   L I S T   O F   S U B R O U T I N E S (functions go with "="):
* MAIN (normally commented-out) - example driving routine
* EOSMAG - EOS of a homogeneous plasma of electrons and one ion species
* CHEMAG8 - EOS of magnetized electron ideal gas (density as input)
* ELECTMAG - EOS of magnetized electron ideal gas (chem.pot. as input)
* MAGNION - EOS of nonrelativistic nondegenerate ion gas in magn.field.
* SPINION - contributions to the EOS of nonrelativistic nondegenerate
*           ion gas due to the spin degeneracy and magnetic moments.
* EOSFIM12 - non-ideal parts of thermodynamic functions in magn.field.
* EXCORM - exchange-correlation contribution for the electron gas
* FHARMAG - 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.
* HLMAG - the same as FHARM12, but only for thermal contributions
* darcth=hyperbolic arc-tangent (inverse dtanh)
************************************************************************
*                           MAIN program:               Version 26.05.12
* 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 (UN_B12=425.438,UN_T6=.3157746)
C%C      data LIQSOL/0/
C%C      write(*,'('' Input Z, A: ''$)')
C%C      read*,Zion,CMI
C%C      write(*,110) Zion,CMI
C%C      write(*,'('' Input B12, T6, and RHO: ''$)')
C%C      read*,B12,T6,RHO
C%C    1 continue
C%C      write(*,111) B12,T6
C%C      TEMP=T6/UN_T6 ! T [au]
C%C      GAMAG=B12*UN_B12
C%C    2 continue
C%C      call EOSMAG(Zion,CMI,RHO,TEMP,GAMAG,
C%C     *  DENS,GAMI,CHI,TPT,LIQSOL,PnkT,UNkT,SNk,CV,CHIR,CHIT)
C%C      Tnk=8.31447d13/CMI*RHO*T6 ! n_i kT [erg/cc]
C%C      P=PnkT*Tnk/1.d12 ! P [Mbar]
C%C*   --------------------   OUTPUT   --------------------------------   *
C%C* Here in the output we have:
C%C* RHO - mass density in g/cc
C%C* P - total pressure in Mbar (i.e. in 1.e12 dyn/cm^2)
C%C* PnkT=P/nkT, where n is the number density of ions, T temperature
C%C* CV - heat capacity at constant volume, divided by number of ions, /k
C%C* CHIT - logarithmic derivative of pressure \chi_T
C%C* CHIR - logarithmic derivative of pressure \chi_\rho
C%C* UNkT - internal energy divided by NkT, N being the number of ions
C%C* SNk - entropy divided by number of ions, /k
C%C* GAMI - ionic Coulomb coupling parameter
C%C* TPT=T_p/T, where T_p is the ion plasma temperature
C%C* CHI - electron chemical potential, divided by kT
C%C* LIQSOL = 0 in the liquid state, = 1 in the solid state
C%C      write(*,111) RHO,P,PnkT,CV,CHIT,CHIR,UNkT,SNk,GAMI,TPT,CHI,
C%C     *  LIQSOL
C%C      write(*,'('' Next RHO ( < 0 to new B and T): ''$)')
C%C      read*,RHO
C%C      if (RHO.gt.0.) goto 2
C%C      write(*,'('' Next B12,T6 ( < 0 to stop): ''$)')
C%C      read*,B12,T6
C%C      if (T6.le.0..or.B12.lt.0.) stop
C%C      write(*,'('' RHO: ''$)')
C%C      read*,RHO
C%C      goto 1
C%C  110 format(' Output of eosmag v.27.05.12'/
C%C     *  '   Zion      CMI'/1P,2E11.4/
C%C     *  ' rho [g/cc]  P [Mbar]  P/(n_i kT)  Cv/N_i',
C%C     *  '     chi_T      chi_r    U/(N_i kT)    S/N_i     Gamma_i',
C%C     *  '   T_p/T   chi_e   liq/sol'/'   B12      T6')
C%C  111 format(1P,11E11.3,I2)
C%C      end

      subroutine EOSMAG(Zion,CMI,RHO,TEMP,GAMAG,
     *  DENS,GAMI,CHI,TPT,LIQSOL,PnkT,UNkT,SNk,CV,CHIR,CHIT)
*                                                       Version 06.07.12
*                                             slightly modified 13.09.12
* EOS of fully ionized electron-ion plasma with magnetic field
* Limitations: inapplicable in the regimes of (1) bound-state formation,
*      (2) quantum liquid, (3) presence of positrons, (4) mixtures
* Input: Zion and CMI - ion charge and mass numbers,
*        RHO - total mass density [g/cc]
*        TEMP - temperature [in a.u.=2Ryd=3.158e5 K]
*        GAMAG - magnetic field [in a.u. (=B/2.3505e9 G)]
*        LIQSOL - regulator of the solid vs liquid regime (see below)
* NB: instead of RHO, a true input is CHI, determined below in CHEMAG8.
*     Hence, disagreement between RHO and DENS is the CHEMAG8 error.
* NB': Photon-gas quantities (radiation pressure etc.) are NOT INCLUDED
*     in this subroutine (if needed, add them externally).
* Output: DENS - electron number density [in a.u.=6.7483346e24 cm^{-3}]
*         GAMI - ion-ion Coulomb coupling constant
*         CHI = mu_e/kT, where mu_e is the electron chem.potential
*         TPT - ionic quantum parameter (T_p/T)
*         LIQSOL - regulator of the solid vs liquid regime (see below)
*         SNk - total dimensionless entropy per 1 ion (assumes spin=1/2)
*         UNkT - internal energy per kT per ion
*         PnkT - pressure / n_i kT, where n_i is the ion number density
*         CV - heat capacity per ion, divided by the Boltzmann constant
*         CHIR - inverse compressibility -(d ln P / d ln V)_T ("\chi_r")
*         CHIT = (d ln P / d ln T)_V ("\chi_T")
* LIQSOL is the liquid/solid switch, which works as follows:
*     if LIQSOL=0 or 1 on the input, then:
*        if (either GAMI or 1/RS is below its critical value) then
*           liquid regime and LIQSOL=0 on the output
*        otherwise solid regime and LIQSOL=1 on the output
*     if LIQSOL=2 or 3 on the input, then:
*        if (LIQSOL=2) then liquid regime and LIQSOL=2 on the output
*        if (LIQSOL=3) then solid regime and LIQSOL=3 on the output
*     if LIQSOL=4 or 5 on the input, then
*         consider non-ideal free energy FC1:
*        if (FC1=min in liquid) then liquid regime, output LIQSOL=4
*        if (FC1=min in solid) then solid regime, output LIQSOL=5
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter(TINY=1.d-7,EPS=1.d-2)
      parameter (PI=3.141592653d0,AUM=1822.888d0, ! a.m.u./m_e
     *   CMN=0.6, ! effective neutron mass m_n* / m_n
     *   GAMIMELT=175., ! OCP value of Gamma_i for melting
     *   GAML=1.d3, ! the bigest allowed GAMI in liquid
     *   GAMS=1.d0, ! the smallest allowed GAMI in solid
     *   BOHR=137.036, ! inverse fine-str.const.(Bohr radius in rel.un.)
     *   RSIMELT=140.) ! ion density parameter of quantum melting
      parameter (PI2=PI**2,BOHR2=BOHR**2,BOHR3=BOHR2*BOHR)
      DENS=RHO/11.20587*Zion/CMI ! number density of electrons [au]
      B=GAMAG/BOHR**2 ! magnetic field in rel.un.
      TEMR=TEMP/BOHR2 ! T [rel.un.]
      Z13=Zion**.33333333d0
      DENR=DENS/BOHR3 ! number density of electrons [rel.un.]
      DENRI=DENR/Zion ! number density of ions (= of cells)
* Plasma parameters:
      RS=(.75d0/PI/DENS)**.33333333d0 ! r_s - electron density parameter
      RSI=RS*CMI*Zion**2*Z13*AUM ! R_S - ion density parameter
      GAME=1.d0/RS/TEMP ! electron Coulomb parameter Gamma_e
      GAMI=Zion**2*GAME/Z13 ! effective Gamma_i - ion Coulomb parameter
      if (LIQSOL.lt.0.or.LIQSOL.gt.5) stop'EOSMAG: invalid LIQSOL'
      if (LIQSOL.eq.0.or.LIQSOL.eq.1) then
        if (GAMI.lt.GAMIMELT.or.RSI.lt.RSIMELT) then
           LIQSOL=0 ! liquid regime
        else
           LIQSOL=1 ! solid regime
        endif
      endif
      TPT2=GAME**2/RS*3.d0/(AUM*CMI)*Zion
* 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.
      if ((LIQSOL.eq.0.or.LIQSOL.eq.2).and.TPT2.gt.18.)
     *   stop'EOSMAG: too strong quantum effects in liquid'
      TPT=dsqrt(TPT2) ! effective T_p/T - ion quantum parameter
* (1) ideal electron gas (including relativity, degeneracy, and m.f.):
      call CHEMAG8(B,TEMR,DENR,
     *  DENR1,CHI,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,DDENR)
      if (abs(DENR/DENR1-1.).gt.EPS) print*,'EOSMAG Warning: DENS1/DENS'
* NB: CHI can be used as true input instead of RHO or DENS. In that case
*     one should define density as "DENS=DENR1*BOHR3"
* (2) ideal ion gas (including relativity, degeneracy, and m.f.):
      GYRO=0.d0
      MULTI=1
      if (dabs(Zion-1.d0).lt.TINY) then ! assume protons
         GYRO=5.5857d0 ! = g_p
         MULTI=2
      endif
      if (LIQSOL.eq.0.or.LIQSOL.eq.2.or.LIQSOL.eq.4.or.LIQSOL.eq.5) then
         call MAGNION(B,TEMR,DENRI,Zion,CMI,GYRO,MULTI,
     *     FIid,PIid,UIid,SIid,CVI,CHITI,CHIRI)
        if (LIQSOL.eq.4.or.LIQSOL.eq.5) then ! remember
           FIidL=FIid
           PIidL=PIid
           UIidL=UIid
           SIidL=SIid
           CVIL=CVI
           CHITIL=CHITI
        endif
      endif
      if (LIQSOL.eq.1.or.LIQSOL.eq.3.or.LIQSOL.eq.4.or.LIQSOL.eq.5) then
         PIid=1.d0
         CHITI=1.d0
         CHIRI=1.d0
         ZETI=Zion/(CMI*AUM)*GAMAG/TEMP ! \hbar\omega_{ci}/kT
         call SPINION(GYRO,ZETI,MULTI,Fspin,Uspin,CVspin)
         FIid=1.5d0*dlog(TPT**2/GAMI)-1.323515d0+Fspin
* 1.3235=1+0.5*ln(6/pi); FIid = F_{id.ion gas}/(N_i kT)
* NB: for a mixture, this FIid must be modified (see FidION in EOSFI7)
         UIid=1.5d0+Uspin
         CVI=1.5d0+CVspin
         SIid=UIid-FIid
      endif
* --- corr.for zero-point Landau vibrations:
      UIZP=GAMAG/(2.d0*CMI*AUM*TEMP)
* (3) Coulomb+xc nonideal contributions:
      if (LIQSOL.eq.0.or.LIQSOL.eq.2) then ! liquid
         call EOSFIM12(0,CMI,Zion,RS,GAMI,TEMP,GAMAG,1,0,
     *     FC1,UC1,PC1,SC1,CV1,PDT1,PDR1,
     *     FC2,UC2,PC2,SC2,CV2,PDT2,PDR2)
      endif
      if (LIQSOL.eq.1.or.LIQSOL.eq.3) then ! solid
         call EOSFIM12(1,CMI,Zion,RS,GAMI,TEMP,GAMAG,1,0,
     *     FC1,UC1,PC1,SC1,CV1,PDT1,PDR1,
     *     FC2,UC2,PC2,SC2,CV2,PDT2,PDR2)
      endif
      if (LIQSOL.eq.4.or.LIQSOL.eq.5) then ! choose liquid or solid
         call EOSFIM12(0,CMI,Zion,RS,GAMI,TEMP,GAMAG,1,0,
     *     FC1L,UC1L,PC1L,SC1L,CV1L,PDT1L,PDR1L,
     *     FC2L,UC2L,PC2L,SC2L,CV2L,PDT2L,PDR2L)
         call EOSFIM12(1,CMI,Zion,RS,GAMI,TEMP,GAMAG,1,0,
     *     FC1,UC1,PC1,SC1,CV1,PDT1,PDR1,
     *     FC2,UC2,PC2,SC2,CV2,PDT2,PDR2)
        if (FC1L.lt.FC1) then ! switch to the liquid; recall it.
          if (GAMI.gt.GAML) stop'EOSMAG: too high GAMI in liquid'
           LIQSOL=4
           FIid=FIidL
           PIid=PIidL
           UIid=UIidL
           SIid=SIidL
           CVI=CVIL
           CHITI=CHITIL
           FC1=FC1L
           UC1=UC1L
           PC1=PC1L
           SC1=SC1L
           CV1=CV1L
           PDT1=PDT1L
           PDR1=PDR1L
           FC2=FC2L
           UC2=UC2L
           PC2=PC2L
           SC2=SC2L
           CV2=CV2L
           PDT2=PDT2L
           PDR2=PDR2L
        else
          if (GAMI.lt.GAMS) stop'EOSMAG: too low GAMI in solid'
           LIQSOL=5
        endif
      endif
* Calculate partial thermodynamic quantities and combine them together:
** Normalization factors:
      DTE=DENS*TEMP ! electron normalization [au]
      DENSI=DENS/Zion ! number density of all ions [au]
      DTI=DENSI*TEMP ! ion normalization [au]
** Electrons:
*** First-order TD functions:
      UINTE=UEid*DTE
      PRESSE=PEid*DTE ! P_e [a.u.]
      StotE=SEid*DENS ! electron entropy per unit volume [a.u.]
*** 2nd-order TD functions:
      CVtotE=CVE*DENS
      PDLTE=PRESSE*CHITE
      PDLRE=PRESSE*CHIRE
** Ions:
      if (LIQSOL.eq.0.or.LIQSOL.eq.2.or.LIQSOL.eq.4) then ! liquid
*** First-order TD functions:
         UINTI=DTI*(UIid+UC1) ! (i+ii+ie+ee)
         PRESSI=DTI*(PIid+PC1)
         StotI=DENSI*(SIid+SC1)
*** 2nd-order TD functions:
         CVtotI=DENSI*(CVI+CV1)
         PDLTI=DTI*(PIid*CHITI+PDT1)
         PDLRI=DTI*(PIid*CHIRI+PDR1)
      else ! solid
*** First-order TD functions:
         UINTI=DTI*UC2
         PRESSI=DTI*PC2
         StotI=DENSI*SC2
*** 2nd-order TD functions:
         CVtotI=DENSI*CV2
         PDLTI=DTI*PDT2
         PDLRI=DTI*PDR2
      endif
** Final sum-up:
      UINT=UINTE+UINTI ! int.en. per unit V (e+i+Coul.) [a.u.]
      PRESS=PRESSE+PRESSI ! pressure (e+i+Coul.) [a.u.]
      Stot=StotE+StotI
      CVtot=CVtotE+CVtotI
      PDLT=PDLTE+PDLTI
      PDLR=PDLRE+PDLRI
* Dimensionless ratios:
** First-order:
      PnkT=PRESS/DTI ! beta P / n_i
      UNkT=UINT/DTI ! beta U / N_i
      SNk=Stot/DENSI ! S / N_i k_B
** 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
      return
      end

*   =========  CHEMICAL POTENTIAL and DENSITY DERIVATIVE  ==========   *
      subroutine CHEMAG8(B,TEMR,DENR,
     *  DENR1,CHI,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,DDENR)
      implicit double precision (A-H), double precision (O-Z)
      save
*                                                       Version 04.03.09
*                                               Cosmetic change 02.09.11
*                                                EPS diminished 26.05.12
* Stems from CHEMPOT v.03.05.07
* Ideal electron-gas EOS
* Input: B - magnetic field [rel.un.] (=B/4.414e13 G)
*        TEMR - temperature [rel.un.] (=T/5.93e9 K)
*        DENR - electron number density n_e [rel.un.]
* Output: DENR1 - the density adjusted to the found CHI
*         CHI=\mu/T, where \mu is the chem.pot-l w/o the rest-energy
*               (= n_e/1.7366e31 cm^{-3})
*         FEid - free energy / N_e kT,
*         PEid - pressure (P_e) / n_e kT
*         UEid - internal energy / N_e kT
*         SEid - entropy, CVE - heat capacity [per 1 electron, in k_B],
*         CHITE=(d ln P_e / d ln T)_V, CHIRE=(d ln P_e / d ln n_e)_T
*         DDENR = (d n / d\mu)_T (rel.un.)
*  All quantities are by default in relativistic units
*  B - magnetic field (=h\omega_B/mc^2)
      parameter (BOHR=137.036,PI=3.141592653d0) ! universal const.
      parameter (PI2=PI**2,BOHR2=BOHR**2,BOHR3=BOHR2*BOHR)
      parameter (NLMAX1=50,PRECLN=9.,EPS=1.D-4,NI=50) ! accuracy param.
      BT=B/TEMR
      if (BT.lt.10.) then ! reduce max.number of Landau levels
         NLMAX=NLMAX1
      else
         NLMAX=NLMAX1/dsqrt(BT/10.)
         NLMAX=max0(10,NLMAX)
      endif
      CNLMAX=NLMAX
      PF0=(29.6088132d0*DENR)**.33333333d0 ! Classical Fermi momentum
      if (B.gt.0.) then
         PFB=PF0**3/(1.5*B)
         PF=dmin1(PF0,PFB) ! magnetic Fermi momentum
        if (PF.gt.1.d-7) then
           TF=dsqrt(1.d0+PF**2)-1.d0 ! Nonmagn.Fermi T
        else
           TF=0.5d0*PF**2
        endif
         DE=TF+PRECLN*TEMR ! (1+DE)=Upper E bound to int.
         CNL0=DE*(1.d0+0.5d0*DE)/B ! Max contributing Landau number
      else
         CNL0=NLMAX+1
      endif
* Preliminary nonmagnetic calculation
      call CHEMFIT7(DENR,TEMR,CHI,CMU1,1,
     *  CMUDENR,CMUDT,CMUDTT)
      DDENR=1.d0/CMUDENR
      TEMP=TEMR*BOHR2
      call ELECT11(TEMP,CHI,
     *  DENSNM,FEidNM,PEidNM,UEidNM,SEidNM,CVENM,CHITENM,CHIRENM,DDNNM,
     *  DlnDT,DlnDHH,DlnDTT,DlnDHT)
      DENR1NM=DENSNM/BOHR3
      if (CNL0.gt.CNLMAX) then ! Nonmagnetic case
         DENR1=DENR1NM
         FEid=FEidNM
         PEid=PEidNM
         UEid=UEidNM
         SEid=SEidNM
         CVE=CVENM
         CHITE=CHITENM
         CHIRE=CHIRENM
         DDN=DDNNM
        return
      else ! Magnetic calculation:
         ZET=(dmin1(1.d0,PF0**2/B/1.5)*PF0)**2
         NL0=CNL0+.5
      endif
*   ----------  Calculation of chemical potential CMU   ------------   *'
      Kup=0.
      Klow=0.
      POWER=.666666667d0
      if (NL0.lt.1) POWER=2.d0
      DENR1=0.
      do 5 I=1,NI
      if (ZET.gt.3.d0*TEMR) then ! Degenerate plasma
         SHIFTMU=0.
      else ! Non-degenerate plasma
         SHIFTMU=1.5d0*TEMR*dlog(3.d0*TEMR/ZET)
         if (NL0.lt.2) SHIFTMU=SHIFTMU/3.d0
      endif
      if (SHIFTMU.lt..5d0) then
         CMU=sqrt(1.+ZET)*(1.-SHIFTMU)
      else
         CMU=sqrt(1.+ZET)-SHIFTMU
      endif
      CHI=(CMU-1.d0)/TEMR
      DENR1old=DENR1
      call ELECTMAG(B,TEMR,CHI,
     *  DENR1,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,DDENR)
      if (DENR1.gt.DENR) then
        if (Kup.eq.0) Zup=ZET
         Kup=1
         Zup=dmin1(Zup,ZET)
      else
        if (Klow.eq.0) Zlow=ZET
         Klow=1
         Zlow=dmax1(Zlow,ZET)
      endif
      ZET1=ZET
      DENRMIN=dmin1(1.d-19,DENR/1.d2)
      DENR1=dmax1(DENRMIN,DENR1)
      ZET=ZET1*(DENR/DENR1)**POWER
      if (Kup.ne.0.and.(ZET.gt.Zup.or.(I.gt.NI/2.and.DENR1.lt.DENR)))
     &   ZET=(ZET1+Zup)/2.d0
      if (Klow.ne.0.and.(ZET.lt.Zlow.or.(I.gt.NI/2.and.DENR1.gt.DENR)))
     &   ZET=(ZET1+Zlow)/2.d0
      if (abs(DENR/DENR1-1.).lt.EPS) goto 6
      if (dabs(DENR1-DENR1old).lt.EPS*.1d0*DENR1.
     .   and.DENR1.gt.DENRMIN*(1.d0+EPS)) goto 55
    5 continue
   55 continue
      if (dabs(DENR1-DENR).gt.1.d1*DENR*EPS) ! improvement of 02.09.11
     *  print*,'Warning CHEMAG8: INCORRECT density:',DENR1,'  i.o.',DENR
    6 continue
* Interpolation in order to exclude jump at nonmagnetic/magnetic switch
      Y=(1.d0-CNL0/CNLMAX)**3
      X=1.d0-Y
      DENR1=X*DENR1NM+Y*DENR1
      FEid=X*FEidNM+Y*FEid
      PEid=X*PEidNM+Y*PEid
      UEid=X*UEidNM+Y*UEid
      SEid=X*SEidNM+Y*SEid
      CVE=X*CVENM+Y*CVE
      CHITE=X*CHITENM+Y*CHITE
      CHIRE=X*CHIRENM+Y*CHIRE
      return
      end

* ==============  IDEAL ELECTRONS  =================================== *
      subroutine ELECTMAG(B,TEMR,CHI,
     *  DENR,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,DDENR)
*                                                       Version 04.03.09
* BLIN8 --> BLIN9                                               24.11.10
* Ideal magnetized electron-gas EOS.
* Input: B - magnetic field [rel.un.] (=B/4.414e13 G)
*        TEMR - temprature [rel.un.] (=T/5.93e9 K)
*        CHI=\mu/T, where \mu is the chem.pot-l w/o the rest-energy
* Output: DENR - electron number density n_e [rel.un.]
*               (= n_e/1.7366e31 cm^{-3})
*         FEid - free energy / N_e kT,
*         PEid - pressure (P_e) / n_e kT
*         UEid - internal energy / N_e kT
*         SEid - entropy, CVE - heat capacity [per 1 electron, in k_B],
*         CHITE=(d ln P_e / d ln T)_V, CHIRE=(d ln P_e / d ln n_e)_T
*         DDENR=(d n_e / d \mu)_T
      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)
      parameter (NLMAX=200,PRECLN=9.,EPST=1.d-17,EPS14=3.d0,EPSR=.1d0)
      if (B.gt.0.) then
         DE=(dmax1(0.d0,CHI)+PRECLN)*TEMR ! (1+DE)=Upper E bound to int.
         NL0=DE*(1.+DE/2.)/B ! Max contributing Landau number to integr.
      else
         NL0=NLMAX+1
      endif
      if (NL0.gt.NLMAX) then ! nonmagnetic calc.
         TEMP=TEMR*BOHR2
         call ELECT11(TEMP,CHI,
     *     DENS,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE,DDN,
     *     DlnDT,DlnDHH,DlnDTT,DlnDHT)
         DENR=DENS/BOHR3
         DDENR=DENR*DDN/TEMR
        return
      endif
      if (CHI.lt.-40.d0) then ! nondegenerate
         CLE=dsqrt(2.d0*PI/TEMR) ! thermal length
         DENR=2.d0*dexp(CHI)/CLE**3
        if (TEMR.gt.EPSR) ! this case is not suitable for this code
     *    stop'ELECTMAG: nondegenerate, magnetic, relativistic case'
         call MAGNION(B,TEMR,DENR,1.d0,5.4858d-4,2.d0,2,
     *     FIid,PIid,UIid,SIid,CVI,CHITI,CHIRI)
         DDENR=DENR/TEMR
        return
      endif
* Magnetic calculation:
      SP=0.
      SN=0.
      SNDH=0.
      SNDT=0.
      SPDH=0.
      SPDT=0.
      SNDHH=0.
      SNDTT=0.
      SNDHT=0.
      SPDHH=0.
      SPDTT=0.
      SPDHT=0.
      do NL=0,NL0
         EN1=B*NL
         EN0=dsqrt(1.d0+2.d0*EN1)
        if (EN1.gt.1.d-7) EN1=EN0-1.d0
         CHIN=CHI-EN1/TEMR
         TN=TEMR/EN0
         SQREN=dsqrt(EN0)
         call BLIN9(TN,CHIN,
     *      W0,W0DX,W0DT,W0DXX,W0DTT,W0DXT,
     *      W1,W1DX,W1DT,W1DXX,W1DTT,W1DXT,
     *      W2,W2DX,W2DT,W2DXX,W2DTT,W2DXT,
     *      W0XXX,W0XTT,W0XXT)
         H14=EPS14
        if (CHIN.lt.14.d0.and.14.d0-CHIN.lt.H14) then ! regularization
           CHINM=14.d0-H14
           CHINP=14.d0
           XM=(CHIN-CHINM)/H14
           XP=1.d0-XM
           call BLIN9(TN,CHINM,
     *        W0M,W0DXM,W0DTM,W0DXXM,W0DTTM,W0DXTM,
     *        W1M,W1DXM,W1DTM,W1DXXM,W1DTTM,W1DXTM,
     *        W2M,W2DXM,W2DTM,W2DXXM,W2DTTM,W2DXTM,
     *        W0XXXM,W0XTTM,W0XXTM)
           call BLIN9(TN,CHINP,
     *        W0P,W0DXP,W0DTP,W0DXXP,W0DTTP,W0DXTP,
     *        W1P,W1DXP,W1DTP,W1DXXP,W1DTTP,W1DXTP,
     *        W2P,W2DXP,W2DTP,W2DXXP,W2DTTP,W2DXTP,
     *        W0XXXP,W0XTTP,W0XXTP)
           W0=XP*W0+XM*W0P
           W0DX=XP*W0DX+XM*W0DXP
           W0DT=XP*W0DT+XM*W0DTP
           W0DXX=XP*W0DXX+XM*W0DXXP
           W0DXT=XP*W0DXT+XM*W0DXTP
           W0DTT=XP*W0DTT+XM*W0DTTP
           W0XXX=XP*W0XXX+XM*W0XXXP
           W0XTT=XP*W0XTT+XM*W0XTTP
           W0XXT=XP*W0XXT+XM*W0XXTP
           W1=XP*W1+XM*W1P
           W1DX=XP*W1DX+XM*W1DXP
           W1DT=XP*W1DT+XM*W1DTP
           W1DXX=XP*W1DXX+XM*W1DXXP
           W1DXT=XP*W1DXT+XM*W1DXTP
           W1DTT=XP*W1DTT+XM*W1DTTP
           W2=XP*W2+XM*W2P
           W2DX=XP*W2DX+XM*W2DXP
           W2DT=XP*W2DT+XM*W2DTP
           W2DXX=XP*W2DXX+XM*W2DXXP
           W2DXT=XP*W2DXT+XM*W2DXTP
           W2DTT=XP*W2DTT+XM*W2DTTP
        endif
         SN=SN+W0DX*SQREN
         SP=SP+SQREN*W0
         SNDH=SNDH+W0DXX*SQREN
         SPDH=SPDH+SQREN*W0DX
         SNDHH=SNDHH+W0XXX*SQREN
         SPDHH=SPDHH+SQREN*W0DXX
         CHINDT=EN1/TEMR**2
         CHINDTT=-2.d0*EN1/TEMR**3
         SNDT=SNDT+(W0DX/(2.d0*TEMR)+W0DXX*CHINDT+W0DXT/EN0)*SQREN
         SPDT=SPDT+SQREN*(1.5d0*W0/TEMR+CHINDT*W0DX+W0DT/EN0)
         SNDTT=SNDTT+(-.25d0*W0DX/TEMR**2-CHINDT/TEMR*W0DXX+
     +     W0DXT/(TEMR*EN0)+CHINDT**2*W0XXX+2.d0*CHINDT/EN0*W0XXT+
     +     W0XTT/EN0**2)*SQREN
         SPDTT=SPDTT+SQREN*(.75d0/TEMR**2*W0+CHINDT/TEMR*W0DX+
     +     3.d0*W0DT/(EN0*TEMR)+CHINDT**2*W0DXX+2.d0*CHINDT/EN0*W0DXT+
     +     W0DTT/EN0**2)
         SNDHT=SNDHT+(.5d0*W0DXX/TEMR+CHINDT*W0XXX+W0XXT/EN0)*SQREN
         SPDHT=SPDHT+SQREN*(1.5d0*W0DX/TEMR+CHINDT*W0DXX+W0DXT/EN0)
         SCHT1=SPDT-SPDH*SNDT/SNDH
        if (NL.eq.0) then
           SN=SN/2.d0
           SP=SP/2.d0
           SNDH=SNDH/2.d0
           SPDH=SPDH/2.d0
           SNDT=SNDT/2.d0
           SPDT=SPDT/2.d0
           SNDHH=SNDHH/2.d0
           SPDHH=SPDHH/2.d0
           SNDTT=SNDTT/2.d0
           SPDTT=SPDTT/2.d0
           SNDHT=SNDHT/2.d0
           SPDHT=SPDHT/2.d0
           SCHT=SCHT1/2.d0
        else
           if(SCHT1.ge.SCHT.or.NL.lt.NL0-1) SCHT=SCHT1
        endif
      enddo ! next NL
    1 continue
      CN=dsqrt(2.d0*TEMR)/PI2*B
      DENR=CN*SN ! n_e [rel.un.]
      PRI=DENR*TEMR ! ideal-gas pressure = nkT = normalization factor
      CNP=CN*TEMR
      PR=CNP*SP ! P [rel.un.]
      PEid=PR/PRI ! P/nkT
      FR=CHI*PRI-PR ! F/V [rel.un.]
      FEid=FR/PRI ! F/NkT
* derivatives over chi at constant T:
      dndH=CN*SNDH ! (d n_e/d\chi)_T
      dPdH=CNP*SPDH ! (d P_e/d\chi)_T
      dndHH=CN*SNDHH ! (d^2 n_e/d\chi^2)_T
      dPdHH=CNP*SPDHH ! (d^2 P_e/d\chi^2)_T
      dFdH=(DENR+CHI*dndH)*TEMR-dPdH ! [d (F/V)/d\chi]_T
      dFdHH=(2.d0*dndH+CHI*dndHH)*TEMR-dPdHH ! [d^2 (F/V)/d\chi^2]_T
* derivatives over T at constant chi:
      dndT=CN*SNDT ! (d n_e/dT)_\chi
      dPdT=CNP*SPDT ! (dP/dT)_\chi
      dFdT=CHI*(DENR+TEMR*dndT)-dPdT
      dndTT=CN*SNDTT
      dPdTT=CNP*SPDTT
      dFdTT=chi*(2.d0*dndT+TEMR*dndTT)-dPdTT
* mixed derivatives
      dndHT=CN*SNDHT
      dPdHT=CNP*SPDHT
      dFdHT=DENR+CHI*dndH+TEMR*dndT+CHI*TEMR*dndHT-dPdHT
* normalized thermodynamic functions:
      G=dndT/dndH ! = - (d\chi/dT)_V
      GDT=(dndTT-dndT*dndHT/dndH)/dndH
      GDH=(dndHT-dndT*dndHH/dndH)/dndH
      S=dFdH*G-dFdT ! S=entropy per volume=-(dF/dT)_V
      SEid=S/DENR ! S/Nk
      UEid=FEid+SEid
      SDT=dFdHT*G+dFdH*GDT-dFdTT
      SDH=dFdHH*G+dFdH*GDH-dFdHT
      CV=TEMR*(SDT-SDH*G) ! C_V per volume = T(dS/dT)_V
      CVE=CV/DENR ! C_V/Nk
C      CHITE1=TEMR/PR*(dPdT-dPdH*G) ! d ln P / d ln T |_V
      CHITE=TEMR/PR*CNP*SCHT ! d ln P / d ln T |_V
      if (dabs(dPdT-dPdH*G).lt.EPST*dPdT) then ! accuracy loss; non-mag.
         TEMP=TEMR*BOHR2
         call ELECT11(TEMP,CHI,
     *     DENS,FEnm,PEnm,UEnm,SEnm,CVEnm,CHITE,CHIREnm,DDNnm,
     *     DlnDT,DlnDHH,DlnDTT,DlnDHT)
         print*,'ELECTMAG: Warning: CHITE switched to non-magnetic'
      endif
      CHIRE=DENR/PR*dPdH/dndH ! d ln P / d ln rho |_T
      DDENR=dndH/TEMR ! (d n / d\mu)_T (rel.un.)
      return
      end

* =============================   IONS   ============================= *
      subroutine MAGNION(B,TEMR,DENRI,Zion,CMI,GYRO,MULTI,
     *   FIid,PIid,UIid,SIid,CVI,CHITI,CHIRI)
*                                                       Version 02.10.12
* Previous version was 04.03.09. The difference: SPINION is separate now
* Thermodynamic functions of an ideal magnetized gas of atomic nuclei
*   (see the "EOS" written notes, page 129)
* Input: B - magnetic field [rel.un.] (=B/4.414e13 G)
*        TEMR - temperature [rel.un.] (=T/5.93e9 K)
*        DENRI - number density of ions [rel.un.]
*        Zion,CMI - atomic charge and mass numbers of the ions
*        GYRO - g-factor
*        MULTI - spin multiplicity (2S+1)
* Output: FIid - free energy per ion per kT,
*         PIid=1 - pressure /(n_i kT)
*         UIid - internal energy per ion per kT,
*         SIid - entropy /k, per ion,
*         CVI - specific heat /k, per ion,
*         CHITI=CHIRI=1 - log.derivatives of pressure.
* NB: zero-point term is not included.
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter (AUM=1822.888d0,PI=3.141592653d0,EPS=1.d-8)
      parameter (TwoPI=2.d0*PI)
      RMI=CMI*AUM ! nucleus-to-electron mass ratio
      CLI=dsqrt(TwoPI/(RMI*TEMR)) ! thermal length
      ZETI=Zion/RMI*B/TEMR ! \hbar\omega_{ci}/kT
      if (ZETI.lt.EPS) then ! non-magnetic case
         FIid=dlog(DENRI*CLI**3)-1.d0
         UIid=1.5d0         
         CVI=1.5d0
      else ! magnetic case
         EZ1=dexp(-ZETI)
         EZ=1.d0-EZ1
         FIid=dlog(TwoPI/B*CLI*DENRI*EZ/Zion)-1.d0+.5d0*ZETI
         ZETI1=ZETI*EZ1/EZ
         UIid=.5d0+ZETI1+.5d0*ZETI
         CVI=.5d0+ZETI1**2
      endif
      if (UIid.lt.FIid) stop'MAGNION: SIid < 0 (degenerate ions)'
      PIid=1.d0
      CHITI=1.d0
      CHIRI=1.d0
      call SPINION(GYRO,ZETI,MULTI,Fspin,Uspin,CVspin)
      FIid=FIid+Fspin
      UIid=UIid+Uspin
      CVI=CVI+CVspin
      SIid=UIid-FIid
      return
      end

      subroutine EOSFIM12(LIQSOL,CMI,Zion,RS,GAMI,TEMP,GAMAG,KWK,KDIR,
     *  FC1,UC1,PC1,SC1,CV1,PDT1,PDR1,
     *  FC2,UC2,PC2,SC2,CV2,PDT2,PDR2)
*                                                       Version 03.06.12
* Stems from EOSFIM9 v.17.01.11
* New (as of 2012): the magnetic effects on phonon gas are included,
*          spin multiplicity in solid is made consistent with the liquid
* Non-ideal parts of thermodynamic functions in the fully ionized plasma
* Input: LIQSOL=0/1(liquid/solid), 
*        Zion,CMI - ion charge and mass numbers,
*        RS=r_s (electronic density parameter),
*        GAMI=Gamma_i (ion coupling),
*        TEMP - temperature in a.u.
*        GAMAG=B/2.35e9 G (magn.field in a.u.),
*        KWK=1 switches quantum diffraction for ion liquid on,
*        KDIR=1 switches an alternative direction of the crystal on.
* 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 WITHOUT magnetic field. 
* This is useful for preventing accuracy loss in some cases
*   (e.g., when SC2 << SC1).
* IDEAL-GAS MAGNETIC-FIELD CONTRIBUTION FOR THE IONS SHOULD BE ADDED
*   EXTERNALLY.
* 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.
      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'EOSFIM12: invalid LIQSOL'
      if (CMI.le..1) stop'EOSFIM12: too small CMI'
      if (Zion.le..1) stop'EOSFIM12: too small Zion'
      if (RS.le..0) stop'EOSFIM12: invalid RS'
      if (GAMI.le..0) stop'EOSFIM12: invalid GAMI'
      GAME=GAMI/Zion**1.666667
      call EXCORM(RS,GAME,GAMAG,
     *  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).
      FWK=0.
      UWK=0.
      PWK=0.
      CVWK=0.
      PDRWK=0.
      PDTWK=0.
      if (LIQSOL.eq.0) then                 ! liquid
         call FITION9(GAMI,FION,UION,PION,CVii,PDTii,PDRii)
         FItot=FION+FidION
         UItot=UION+1.5d0
         PItot=PION+1.d0
         CVItot=CVii+1.5d0
         SCItot=UItot-FItot
         PDTi=PDTii+1.d0
         PDRi=PDRii+1.d0
* Wigner-Kirkwood (quantum diffraction) terms:
        if (KWK.eq.1) then
           ZETI=Zion/(CMI*AUM)*GAMAG/TEMP ! \hbar\omega_{ci}/kT
           call QDMAG(TPT,ZETI,FWK,UWK,PWK,CVWK,PDTWK,PDRWK)
        endif
      else ! solid
         BP=Zion/(CMI*AUM)*GAMAG/TEMP/TPT ! \hbar\omega_{ci}/kT_p
         call FHARMAG(GAMI,TPT,BP,
     *     Fharm,Uharm,Pharm,CVharm,Sharm,PDTharm,PDRharm,EVAR)
        if (KDIR.eq.1) then ! less favorable orientation of the crystal
           Fharm=Fharm+EVAR
           Uharm=Uharm+EVAR
        endif
         call ANHARM8(GAMI,TPT,Fah,Uah,Pah,CVah,PDTah,PDRah) ! anharm.
         FItot=Fharm+Fah ! spin term -.6931=-ln2 is excluded 28.04.12
         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 ! spin term .6931=ln2 is excluded 28.04.12
         CVii=CVItot-1.5d0 ! minus 1.5=ideal-gas
      endif
* Calculate "ie" part: NB: here the ie-scaling is abandoned!
      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 EXCHANGE AND CORRELATION   ================ *
      subroutine EXCORM(RS,GAME,GAMAG,
     *  FXC,UXC,PXC,CVXC,SXC,PDTXC,PDRXC) ! "ee"="xc"
*                                                       Version 06.05.08
*                          "d0" is added at some exact numerals 27.04.12
* Exchange-correlation contribution for the electron gas
* Stems from EXCOR7 v.09.06.07 and TANAKMAG v.27.04.98
* Input: RS - electron density parameter =electron-sphere radius in a.u.
*        GAME - electron Coulomb coupling parameter
*        GAMAG - atomic magnetic-field 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
      parameter(TINY=1.d-99)
      if (RS.le.0.) stop'wrong RS'
      if (GAME.le.0.) stop'wrong GAME'
      THETA1=.543*RS/GAME ! non-relativistic non-mag.degeneracy par.
*  ------   06.04.98:  magnetic scale factor on RS at high T   ------  *
      U=GAMAG/2.d0*GAME*RS ! hbar omega_c/(2T)...
      SCALE=1.d0
      SCADU=0.
      SCADUU=0.
      if (U.gt.1.d-4) then
         THU=dtanh(U)
         TU=THU/U
         S=dsqrt(1.d0-TU)
         CH2=2.d0
         TDU=-TU/U
         CDU=0.
         CDUU=0.
        if (U.lt.20.) then
           CHU2=dcosh(U)**2
           CH2=dcosh(2.d0*U)/CHU2
           VU=darcth(S)/S
           TDU=1.d0/(U*CHU2)+TDU
           TDUU=2.d0/U**2*TU-2./(U**2*CHU2)-2.d0*TU/CHU2
           CDU=2.d0*THU/CHU2
           CDUU=(6.d0/CHU2-4.d0)/CHU2
        else
           TDUU=2.d0/U**2*TU
           VU=dlog(4.d0*U-1.d0)/(2.d0-1.d0/U)
        endif
         VDU=.5d0*TDU/S**2*(VU-1.d0/TU)
         VDUU=VDU*TDUU/TDU+(VDU*TDU*1.5d0+(TDU/TU)**2*.5d0)/S**2
         SCALE=CH2*TU*VU ! =f_2 of Steinberg=f_1 of PCS99
         SCADU=CDU*TU*VU+CH2*TDU*VU+CH2*TU*VDU
         SCADUU=CDUU*TU*VU+CH2*TDUU*VU+CH2*TU*VDUU+
     +     2.d0*(CDU*TDU*VU+CDU*TU*VDU+CH2*TDU*VDU)
      endif
      UDG=GAMAG*THETA1*GAME/.543
      SCADH=SCADU*GAMAG/1.086*GAME**2
      SCADG=SCADU*UDG
      SCADHH=SCADUU*(GAMAG*GAME**2/1.086)**2
      SCADGG=SCADUU*UDG**2+SCADU*GAMAG*THETA1/.543
      SCADHG=GAMAG*GAME/.543*(SCADU+SCADUU*GAME/2.d0*UDG)
*  ------ Finished calculation of SCALE and its derivatives.
      GR4=(GAMAG*RS**2)**2+TINY
*   ----------------------   scaled theta and its derivatives:
      THETAM=.090063*GR4*RS/GAME ! magn.theta; .09=8/(9\pi^2)
      TMT=.16586*GR4 ! theta_m/theta_0
      TMTDH=4.d0*TMT/THETA1
      TMTDG=4.d0*TMT/GAME
      TMTDHH=12.d0*TMT/THETA1**2
      TMTDGG=12.d0*TMT/GAME**2
      TMTDHG=16.d0*TMT/(THETA1*GAME)
      SCUP=1.d0+TMT
      E=dexp(-1.d0/THETAM)
      EDH=5.d0*E/(THETA1*THETAM)
      EDG=4.d0*E/(THETAM*GAME)
      EDHH=EDH/THETA1*(5.d0/THETAM-6.d0)
      EDGG=EDG/GAME*(4.d0/THETAM-5.d0)
      EDHG=EDH/GAME*(4.d0/THETAM-4.d0)
      PM=SCALE*TMT*E
      PMDH=TMTDH*SCALE*E+TMT*E*SCADH+TMT*EDH*SCALE
      PMDG=TMTDG*SCALE*E+TMT*E*SCADG+TMT*EDG*SCALE
      PMDHH=TMTDHH*SCALE*E+TMT*E*SCADHH+TMT*EDHH*SCALE+
     +  2.d0*(TMTDH*SCADH*E+TMTDH*SCALE*EDH+TMT*SCADH*EDH)
      PMDGG=TMTDGG*SCALE*E+TMT*E*SCADGG+TMT*EDGG*SCALE+
     +  2.d0*(TMTDG*SCADG*E+TMTDG*SCALE*EDG+TMT*SCADG*EDG)
      PMDHG=TMTDHG*SCALE*E+TMT*E*SCADHG+TMT*EDHG*SCALE+
     +  E*(TMTDH*SCADG+TMTDG*SCADH)+SCALE*(TMTDH*EDG+TMTDG*EDH)+
     +  TMT*(SCADH*EDG+SCADG*EDH)
      SCDN=1.d0+PM
      THETA=THETA1*SCUP/SCDN ! SCALED degeneracy parameter
      THDH=(SCUP+THETA1*TMTDH-THETA*PMDH)/SCDN
      THDG=(THETA1*TMTDG-THETA*PMDG)/SCDN
      THDHH=(2.d0*(TMTDH+(THETA*PMDH-SCUP)/SCDN*PMDH)+
     +  THETA1*(TMTDHH-2.d0*TMTDH*PMDH/SCDN)-THETA*PMDHH)/SCDN
      THDGG=(THETA1*TMTDGG+
     +  2.d0*(THETA*PMDG-THETA1*TMTDG)*PMDG/SCDN-THETA*PMDGG)/SCDN
      THDHG=(TMTDG+THETA1*TMTDHG-THETA*PMDHG+
     +  ((2.d0*THETA*PMDH-SCUP)*PMDG-
     -  THETA1*(TMTDG*PMDH+TMTDH*PMDG))/SCDN)/SCDN
*   ----------------------------------------------------------------   *
      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./THETA)
         T2=SHT2/CHT2 ! dtanh(1./dsqrt(THETA))
         T1DH=-1.d0/(THETA*CHT1)**2 ! d T1 / d\theta
         T1DHH=2.d0/(THETA*CHT1)**3*(CHT1-SHT1/THETA)
         T2DH=-.5d0*SQTH/(THETA*CHT2)**2
         T2DHH=(.75d0*SQTH*CHT2-.5d0*SHT2)/(THETA*CHT2)**3
      else
         T1=1.d0
         T2=1.d0
         T1DH=0.
         T2DH=0.
         T1DHH=0.
         T2DHH=0.
      endif
      A0=.75d0+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.d0+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.d0+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=.5d0/THETA+T2DH/T2+B0DH/B0-B1DH/B1
      BDH=B*BH
      BDHH=BDH*BH+B*(-.5d0/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.d0+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=.5d0/THETA+T2DH/T2+D0DH/D0-D1DH/D1
      DDH=D*DH
      DDHH=DDH*DH+D*(-.5d0/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.d0+2.555501*THETA2+.146319*THETA4
      E1DH=5.111002*THETA+.585276*THETA3
      E1DHH=5.111002+1.755828*THETA2
      E=THETA*T1*E0/E1
      EH=1.d0/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.d0/THETA2)
      EXP1TH=dexp(-1.d0/THETA)
      C=(.872496+.025248*EXP1TH)*E
      CDH=.025248*EXP1TH/THETA2*E+C*EDH/E
      CDHH=.025248*EXP1TH/THETA2*(EDH+(1.d0-2.d0*THETA)/THETA2*E)+
     +  CDH*EDH/E+C*EDHH/E-C*(EDH/E)**2
      DISCR=dsqrt(4.d0*E-D**2)
      DIDH=.5d0/DISCR*(4.d0*EDH-2.d0*D*DDH)
      DIDHH=(-((2.d0*EDH-D*DDH)/DISCR)**2+2.d0*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.d0*CDH*DDH+C*DDHH)/E+
     +  (2.d0*(CDH*D+C*DDH-C*D*EDH/E)*EDH+C*D*EDHH)/E**2
      SQGE=dsqrt(GAME)
      S2=-2.d0/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=.5d0*S2/GAME
      S2DGG=-.5d0*S2DG/GAME
      S2DHG=.5d0*S2DH/GAME
      R3=E*GAME+D*SQGE+1.d0
      R3DH=EDH*GAME+DDH*SQGE
      R3DHH=EDHH*GAME+DDHH*SQGE
      R3DG=E+.5d0*D/SQGE
      R3DGG=-.25d0*D/(GAME*SQGE)
      R3DHG=EDH+.5d0*DDH/SQGE
      B3=A-C/E
      B3DH=ADH-CDH/E+C*EDH/E**2
      B3DHH=ADHH-CDHH/E+(2.d0*CDH*EDH+C*EDHH)/E**2-2.d0*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.d0*D*B2*EDH/E**3-B3DH/E
      C3DHH=(-B3DHH+
     +  (DDHH*B2+2.d0*DDH*B2DH+D*B2DHH+B3DH*EDH+B3*EDHH+B3DH*EDH)/E-
     -  2.d0*((DDH*B2+D*B2DH+B3*EDH+DDH*B2+D*B2DH)*EDH+D*B2*EDHH)/E**2+
     +  6.d0*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.d0-D**2/E
      B4DH=EDH*(D/E)**2-2.d0*D*DDH/E
      B4DHH=EDHH*(D/E)**2+2.d0*EDH*(D/E)**2*(DDH/D-EDH/E)-
     -  2.d0*(DDH**2+D*DDHH)/E+2.d0*D*DDH*EDH/E**2
      C4=2.d0*E*SQGE+D
      C4DH=2.d0*EDH*SQGE+DDH
      C4DHH=2.d0*EDHH*SQGE+DDHH
      C4DG=E/SQGE
      C4DGG=-.5d0*E/(GAME*SQGE)
      C4DHG=EDH/SQGE
      S4A=2.d0/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.d0*DDH*B3DH+D*B3DHH+
     +  B4DHH*B2+2.d0*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.d0*(DISCR*DIDH+C4*C4DH)/DN1**2-
     -  (DDHH*DISCR-D*DIDHH)/DN2+UP2*2.d0*(DISCR*DIDH+D*DDH)/DN2**2
      S4CDG=C4DG*DISCR/DN1
      S4CDGG=C4DGG*DISCR/DN1-2.d0*C4*DISCR*(C4DG/DN1)**2
      S4CDHG=(C4DHG*DISCR+C4DG*DIDH-
     -  C4DG*DISCR/DN1*2.d0*(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.d0*(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
      FXCDH1=S1DH+S2DH+S3DH+S4DH
      FXCDG1=S1DG+S2DG+S3DG+S4DG
      FXCDHH1=S1DHH+S2DHH+S3DHH+S4DHH
      FXCDGG1=S2DGG+S3DGG+S4DGG
      FXCDHG1=S1DHG+S2DHG+S3DHG+S4DHG
* -------- Magnetic correction of derivatives
       FXCDH=FXCDH1*THDH
       FXCDG=FXCDG1+FXCDH1*THDG
       FXCDHH=FXCDHH1*THDH**2+FXCDH1*THDHH
       FXCDGG=FXCDGG1+FXCDH1*THDGG+(2.d0*FXCDHG1+FXCDHH1*THDG)*THDG
       FXCDHG=FXCDHG1*THDH+FXCDH1*THDHG+FXCDHH1*THDH*THDG
* --------
      PXC=(GAME*FXCDG-2.d0*THETA1*FXCDH)/3.d0
      UXC=GAME*FXCDG-THETA1*FXCDH
      SXC=(GAME*S2DG-S2+GAME*S3DG-S3+S4A*S4B*(GAME*S4CDG-S4C))-
     -  THETA1*FXCDH
      if (dabs(SXC).lt.1.d-9*dabs(THETA1*FXCDH)) SXC=0. ! accuracy loss
      CVXC=2.d0*THETA1*(GAME*FXCDHG-FXCDH)-
     -  THETA1**2*FXCDHH-GAME**2*FXCDGG
      if (dabs(CVXC).lt.1.d-9*dabs(GAME**2*FXCDGG)) CVXC=0. ! accuracy
      PDLH=THETA1*(GAME*FXCDHG-2.d0*FXCDH-2.d0*THETA1*FXCDHH)/3.d0
      PDLG=GAME*(FXCDG+GAME*FXCDGG-2.d0*THETA1*FXCDHG)/3.d0
      PDRXC=PXC+(PDLG-2.d0*PDLH)/3.d0
      PDTXC=GAME*(THETA1*FXCDHG-GAME*FXCDGG/3.d0)-
     -  THETA1*(FXCDH/.75d0+THETA1*FXCDHH/1.5d0)
      return
      end

      function darcth(X)
*                                                       Version 26.03.98
      implicit double precision (A-H), double precision (O-Z)
      if (dabs(X).ge.1.d0) stop'darcth: X out of range.'
      if (dabs(X).lt.1.d-9) then
         D=X
      else
         D=.5d0*dlog((1.d0+X)/(1.d0-X))
      endif
      darcth=D
      return
      end

* ==============   SUBROUTINES FOR THE SOLID STATE   ================= *
      subroutine HLMAG(TPT,BP,UMth,CMth,SMth,PMth,PDTth,PDRth,
     *  U1m,dLNu1,d2LNu1,DU1,UM,SM,CMU,CMS)
*                                                       Version 09.07.12
* Input: TPT=T_p/T
*        BP=\omega_c/T_p
* Output: UMth - thermal energy / NkT
*         CMth - specific heat/Nk
*         SMth - entropy/Nk
* Auxiliary output:
*         U1m=u1=<\omega_{phonon}/\omega_p> - normalized 1st moment
*             NB: zero-point cyclotron oscillations are excluded,
*             i.e., \omega_{phonon}=\omega_{phonon,total}-\omega_c/3.
*         dLNu1=d\ln(u1)/d\ln(BP)
*         d2LNu1=d^2\ln(u1)/[d\ln(BP)]^2
*         DU1 - variation of u1 with changing direction of the field
* Additional provisional output:
*         UM - a simpler fit to thermal energy / NkT
*         SM - a simpler fit to entropy / Nk
*         CMU - specific heat derived from UM
*         CMS - specific heat derived from SM, = d(SM)/dt,
*              where t=ln(TTP)=-ln(TPT)
* Useful internal variables:
*         dUMdt=d(UM)/dt, dUMdb=d(UM)/db, where b=ln(BP),
*         dUMdtt=d^2(UM)/dt^2, dUMdbb=d^2(UM)/db^2, dUMdtb=d^2(UM)/dtdb,
*         dSMdb=d(SM)/db, dCMSd=d(CMS)/dt=d^2(SM)/dt^2,
*         dCMSdt=d^2(SM)/dt^2, dSMdbb=d^2(SM)/db^2, dSMdtb=d^2(SM)/dtdb
      implicit double precision (A-H), double precision (O-Z)
      parameter(TINY=1.d-19)
      parameter(Zb=12.5d0,A1=24.d0,A3=15.d0,A4=119.d0,A6=0.4d0)
      parameter(A4S=9.8d0,A3S=7.9d0,A1S=42.d0)
      parameter(U1inf=0.1801,CU1=1.27)
      if (TPT.lt.0.) stop'HLMAG: TPT < 0'
      TTP=1.d0/(TPT+TINY)
      call HLfit12(TPT,F,U,CVth,Sth,U1,CWth,1)
      if (BP.lt.0.) stop'HLMAG: BP < 0'
      if (BP.eq.0.) then
         UMth=U
         CMth=CVth
         SMth=Sth
         UM=U
         CMU=CVth
         CMS=CVth
         SM=Sth
         PMth=U/2.d0
         PDTth=CVth/2.d0
         PDRth=.75d0*U-.25d0*CVth
         U1m=U1
         dLNu1=0.
         d2LNu1=0.
         DU1=0.
      else
         TA=TTP/BP
         TB=TTP*BP
         SqrtTB=sqrt(TB)
* U:
         TA3=3.d0*TA
         TQ3=(TA3)**2*dsqrt(TA3)
         Q=.5d0/(1.d0+TQ3)**A6
         UMH=U/(1.d0+Q)
         TV=(Zb*SqrtTB+A4*TB)*TB
         DV=1.d0+TV+A3*TB
         TW=A1*TTP**2
         DW=1.d0+TW
         UML=TV/DV/dsqrt(DW)
         UM=UML+UMH
* S:
         TVS=Zb*(SqrtTB/.6d0+A4S*TB)*TB
         DVS=1.d0+A3S*TB
         TWS=A1S*TTP**2
         DWS=1.d0+TWS
         DQS=(1.+TA)**4
         TZ=TVS/DVS/DWS/DQS
         SML=dlog(1.+TZ)
         SM=Sth+SML
* Derivatives for UM:
** over t=ln(TTP) -
         dTQ3dt=TQ3*2.5d0
         dQdt=-Q*dTQ3dt*A6/(1.d0+TQ3)
         dUdt=CVth-U
         dUMHdt=dUdt/(1.d0+Q)-U/(1.d0+Q)**2*dQdt
         dTVdt=(Zb*1.5d0*SqrtTB+2.d0*A4*TB)*TB
         dDVdt=dTVdt+A3*TB
         dDWdt=2.d0*TW
         DML=dTVdt/TV-dDVdt/DV-.5d0*dDWdt/DW
         dUMLdt=UML*DML
         dUMdt=dUMLdt+dUMHdt
         CMU=dUMdt+UM
** over b=ln(BP) -
         dQdb=-dQdt ! because Q=Q(TTP/BP):=Q(t-b)
         dUMHdb=-U/(1.d0+Q)**2*dQdb
         dTVdb=dTVdt
         dDVdb=dDVdt ! because DV=DV(TTP*BP):=DV(t+b)
         DMLb=dTVdb/TV-dDVdb/DV ! because DW=DW(TTP) => dDW/db=0
         dUMLdb=UML*DMLb
         dUMdb=dUMLdb+dUMHdb
* Derivatives for SM:
** over t=ln(TTP) -
         dTVSdt=Zb*(SqrtTB*2.5+A4S*TB*2.d0)*TB
         dDWSdt=2.d0*TWS
         DlnDQSdt=TA*4.d0/(1.+TA)
         DTZ=dTVSdt/TVS-A3S*TB/DVS-dDWSdt/DWS-DlnDQSdt
         dTZdt=TZ*DTZ
         CVL=dTZdt/(1.d0+TZ)
         CMS=CVL+CVth
** over b=ln(BP) -
         dTVSdb=dTVSdt
         DlnDQSdb=-DlnDQSdt
         DTZb=dTVSdb/TVS-A3S*TB/DVS-DlnDQSdb ! because dDWS/db=0
         dTZdb=TZ*DTZb
         dSMdb=dTZdb/(1.d0+TZ) ! because dSth/db=0
* Second derivatives:
** for U:
*** over tt
         dQdtt=-dQdt*dTQ3dt*A6/(1.d0+TQ3)+2.5d0*dQdt+
     +     Q*A6*(dTQ3dt/(1.d0+TQ3))**2
         dTVdtt=(Zb*2.25d0*SqrtTB+4.d0*A4*TB)*TB
         dDVdtt=dTVdtt+A3*TB
         dDWdtt=2.d0*dDWdt
         dDMLdt=dTVdtt/TV-dDVdtt/DV-.5d0*dDWdtt/DW-
     -     (dTVdt/TV)**2+(dDVdt/DV)**2+.5d0*(dDWdt/DW)**2
         dUMLdtt=dUMLdt*DML+dDMLdt*UML
         dUdtt=CWth-dUdt
         dUMHdtt=dUdtt/(1.d0+Q)-2.d0*dUdt/(1.d0+Q)**2*dQdt-
     -     U/(1.d0+Q)**2*dQdtt+2.d0*U/(1.d0+Q)**3*dQdt**2
         dUMdtt=dUMLdtt+dUMHdtt
         dCMUdt=dUMdtt+dUMdt
*** over bb
         dQdbb=dQdtt ! Q(t-b)
         dTVdbb=dTVdtt
         dDVdbb=dDVdtt ! DV(t+b)
         dDMLbdb=dTVdbb/TV-dDVdbb/DV-(dTVdb/TV)**2+(dDVdb/DV)**2
         dUMLdbb=dUMLdb*DMLb+dDMLbdb*UML
         dUMHdbb=-U/(1.d0+Q)**2*dQdbb+2.d0*U/(1.d0+Q)**3*dQdb**2
         dUMdbb=dUMLdbb+dUMHdbb
*** over tb
         dQdtb=-dQdtt ! Q(t-b)
         dTVdtb=dTVdtt
         dDVdtb=dDVdtt ! DV(t+b)
         dDMLbdt=dTVdtb/TV-dDVdtb/DV-dTVdt*dTVdb/TV**2+dDVdt*dDVdb/DV**2
         dUMLdtb=dUMLdt*DMLb+dDMLbdt*UML
         dUMHdtb=-(dUdt*dQdb+U*dQdtb)/(1.d0+Q)**2+
     +     2.d0*U/(1.d0+Q)**3*dQdb*dQdt
         dUMdtb=dUMLdtb+dUMHdtb
** for S:
*** over tt
         dTVSdtt=Zb*(SqrtTB*3.75d0+A4S*TB*4.d0)*TB
         dDWSdtt=4.d0*TWS
         dDTZdt=dTVSdtt/TVS-A3S*TB/DVS-dDWSdtt/DWS-DlnDQSdt-
     -  (dTVSdt/TVS)**2+(A3S*TB/DVS)**2+(dDWSdt/DWS)**2+DlnDQSdt**2/4.d0
         dTZdtt=dTZdt*DTZ+TZ*dDTZdt
         dCVLdt=dTZdtt/(1.d0+TZ)-(dTZdt/(1.d0+TZ))**2
         dCMSdt=dCVLdt+CWth
*** over bb
         dTVSdbb=dTVSdtt
         dDTZbdb=dTVSdbb/TVS-A3S*TB/DVS+DlnDQSdb-
     -     (dTVSdb/TVS)**2+(A3S*TB/DVS)**2+DlnDQSdb**2/4.d0
         dTZdbb=dTZdb*DTZb+TZ*dDTZbdb
         dSMdbb=dTZdbb/(1.d0+TZ)-(dTZdb/(1.d0+TZ))**2
*** over tb
         dTVSdtb=dTVSdtt
         dDTZbdt=dTVSdtb/TVS-A3S*TB/DVS-DlnDQSdb-
     -     dTVSdb*dTVSdt/TVS**2+(A3S*TB/DVS)**2-DlnDQSdb**2/4.d0
         dTZdtb=dTZdt*DTZb+TZ*dDTZbdt
         dSMdtb=dTZdtb/(1.d0+TZ)-dTZdt*dTZdb/(1.d0+TZ)**2
* Corrected values:
         SMth=SM+CMS-CMU
         UMth=UM+CMS-CMU
         CMth=CMS+dCMSdt-dCMUdt
         PMth=.5d0*(CMS+dSMdb-dUMdt-dUMdb)
         PDTth=PMth+.5d0*(dCMSdt+dSMdtb-dUMdtt-dUMdtb)
         PDRth=PMth+.25d0*(dUMdtt+2.d0*dUMdtb+dUMdbb-
     -     dCMSdt-2.d0*dSMdtb-dSMdbb)
* First moment of phonon frequencies and its variation:
         BC=CU1*BP*dsqrt(dsqrt(dsqrt(BP)))
         U1up=U1+BC*U1inf
         U1dn=1.d0+BC
         U1m=U1up/U1dn
         dLNu1=BC*1.125d0*(U1inf/U1up-1.d0/U1dn)
         d2LNu1=BC*CU1*1.265625d0*(U1*U1inf/U1up**2-1.d0/U1dn**2)
         B32=BP*dsqrt(BP)
         DU1=.0064*B32**2*dsqrt(B32)/((1.d0+B32)**2*dsqrt(1.d0+B32))
      endif
      return
      end

      subroutine FHARMAG(GAMI,TPT,BP,
     *     FMHL,UMHL,PMHL,CVMHL,SMHL,PDTMHL,PDRMHL,EVAR)
* Thermodynamic functions of a harmonic crystal in magnetic field
* 
*                                                       Version 01.06.12
* Input: GAMI - ionic Gamma, TPT=T_{p,i}/T
* Output: FMHL=F/(N_i T), UMHL=U/(N_i T), PMHL=P/(n_i T),
* CVMHL=C_V/N_i, SMHL=S/N_i
* PDTMHL = PMHL + d PMHL / d ln T, PDRMHL = PMHL + d PMHL/d ln\rho
* EVAR - energy variation with changing field direction
*  NB: zero-point cyclotron oscillations are excluded by U1m definition,
*  i.e., E0=E0(total)-\hbar\omega_c/(2kT)
      implicit double precision (A-H), double precision (O-Z)
      save
      parameter(CM=.895929256d0) ! Madelung
      call HLMAG(TPT,BP,UMth,CVMHL,SMHL,PMth,PDTMHL,PDRth,
     *  U1m,dLNu1,d2LNu1,DU1,UM,SM,CMU,CMS)
      U0=-CM*GAMI ! perfect lattice
      E0=1.5d0*U1m*TPT ! zero-point energy
      P0=.5d0*E0*(1.d0-dLNu1) ! zero-point pressure
      PDR0=E0*(.75d0-dLNu1+.25d0*dLNu1**2+.25d0*d2LNu1)
      EVAR=1.5d0*DU1*TPT
      UMph=UMth+E0 ! phonons: thermal + zero-point
      FMph=UMph-SMHL
      PMph=PMth+P0
      PDRph=PDRth+PDR0
      UMHL=U0+UMph
      FMHL=U0+FMph
      PMHL=U0/3.d0+PMph
      PDRMHL=U0/2.25d0+PDRph
      return
      end

      subroutine SPINION(GYRO,ZETI,MULTI,Fspin,Uspin,CVspin)
*                                                       Version 28.04.12
* Spin-degeneracy part of the entropy.
* Input: GYRO - gyrofactor
*        ZETI=\hbar\omega_c/kT
*        MULTI - spin multiplicity:
*                if MULTI=2, then assume fermions, otherwise bosons.
* Output: Fspin - spin free energy contribution / NkT
*         Uspin - spin internal energy contribution / NkT
*         CVspin - spin specific heat contribution / Nk
      implicit double precision (A-H), double precision (O-Z)
      parameter(TINY=1.d-7)
      CZ4=GYRO*ZETI*.25d0
      if (MULTI.gt.1.and.dabs(GYRO).gt.TINY) then ! add spin terms
        if (MULTI.eq.2) then
           CH=dcosh(GZ4)
           Fspin=-dlog(2.d0*CH)
           Uspin=-GZ4*dtanh(GZ4)
           CVspin=(GZ4/CH)**2
        else
           GZM4=GZ4*MULTI
           Fspin=-dlog(dsinh(GZM4)/dsinh(GZ4))
           Uspin=GZ4/dtanh(GZ4)-GZM4/dtanh(GZM4)
           CVspin=(GZ4/dsinh(GZ4))**2-(GZM4/dsinh(GZM4))**2
        endif
      else
         Fspin=0.
         Uspin=0.
         CVspin=0.
      endif
      return
      end

      subroutine QDMAG(TPT,ZETI,FWK,UWK,PWK,CVWK,PDTWK,PDRWK)
*                                                       Version 24.05.12
* Wigner-Kirkwood (quantum diffraction) corrections in a magnetic field
* Ref.: Alastuey A. & Jancovici B. (1980) Physica A 102, 327.
* Input: TPT - quantum plasma parameter (T_p/T)
*        ZETI=\hbar\omega_c/kT - cyclotron-to-thermal energy ratio
* Output: FWK, UWK, PWK - Wigner-Kirkwood contributions to the free and
*           internal energies and pressure, divuded by (n k T)
*         CVWK - contribution to the specific heat divided by (N k)
*         PDTWK=(1/nkT)*(d PWK/d ln T)_V
*         PDRWK=(1/nkT)*(d PWK/d ln\rho)_T
      implicit double precision (A-H), double precision (O-Z)
      if (ZETI.lt.0.) stop'WK-mag: ZETI < 0'
      if (ZETI.lt.1.d-2) then
         F1=-ZETI**2/11.25d0 ! u*f'(u)
         F0=1.d0+F1/2.d0 ! f(u)
         F2=F1 ! u^2*f''(u)
      elseif (ZETI.lt.20.) then
         CH=dcosh(ZETI)
         SH=dsinh(ZETI)
         F0=2.d0/ZETI*(CH/SH-1.d0/ZETI)+1.d0/3.d0 ! f(u)
         F1=4.d0/ZETI**2-2.d0/SH*(1.d0/SH+CH/ZETI) ! u*f'(u)
         F2=4.d0/SH*(CH/ZETI+ZETI*CH/SH**2+1.d0/SH)-12.d0/ZETI**2
      else ! u > 20, CH=SH\to\infty, CH/SH=1.
         F0=1.d0/3.d0+2.d0/ZETI-2.d0/ZETI**2 ! f(u)
         F1=4.d0/ZETI**2-2.d0/ZETI ! u*f'(u)
         F2=4.d0/ZETI-12.d0/ZETI**2 ! u^2*f''(u)
      endif
      WK=TPT**2/24.d0
      FWK=WK*F0
      UWK=WK*(2.d0*F0+F1)
      CVWK=-WK*(2.d0*F0+4.*F1+F2)
      PWK=FWK
      PDRWK=2.d0*PWK
      PDTWK=-WK*(F0+F1)
      return
      end
