* Fitting functions reproducing EOSs BSk19,20,21.
* Last change: 17.02.2013.
* Contact: Alexander Potekhin <palex@astro.ioffe.ru>
* List of subroutines:
*   BSKfit returns pressure P and log.derivative d ln P / d ln\rho
*     at a given mass density \rho (at \rho > 10^6 g/cc);
*   BSKfitH returns \rho (at \rho > 10^6 g/cc) as a function of
*     H1 = e^H-1, where H is the dimensinoless pseudo-enthalpy;
*   BSkRofN returns mass density \rho as a function of n_b;
*   BSkNofR returns baryon number density n_b as a function of \rho;
*   FRACORE returns fractions of electrons Y_e and muons Y_\mu
*     as functions of n_b in the stellar core;
*   INCRUST returns the number of bound nucleons in a nucleus A,
*     the total number of nucleons per one nucleus A',
*     nuclear-shape parameters a_p, a_n, C_p, C_n, 
*     as well as additional parameters n_{0p} and n_{0n} - the heights
*     of the bumps of proton and neutron densities near the center of
*     a Wigner-Seitz cell, and the ratios xnuc and xnucn
*     of the effective proton and neutron radii of these bumps,
*     respectively, to the cell radius.
 
      subroutine BSKfit(KEOS,RLG,PLG,CHIrho,Gamma)
* Fit of P(\rho) for BSk19,20,21 at rho > 10^6 g/cc
* average error 1.%; max.error 3.7% at RLG=9.51
* Input: KEOS - number of BSK model (19,20,or 21)
*        RLG = lg(rho[g/cc])
* Output: PLG = lg(P[dyn/cm^2]),
*         CHIrho = d log P / d log rho
*         Gamma = d log P / d log n_b
*                                                          Version 12/11
      implicit double precision (A-H), double precision (O-Z)
      parameter(c2=2.99792458d10**2)
      dimension A(23,3)
      save A
      data A/
* BSk19:
     *  3.916,7.701,.00858,.22114,3.269,11.964,13.349,1.3683,
     *    3.254,-12.953,0.9237,6.20,14.383,16.693,-1.0514,2.486,
     *    15.362,.085,11.68,6.23,-.029,14.19,20.1,
* BSk20:
     *  4.078,7.587,.00839,.21695,3.614,11.942,13.751,1.3373,
     *    3.606,-22.996,1.6229,4.88,14.274,23.560,-1.5564,2.095,
     *    15.294,.084,11.67,6.36,-.042,14.18,14.8,
* BSk21:
     *  4.857,6.981,.00706,.19351,4.085,12.065,10.521,1.5905,
     *    4.104,-28.726,2.0845,4.89,14.302,22.881,-1.7690,0.989,
     *    15.313,.091,11.65,4.68,-.086,14.15,10.0/
      if (RLG.lt.-1.) stop'BSKfit: too low RLG'
      if (RLG.gt.16.5) stop'BSKfit: too high RLG'
      K=KEOS-18
      if (K.lt.1.or.K.gt.3) stop'BSKfit: unknown EOS'
      X=RLG
      D4=1.d0+A(4,K)*X
      P01=(A(1,K)+A(2,K)*X+A(3,K)*X**3)/D4
      F1=FERMI(A(5,K)*(X-A(6,K)))
      P02=A(7,K)+A(8,K)*X
      F2=FERMI(A(9,K)*(A(6,K)-X))
      P03=A(10,K)+A(11,K)*X
      F3=FERMI(A(12,K)*(A(13,K)-X))
      P04=A(14,K)+A(15,K)*X
      F4=FERMI(A(16,K)*(A(17,K)-X))
      G05=(X-A(19,K))*A(20,K)
      P05=A(18,K)/(1.d0+G05**2)
      G06=(X-A(22,K))*A(23,K)
      P06=A(21,K)/(1.d0+G06**2)
      PLG=P01*F1+P02*F2+P03*F3+P04*F4+P05+P06
      CHIrho=F1/D4* ! d\zeta/dX=d ln P/d ln rho
     *  ((A(2,K)-A(1,K)*A(4,K)+3.*A(3,K)*X**2+
     +    2.*A(3,K)*A(4,K)*X**3)/D4-
     -  A(5,K)*(A(1,K)+A(2,K)*X+A(3,K)*X**3)*(1.d0-F1)) +
     +  F2*(A(8,K)+A(9,K)*(1.d0-F2)*P02) +
     +  F3*(A(11,K)+A(12,K)*(1.d0-F3)*P03)+
     +  F4*(A(15,K)+A(16,K)*(1.d0-F4)*P04)-
     -  2.d0*(P05*G05*A(20,K)/(1.d0+G05**2)+
     +        P06*G06*A(23,K)/(1.d0+G06**2))
      Gamma=(1.d0+10.d0**(PLG-RLG)/c2)*CHIrho
      return
      end

      subroutine BSKfitH(KEOS,H1,RLG)
