* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ** * ANALYTICAL REPRESENTATIONS OF NEUTRON-STAR EQUATIONS OF STATE * * Remarks and suggestions are welcome. Please send them to * * Alexander Potekhin * * For theoretical background and references see * * P.Haensel & A.Y. Potekhin, A&A, *** (submitted) *** (2004) * * and http://www.ioffe.ru/astro/NSG/NSEOS/ * * Last update: 27.08.04 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ** * ---------------------- 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). * * ----------------------------------------------------------------- * C%C implicit double precision (a-h,o-z) C%C write(*,'('' Arguments:''/ C%C * '' (1) number density of baryons n (in fm^{-3}),''/ C%C * '' (2) mass density of baryons rho (in g/cc),''/ C%C * '' (3) excess enthalpy per baryon h1=h/h0-1,'' C%C * '' where h0=m_0*c^2,''/ C%C * '' (4) pressure P (in dyn/cm^2),''/ C%C * '' (5) adiabatic exponent Gamma.''/ C%C * '' For one input argument (1, 2, or 3), other '', C%C * ''arguments are fitted.''/ C%C * '' Input mode: 1 (n), 2 (rho), or 3 (h1) ? ''$)') C%C read*,MODE C%C 1 write(*,'('' EOS: 1 (SLy), 2 (FPS), or 3 (SLy+APR)'', C%C * '' [enter 0 to stop] ? ''$)') C%C read*,KEOS C%C if (KEOS.eq.0) stop C%C if (KEOS.lt.0.or.KEOS.gt.3) goto 1 C%C write(*,110) C%C IRUN=0 C%C 10 IRUN=IRUN+1 C%C if (MODE.eq.1) then C%C write(*,'('' n''$)') C%C elseif (MODE.eq.2) then C%C write(*,'('' rho''$)') C%C elseif (MODE.eq.3) then C%C write(*,'('' h1''$)') C%C endif C%C if (IRUN.eq.1) write(*,'('' ( < 0 to terminate)''$)') C%C write(*,'('': ''$)') C%C if (MODE.eq.1) then C%C read*,XN C%C if (XN.le.0.) goto 1 C%C elseif (MODE.eq.2) then C%C read*,RHO C%C if (RHO.le.0.) goto 1 C%C elseif (MODE.eq.3) then C%C read*,H1 C%C if (H1.le.0.) goto 1 C%C else C%C stop'unknown MODE' C%C endif C%C call NSEOSFIT(MODE,KEOS,XN,RHO,H1,P,Gamma) C%C write (*,111) MODE,KEOS,XN,RHO,H1,P,Gamma C%C goto 10 C%C 110 format(' mode EOS# n rho h1',9X,'P',9X,'Gamma') C%C 111 format(2I4,1P10E11.3) C%C end * ---------------- END of MAIN block --------------------------- * subroutine NSEOSFIT(MODE,KEOS,XN,RHO,H1,P,Gamma) * Version 07.04.04 * Neutron-star EOS fitting. * Arguments: * MODE regulates input/output (see below), * KEOS=1 for SLy EOS, 2 for FPS EOS, * 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, respectively; * 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 FitRofN(KEOS,XN,RHO) elseif (MODE.eq.2) then call FitNofR(KEOS,RHO,XN) elseif (MODE.eq.3) then call FitHEOS(KEOS,H1,RLG) RHO=10.d0**RLG else stop'NSEOSFIT: unknown MODE' endif if (KEOS.lt.1.or.KEOS.gt.3) stop'NSEOSFIT: unknown EOS' RLG=dlog10(RHO) call FitGammaP(KEOS,RLG,PLG,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 subroutine FitRofN(KEOS,XN,RHO) * Version 07.04.04 * Fit of (rho(n_b)-m_0*n_b) consistent with FitEOS to within * 2.1/5.4% for SLy4, 2.0/6.6% for FPS, 3.1/7.6% for SLy+APR * Error in rho within 0.47% for SLy4, 0.35% for FPS, 0.8% for SLy+APR * Input: KEOS; XN=n_b(fm^{-3}) * Output: RHO = rho[g/cc] * KEOS=1: SLy, KEOS=2: FPS, KEOS=3: SLy+APR implicit double precision (A-H), double precision (O-Z) dimension A(7,3) data A/.423,2.42,.031, .78,.238,.912,3.674, * .320,2.17,.173,3.01,.540,.847,3.581, * .342,2.23,2.2d-6,10.92,0.,.839,3.68/ if (KEOS.lt.1.or.KEOS.gt.3) stop'FitEOS: KEOS' if (XN.lt.1.d-16) stop'FitRofN: too low XN' if (XN.gt.5.) stop'FitRofN: too high XN' F1=(A(1,KEOS)*XN**A(2,KEOS)+A(3,KEOS)*XN**A(4,KEOS))/ / (1.d0+A(5,KEOS)*XN)**2 F2=XN/(8.d-6+XN**.585*2.1) XLG=dlog10(XN) FRM=FERMI(A(6,KEOS)*(XLG+A(7,KEOS))) F=FRM*F2+(1.-FRM)*F1 RHO=XN*1.66d15*(1.d0+F) return end subroutine FitNofR(KEOS,RHO,XN) * Version 07.04.04 * Fit of (rho-m_0*n_b(rho)) consistent with FitEOS to within * 2.1/5.4% for SLy4, 2.1/6.6% for FPS, 2.3/7.3% for SLy+APR * Error in n_b within .56% for SLy, .74% for FPS, 1.4% for SLy+APR * Input: KEOS; RHO = rho[g/cc] * Output: XN=n_b(fm^{-3}) * KEOS=1: SLy4, KEOS=2: FPS, KEOS=3: SLy+APR implicit double precision (A-H), double precision (O-Z) dimension A(7,3) save data A/.183,1.26,6.88,3.612,2.248,.911,11.56, * .608,2.41,2.39,3.581,1.681,.850,11.64, * .173,1.18,9.97,3.787,2.634,.917,11.56/ if (KEOS.lt.1.or.KEOS.gt.3) stop'FitNofR: KEOS' if (RHO.lt..1) stop'FitNofR: too low RHO' if (RHO.gt.1.d16) stop'FitNofR: too high RHO' X=RHO/1.66d15 F1=(A(1,KEOS)*X**A(2,KEOS)+A(3,KEOS)*X**A(4,KEOS))/ / (1.d0+A(5,KEOS)*X)**3 F2=X/(8.d-6+2.1*X**.585) RLG=dlog10(RHO) FRM=FERMI(A(6,KEOS)*(RLG-A(7,KEOS))) F=FRM*F2+(1.-FRM)*F1 XN=X/(1.d0+F) return end subroutine FitGammaP(KEOS,RLG,PLG,Gamma) * Version 05.04.04 * Adiabatic index - formula consistent with FitEOS * Input: KEOS; RLG = lg(rho[g/cc]) * Output: PLG=lg(P[dyn/cm^2]), Gamma=d log P / d log n_b * KEOS=1: SLy4 * KEOS=2: FPS * KEOS=3: SLy4+APR implicit double precision (A-H), double precision (O-Z) save parameter(c2=2.99792458d10**2) dimension A(18,3) data A/6.22,6.121,.005925,.16326,6.48,11.4971,19.105,.8938, *6.54,11.4950,-22.775,1.5707,4.3,14.08,27.80,-1.653,1.50,14.67, * 6.22,6.121,.006004,.16345,6.50,11.8440,17.24,1.065, *6.54,11.8421,-22.003,1.5552,9.3,14.19,23.73,-1.508,1.79,15.13, * 6.22,6.121,.006035,.16354,4.73,11.5831,12.589,1.4365, *4.75,11.5756,-42.489,3.8175,2.3,14.81,29.80,-2.976,1.99,14.93/ if (KEOS.lt.1.or.KEOS.gt.3) stop'FitEOS: KEOS' if (RLG.lt.-1.) stop'FitGammaP: too low RLG' if (RLG.gt.16.) stop'FitGammaP: too high RLG' X=RLG D4=(1.+A(4,KEOS)*X) P01=(A(1,KEOS)+A(2,KEOS)*X+A(3,KEOS)*X**3)/D4 P02=A(7,KEOS)+A(8,KEOS)*X P03=A(11,KEOS)+A(12,KEOS)*X P04=A(15,KEOS)+A(16,KEOS)*X F1=FERMI(A(5,KEOS)*(X-A(6,KEOS))) F2=FERMI(A(9,KEOS)*(A(10,KEOS)-X)) F3=FERMI(A(13,KEOS)*(A(14,KEOS)-X)) F4=FERMI(A(17,KEOS)*(A(18,KEOS)-X)) PLG=P01*F1+P02*F2+P03*F3+P04*F4 DZDX=F1/D4* ! d\zeta/dX=d ln P/d ln rho * ((A(2,KEOS)-A(1,KEOS)*A(4,KEOS)+3.*A(3,KEOS)*X**2+ + 2.*A(3,KEOS)*A(4,KEOS)*X**3)/D4- - A(5,KEOS)*(A(1,KEOS)+A(2,KEOS)*X+A(3,KEOS)*X**3)* * FERMI(A(5,KEOS)*(A(6,KEOS)-X))) + + F2*(A(8,KEOS)+A(9,KEOS)*FERMI(A(9,KEOS)*(X-A(10,KEOS)))*P02) + +F3*(A(12,KEOS)+A(13,KEOS)*FERMI(A(13,KEOS)*(X-A(14,KEOS)))*P03)+ +F4*(A(16,KEOS)+A(17,KEOS)*FERMI(A(17,KEOS)*(X-A(18,KEOS)))*P04) Gamma=(1.d0+10.d0**(PLG-RLG)/c2)*DZDX return end subroutine FitHEOS(KEOS,H1,RLG) * Version 04.03.04 * Fit of NS EOS for rotating configurations * Input: KEOS; H1=H/H0-1, H0=H(RHO=7.86 g/cc) \approx 1 * Output: RLG = lg(rho[g/cc]) * KEOS=1: SLy4, rho=1.e4-4.e15 g/cc, * average error 1.3%, max.error 4.0% at H1=0.0095 (rho=3.69e11 g/cc) * KEOS=2: FPS, rho=1.e4-1.e17 g/cc, * average error 0.95%, max.error 4.2% at H1=0.0225 (rho=1.23e14 g/cc) * KEOS=3: APR(+SLy), rho=1.e4-1.e17 g/cc, * average error 0.86%, max.error 2.3% at H1=0.0098 (rho=5.77e11 g/cc) implicit double precision (A-H), double precision (O-Z) save dimension A(16,3) data A/5.926,.4704,20.13,.2347,3.07,97.8,-2.012,89.85, * 34.96,15.328,.621,63.1,68.5,2.518,2.6,1.363, * 5.926,.4704,19.92,.2333,2.63,54.7,-1.926,36.89, * 11.97,15.432,.6731,49.40,11.47,1.425,3.0,.913, * 5.926,.4704,19.927,.2332,2.90,90.2,-2.002,37.254, * 11.95,15.259,.5080,31.18,11.35,.678,3.0,.762/ if (KEOS.lt.1.and.KEOS.gt.3) stop'FitEOS: KEOS' if (H1.lt.0.) stop'FitHEOS: negative H1' if (H1.gt.5.) stop'FitHEOS: too high H1' X=dlog10(H1) RLG1=A(1,KEOS)+A(2,KEOS)*X+A(3,KEOS)*H1**A(4,KEOS)/ / (1.+A(5,KEOS)*H1) RLG3=(A(10,KEOS)+A(11,KEOS)*X) G=(A(12,KEOS)*H1)**7 RLG2=(A(8,KEOS)+A(9,KEOS)*X+G*RLG3)/ / (1.+A(13,KEOS)*H1+G) C1=A(6,KEOS) C2=A(7,KEOS) RLG=RLG1*FERMI(C1*(X-C2))+RLG2*FERMI(C1*(C2-X))+ + FERMI(A(15,KEOS)*(A(16,KEOS)-X))*A(14,KEOS) 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