* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * **
*    ANALYTICAL REPRESENTATIONS OF NEUTRON-STAR EQUATIONS OF STATE     *
*    This is a supplement to bskfit.f: auxiliary driving routines      *
*            ( !!! TO BE LINKED WITH bskfit.f !!! )                    *
*      Remarks and suggestions are welcome. Please send them to        *
*         Alexander Potekhin <palex@astro.ioffe.ru>                    *
*   For theoretical background and references see                      *
*           http://www.ioffe.ru/astro/NSG/BSk/                         *
*   Last update: 20.07.12                                              *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * **

*   ----------------------   MAIN block   ---------------------------  *
*       This is auxiliary MAIN program for input/output purposes.      *
*                You can change it or write your own.                  *
*  Alternatively, you can just delete this MAIN block and link the set *
*           of the remaining subroutines with yours.                   *
*     Calculations are performed in the subroutine NSEOSFIT (below).   *
*   -----------------------------------------------------------------  *
      implicit double precision (a-h,o-z)
      write(*,'('' Arguments:''/
     *  ''  (1) number density of baryons n (in fm^{-3}),''/
     *  ''  (2) mass density of baryons rho (in g/cc),''/
     *  ''  (3) excess enthalpy per baryon h1=h/h0-1,''
     *  '' where h0=m_0*c^2,''/
     *  ''  (4) pressure P (in dyn/cm^2),''/
     *  ''  (5) adiabatic exponent Gamma.''/
     *  '' For one input argument (1, 2, or 3), other '',
     *  ''arguments are fitted.''/
     *  '' Input mode: 1 (n), 2 (rho), or 3 (h1) ? ''$)')
      read*,MODE
    1 write(*,'('' Enter EOS number: 19, 20, or 21'',
     *  '' (for BSk19, BSk20, BSk21) [enter 0 to stop] ? ''$)')
      read*,KEOS
      if (KEOS.eq.0) stop
      if (KEOS.ne.19.and.KEOS.ne.20.and.KEOS.ne.21) goto 1
      write(*,110)
      IRUN=0
   10 IRUN=IRUN+1
      if (MODE.eq.1) then
         write(*,'('' n''$)')
      elseif (MODE.eq.2) then
         write(*,'('' rho''$)')
      elseif (MODE.eq.3) then
         write(*,'('' h1''$)')
      endif
      if (IRUN.eq.1) write(*,'(''  ( < 0 to terminate)''$)')
      write(*,'('': ''$)')
      if (MODE.eq.1) then
         read*,XN
        if (XN.le.0.) goto 1
      elseif (MODE.eq.2) then
         read*,RHO
        if (RHO.le.0.) goto 1
      elseif (MODE.eq.3) then
         read*,H1
        if (H1.le.0.) goto 1
      else
         stop'unknown MODE'
      endif
      call BSKEOS(MODE,KEOS,XN,RHO,H1,P,Gamma)
      write (*,111) MODE,KEOS,XN,RHO,H1,P,Gamma
      goto 10
  110 format(' mode EOS#   n         rho        h1',9X,'P',9X,'Gamma')
  111 format(2I4,1P10E11.3)
      end
*   ----------------  END of MAIN block   ---------------------------  *
      
      subroutine BSKEOS(MODE,KEOS,XN,RHO,H1,P,Gamma)
*                                                       Version 20.07.12
* Neutron-star EOS fitting.
* Arguments:
* MODE regulates input/output (see below),
* KEOS=19, 20, 21 for BSk19, BSk20, BSk21, respectively,
* XN  - number density of baryons n (in fm^{-3}),
* RHO - mass density of baryons rho (in g/cc),
* H1  - excess enthalpy per baryon h1=h/h0-1, where h0=m_0*c^2,
* P   - pressure (in dyn/cm^2),
* Gamma - adiabatic exponent, Gamma= d ln P / d ln n.
* If MODE=1, 2, or 3, then an input argument is XN, RHO, or H1, resp.;
*    the other 4 arguments are fitted on output.
      implicit double precision (a-h,o-z)
      parameter(c2=2.99792458d10**2)
      if (MODE.eq.1) then
         call BSkRofN(KEOS,XN,RHO)
      elseif (MODE.eq.2) then
         call BSkNofR(KEOS,RHO,XN)
      elseif (MODE.eq.3) then
         call BSKfitH(KEOS,H1,RLG)
         RHO=10.d0**RLG
      else
         stop'BSkEOS: unknown MODE'
      endif
      RLG=dlog10(RHO)
      call BSKfit(KEOS,RLG,PLG,CHIrho,Gamma)
      P=10.**PLG+3.5d14*RHO
      if (MODE.eq.3) then
         XN=(RHO+P/c2)/(H1+1.d0)/1.66d15
      else
         H1=(RHO+P/c2)/XN/1.66d15-1.d0
      endif
      return
      end