* Fit of rho(H) for BSk19,20,21 at rho > 10^6 g/cc (i.e.,e^H-1 > 5.9e-5)
* average error 1.3%; max.error < 5% (at both inner-crust boundaries)
* Input: KEOS - number of BSK model (19,20,or 21)
*        H1 = e^H-1, where H = pseudo-enthalpy [Eq.(7) of Haens.Pot.04]
* Output: RLG = lg(rho[g/cc])
*                                                          Version 02/12
      implicit double precision (A-H), double precision (O-Z)
      save A
      dimension A(10,3)
      data A/
     *  93.650,36.893,15.450,.672,61.240,68.97,.292,5.2,.48,6.8, ! BSk19
     *  81.631,31.583,15.310,.594,58.890,56.74,.449,4.5,.58,7.5, ! BSk20
     *  63.150,23.484,15.226,.571,54.174,37.15,.596,3.6,.51,10.4/! BSk21
      X=dlog10(H1)
      K=KEOS-18
      if (K.lt.1.or.K.gt.3) stop'BSKfitH: unknown EOS'
      RLG1=2.367+21.84*H1**.1843/(1.+.7*H1)
      RLG3=(A(3,K)+A(4,K)*X)
      G=(A(5,K)*H1)**A(10,K)
      RLG2=(A(1,K)+A(2,K)*X+G*RLG3)/(1.+A(6,K)*H1+G)
      Y=84.*(X+1.99)
      FY=FERMI(Y)
      RLG=RLG1*FY+RLG2*dim(1.d0,FY)+FERMI(A(8,K)*(A(9,K)-X))*A(7,K)
      return
      end

      subroutine BSkRofN(KEOS,XN,RHO)
*                                                          Version 02/12
* Fit of (rho(n_b)-m_0*n_b) consistent with BSKfit
*   for BSk19,20,21 at rho > 10^6 g/cc (i.e., n_b > 6.e-10 fm^{-3})
* Criterion is F = \rho/(m_0 n_b)-1 = e^H - 1 - P/(m_0 n_b c^2);
* average error in F is 1.7%; max.error 4.2% at rho=10^6 g/cc = min(rho)
* Input: KEOS - number of BSK model (19,20,or 21)
*        XN = n_b - baryon number density in fm^{-3}
* Output: RHO - mass-energy density in [g/cc]
      implicit double precision (A-H), double precision (O-Z)
      save A
      dimension A(8,3)
      data A/
     *  .259,2.30,.0339,.1527,3.170,1.085E-05,.6165,3.50, ! BSk19
     *  .632,2.71,.0352,.383,3.165,1.087E-05,.6167,3.51,  ! BSk20
     *  3.85,3.19,.0436,1.99,3.210,1.075E-05,.6154,3.44/  ! BSk21
      K=KEOS-18
      if (K.lt.1.or.K.gt.3) stop'BSkRofN: unknown EOS'
      F1=(A(1,K)*XN**A(2,K)+A(3,K)*dsqrt(XN))/(1.d0+A(4,K)*XN)**2
      F2=XN/(A(6,K)+XN**A(7,K)*A(8,K))
      XLG=dlog10(XN)
      FRM=FERMI(1.1*(XLG+A(5,K)))
      F=FRM*F2+(1.d0-FRM)*F1
      RHO=XN*1.66d15*(1.d0+F)
      return
      end

      subroutine BSkNofR(KEOS,RHO,XN)
