* Effective nucleon masses.
* Contact: Alexander Potekhin <palex@astro.ioffe.ru>
* Last change: 27.08.13
      implicit double precision (A-H), double precision (O-Z)
      write(*,'('' Input BSk EOS number (19, 20, or 21): ''$)')
      read*,KEOS
      CMS=.8
      if (KEOS.eq.19) then
         CMV=.61
      elseif (KEOS.eq.20) then
         CMV=.65
      elseif (KEOS.eq.21) then
         CMV=.71
      else
         stop'unknown KEOS'
      endif
      write(*,112)
    1 continue
      write(*,'('' Input n_b [fm^{-3}] ( < 0 to stop): ''$)')
      read*,XN
      if(XN.le.0.) stop
      call FRACORE(KEOS,XN,Ye,Ymu)
      Yp=Ye+Ymu
      Yn=dim(1.d0,Yp)
      CMn=1.d0/(2.*Yn/CMS+(1.-2.*Yn)/CMV)
      CMp=1.d0/(2.*Yp/CMS+(1.-2.*Yp)/CMV)
      call EFFMASS(KEOS,XN,Yp,EFMp,EFMn)
      write(*,111) Yp,EFMp,EFMn
      goto 1
  111 format(1P10E11.3)
  112 format('Output of meff v.11.04.2013: effective nucleon masses'/
     /  '    Y_p        m*_p      m*_n')
      end

      subroutine EFFMASS(KEOS,XN,Yp,EFMp,EFMn)
*                                                       Version 11.04.13
* Proton and neutron effective masses in the nuclear matter.
* Calculation according to Eq.(A10) of Chamel, Goriely, & Pearson (2009)
*     with parameters from Table IV of Goriely, Chamel, & Pearson (2010)
* The code by A.Y.Potekhin <palex@astro.ioffe.ru>
* Input: KEOS - number of BSK model (19,20,or 21)
*        XN = n_b - baryon number density in fm^{-3}
*        Yp - proton fraction
* Output: EFMp, EFMn - ratios of the proton and neutron effective masses
*        to the masses of an isolated proton and neutron, respectively.
      implicit double precision (A-H), double precision (O-Z)
      save AT1,AT2X2,AT4,AT5,AX1,AX4,AX5
      dimension AT1(3),AT2X2(3),AT4(3),AT5(3),AX1(3),AX4(3),AX5(3)
      parameter (HbarC=197.3269718) ! \hbar*c [Mev*fm]
      parameter (Ep=938.272046,En=939.565379) ! m_p*c^2, m_n*c^2 [MeV]
      data AT1/403.072,438.219,396.131/,
     *  AT2X2/-1055.55,-1147.64,-1390.38/,
     *  AT4/-60.,-100.,-100./
     *  AT5/-90.,-120.,-150./,
     *  AX1/-0.137960,-0.392047,-0.0648452/,
     *  AX4/-6.,-3.,2./,
     *  AX5/-13.,-11.,-11./
      K=KEOS-18
      Yn=1.d0-Yp
      if (K.eq.1) then
         BETA=1.d0/3.d0
      elseif (K.eq.2) then
         BETA=1.d0/6.d0
      elseif (K.eq.3) then
         BETA=.5d0
      else
        stop'EFFMASS: unknown EOS'
      endif
      GAMMA=1.d0/12.d0
      DFMp=AT1(K)*((1.d0+.5d0*AX1(K))-(.5d0+AX1(K))*Yp)+
     +  AT2X2(K)*(.5d0+Yp)+
     +  AT4(K)*((1.d0+.5d0*AX4(K))-(.5d0+AX4(K))*Yp)*XN**BETA+
     +  AT5(K)*((1.d0+.5d0*AX5(K))+(.5d0+AX5(K))*Yp)*XN**GAMMA
      DFMn=AT1(K)*((1.d0+.5d0*AX1(K))-(.5d0+AX1(K))*Yn)+
     +  AT2X2(K)*(.5d0+Yn)+
     +  AT4(K)*((1.d0+.5d0*AX4(K))-(.5d0+AX4(K))*Yn)*XN**BETA+
     +  AT5(K)*((1.d0+.5d0*AX5(K))+(.5d0+AX5(K))*Yn)*XN**GAMMA
      EFMp=1.d0/(1.d0+.5d0*DFMp*XN*Ep/HbarC**2)
      EFMn=1.d0/(1.d0+.5d0*DFMn*XN*En/HbarC**2)
      return
      end