*                                                          Version 02/12
* Fit of (rho-m_0*n_b(rho)) consistent with BSKfit at rho > 10^6 g/cc.
* Criterion is F = \rho/(m_0 n_b)-1 = e^H - 1 - P/(m_0 n_b c^2);
* average error in F is 1.8% for BSk19, 1.5% for BSk20 and 21,
*   max.error in F is 3.9% for BSk19,20 and 3.7% for BSk21 (at min\rho)
* Input: KEOS - number of BSK model (19,20,or 21)
*        RHO - mass-energy density in [g/cc]
* Output: XN = n_b - baryon number density in fm^{-3}
      implicit double precision (A-H), double precision (O-Z)
      save A
      dimension A(9,3)
      data A/
* BSk19:
     *  .378,1.28,17.20,3.844,3.778,12.0741,1.071E-05,3.287,.613,
* BSk20:
     *  .152,1.02,10.26,3.691,2.586,12.0570,1.067E-05,3.255,.6123,
* BSk21:
     *  .085,.802,16.35,3.634,2.931,12.0190,1.050E-05,3.187,.611/
      K=KEOS-18
      if (K.lt.1.or.K.gt.3) stop'BSkNofR: unknown EOS'
      X=RHO/1.66d15
      F1=(A(1,K)*X**A(2,K)+A(3,K)*X**A(4,K))/(1.d0+A(5,K)*X)**3
      F2=X/(A(7,K)+A(8,K)*X**A(9,K))
      RLG=dlog10(RHO)
      FRM=FERMI(RLG-A(6,K))
      F=FRM*F2+(1.-FRM)*F1
      XN=X/(1.d0+F)     
      return
      end

      function FERMI(X)
      implicit double precision (A-H), double precision (O-Z)
      if (X.gt.40.d0) then
         F=0.d0
         goto 50
      endif
      if (X.lt.-40.d0) then
         F=1.d0
         goto 50
      endif
      F=1.d0/(dexp(X)+1.d0)
   50 FERMI=F
      return
      end

      subroutine FRACORE(KEOS,XN,Ye,Ymu)
*                                                          Version 05/12
* Fit of particle fractions in the core.
* Input: KEOS - number of BSK model (19,20,or 21)
*        XN = n_b - baryon number density in fm^{-3}
* Output: Ye - number of electrons divided by number of nucleons
*         Ymu - number of muons divided by number of nucleons
* NB: proton fraction is Yp=Ye+Ymu, neutron fraction is Yn=1-Yp.
* average error 0.00011 (Ye, BSk19), 0.00013 (Ymu, BSk19),
*               0.00010 (Ye, BSk20), 0.00012 (Ymu, BSk20),
*               0.00023 (Ye, BSk21), 0.00026 (Ymu, BSk21),
*     max.error 0.00024 (Ye, BSk19), 0.00034 (Ymu, BSk19),
*               0.00020 (Ye, BSk20), 0.00070 (Ymu, BSk20),
*               0.00042 (Ye, BSk21), 0.00058 (Ymu, BSk21),
      implicit double precision (A-H), double precision (O-Z)
      save A
      dimension A(12,3)
      data A/
* BSk19:
     *  -.0157,.9063,0.,26.97,106.5,4.82,
     *  -.0315,.25,0.,12.42,72.4,19.5,
* BSk20:
     *  -.0078,.745,.508,22.888,.449,.00323,
     *  -.0364,.2748,.2603,12.99,.0767,.00413,
* BSk21:
     *  .00575,.4983,9.673,16.31,38.364,0.,
     *  -.0365,.247,11.49,24.55,48.544,0./
      K=KEOS-18
      if (K.lt.1.or.K.gt.3) stop'FRACORE: unknown EOS'
      Ye=(A(1,K)+A(2,K)*XN+A(3,K)*XN**4)/
     /  (1.+A(4,K)*XN*dsqrt(XN)+A(5,K)*XN**4)*
     *  exp(-A(6,K)*XN**5)
      if (Ye.lt.0.) Ye=0.
      Ymu=(A(7,K)+A(8,K)*XN+A(9,K)*XN**4)/
     /  (1.+A(10,K)*XN*dsqrt(XN)+A(11,K)*XN**4)*
     *  exp(-A(12,K)*XN**5)
      if (Ymu.lt.0.) Ymu=0.
      return
      end

      subroutine INCRUST(KEOS,XN,XND,XNC,CMI,CMI1,Zcell,Zclust,
     *  ap,an,Cp,Cn,Dn_p,Dn_n,xnuc,xnucn)
*                                                       Version 17.02.13
* Fit of the nuclear-shape parameters.
* Input: KEOS - number of BSK model (19,20,or 21)
*        XN = n_b - baryon number density in fm^{-3}
* Output: XND - neutron-drip number density of baryons in fm^{-3} 
*              = the inner/outer crust interface (Pearson et al. 2012),
*         XNC - number density of baryons in fm^{-3} at the crust/core
*              interface (Pearson et al. 2012),
*         CMI - average number of bound nucleons in a nucleus,
*         CMI1 - average total number of nucleons per 1 nucleus,
*         Zcell - total number of protons per a nucleus,
*         Zclust - number of clusterized protons in a nucleus,
*         a_q,C_q,Dn_q (q=p,n) - shape parameters (Onsi et al. 2008),
*         xnuc, xnucn - proton and neutron equivalent size parameters,
*          defined as (Gnedin et al. 2003) xnuc = sqrt{5<r^2>/3}
* maximal fractional errors [measured for n_b up to max(n_b[fm^{-3}])]:
* BSk:          19                    20                   21   
* ap:     1.1%  [.065]            1.2%  [.066]           0.34% [.062]
* an:  5.9%[.082]/1.0%[.047] 6.0%[.081]/0.85%[.046]       2.0% [.065]
* Cp: 0.43%[.048]/3%[.0815]    1%[.055]/3.3%[.08]    1%[.052]/9%[.078]
* Cn: 1%[.053]/3%[.0815] 1%[.056]/3%[.0795] 0.7%[.061]/2%[.075]/7%[.082]
* Dn_p: 1%[.056]/7%[.0745]     1%[.059]/9%[.0815]    2%[.056]/9%[.078]
* Dn_n: 0.7% [.065]/7%[.081] 0.6%[.066]/7%[.079]   0.4%[.065]/7%[.079]
* xnuc:    5% [.0815]               4% [.080]      2.5%[.065]/6.3%[.078]
* xnucn:  4.6% [.086]               4% [.084]              2.5% [.082]
* CMI:  0.8%[.063]/7%[.081]  0.7%[.063]/6%[.079]   0.8%[.059]/8%[.079]
* CMI1: 0.6%[.054]/4%[.070]        1.5% [0.1]      0.6%[.065]/2%[.079]
* The listed max.errors are reached near the lower and/or higher ends
* of the XN range. The second (larger) values (after slash) are probably
* due to inaccuracies in the numerical data.
* Typical errors are 10 times smaller than the first (lower) max.errors.
      implicit double precision (A-H), double precision (O-Z)
* Fit parameters for shape parameters a_q,C_q,xnuc_q,dn_q (q=p,n):
      save PARap,PARan,PARCp,PARCn,PARxp,PARxn,PARdp,PARdn,PARNi,A1,LZ
      dimension PARap(3,3),PARan(4,3),PARCp(4,3),PARCn(4,3),
     *  PARxp(4,3),PARxn(4,3),PARdp(4,3),PARdn(3,3),
     *  PARNi(6,3),A1(6,3),LZ(2,3)
      data PARap/
     *  .4377,4.360,1084., ! BSk19
     *  .4353,4.440,1154., ! BSk20
     *  .4316,4.704,1253./ ! BSk21
      data PARan/
     *  .639,1.461,.457,1137., ! BSk19
     *  .632,1.98,0.514,1122., ! BSk20
     *  .636,5.32,0.739,624./  ! BSk21
      data PARCp/
     *  5.500,11.7,.643,472., ! BSk19
     *  5.493,12.8,.636,484., ! BSk20
     *  5.457,14.2,.601,566./ ! BSk21
      data PARCn/
     *  5.714,14.05,.642,182., ! BSk19
     *  5.714,16.3,0.645,175., ! BSk20
     *  5.728,22.2,0.663,144./ ! BSk21
      data PARxp/
     *  .1120,2.06,.633,507., ! BSk19
     *  .1094,2.04,.613,509., ! BSk20
     *  .1045,2.09,.586,513./ ! BSk21
      data PARxn/
     *  .122,2.27,.618,193., ! BSk19
     *  .119,2.30,.603,182., ! BSk20
     *  .114,2.56,.595,107./ ! BSk21
      data PARdp/
     *  .05509,.1589,  0., .4917, ! BSk19
     *  .05382,.1400,.01715,.566, ! BSk20
     *  .05273,.1107,.02218,.6872/ ! BSk21
      data PARdn/
     *  .10336,.0772,1.129, ! BSk19
     *  .10283,.0825,1.189, ! BSk20
     *  .10085,.0942,1.279/ ! BSk21
      data PARNi/ ! for number of bound neutrons
     *  93.0,11.90,1.490,.334,5.05, 8.40, ! BSk19
     *  92.8,12.95,1.493,.354,7.57, 9.30, ! BSk20
     *  92.3,13.80,1.625,.3874,13.8,10.8/ ! BSk21
      data A1/
     *  134.7,183.7,308.7,.3814,-.00058,.00049, ! BSk19
     *  134.7,188.2,275.6,.4346,.00163,.00149,  ! BSk20
     *  132.6,187.6,229.2,.5202,.00637,.00151/  ! BSk21
      data LZ/19,16, 20,19, 27,17/
      K=KEOS-18
      if (KEOS.eq.19) then
         XND=2.63464d-4
         XNC=.0885
      elseif (KEOS.eq.20) then
         XND=2.62873d-4
         XNC=.0854
      elseif (KEOS.eq.21) then
         XND=2.57541d-4
         XNC=.0809
      else
         stop'INCRUST: unknown EOS'
      endif
      X=XN/XND
      Xlg=dmax1(dlog10(X),0.d0)
      XD=XN/XNC
      Zcell=40.d0
      Zclust=Zcell*dexp(-XD**LZ(1,K))*dim(1.d0,XD**LZ(2,K))
      P3Xi=PARNi(3,K)*Xlg
      CNI=(PARNi(1,K)+PARNi(2,K)*Xlg+P3Xi**3*dsqrt(P3Xi))/
     /  (1.d0+(PARNi(4,K)*Xlg)**PARNi(5,K))*
     *  (dim(1.d0,XD**PARNi(6,K)))
      CMI=Zclust+CNI
      CMI1=(A1(1,K)+A1(2,K)*Xlg+A1(3,K)*Xlg**2)/(1+(A1(4,K)*Xlg)**4)*
     *  (1+A1(5,K)*X)*dim(1.d0,(A1(6,K)*X)**2)
      ap=(PARap(1,K)+PARap(2,K)*XN)/(1.d0-PARap(3,K)*XN**3)
      an=(PARan(1,K)+PARan(2,K)*XN**PARan(3,K))/(1.d0-PARan(4,K)*XN**3)
      Cp=(PARCp(1,K)+PARCp(2,K)*XN**PARCp(3,K))/(1.d0-PARCp(4,K)*XN**3)
      Cn=(PARCn(1,K)+PARCn(2,K)*XN**PARCn(3,K))/(1.d0-PARCn(4,K)*XN**3)
      xnuc=(PARxp(1,K)+PARxp(2,K)*XN**PARxp(3,K))/
     /  (1.d0-PARxp(4,K)*XN**3)
      xnucn=(PARxn(1,K)+PARxn(2,K)*XN**PARxn(3,K))/
     /  (1.d0-PARxn(4,K)*XN**3)
      Dn_p=PARdp(1,K)-PARdp(2,K)*XN/(PARdp(3,K)+XN**PARdp(4,K))
      Dn_p=Dn_p*dim(1.d0,XD**9)
      Dn_n=PARdn(1,K)+PARdn(2,K)*XN**.5-PARdn(3,K)*XN
      Dn_n=Dn_n*dim(1.d0,XD**16)
      return
      end
