* Unified smooth fit to the EOS of liquid/plasma/superionic water * Reference: S. Mazevet, A. Licari, G. Chabrier, & A. Y. Potekhin, * Equation of state of dense water for planetary and exoplanetary modeling * Astronomy and Astrophysics, 2018 (submitted) * URL http://www.ioffe.ru/astro/H2O/ * Contact: Alexander Potekhin ** L I S T O F S U B R O U T I N E S : * MAIN (normally commented-out) - example driving routine * H2OFIT - implementation of the analytic fit to the equation of state * ELECNR - equation of state of the ideal electron Fermi gas (fit) * FERINT - analytic approximations to the Fermi integrals * FINVER - analytic approximations to the inverse Fermi integrals ************************************************************************ * MAIN program: Version 26.08.18 * 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 write(*,'('' Input T (K): ''$)') C%C read*,T C%C write(*,112) C%C 10 continue C%C write(*,'('' Input RHO [ 0 to stop]: ''$)') C%C read*,RHO C%C if (RHO.le.0.) stop C%C call H2OFIT(RHO,T,PnkT,FNkT,UNkT,CV,CHIT,CHIR,PMbar,USPEC) C%C SNk=UNkT-FNkT C%C write(*,111) RHO,T,PMbar*100,PnkT,FNkT,UNkT,CV,CHIT,CHIR,SNk C%C goto 10 C%C stop C%C 112 format(' rho[g/cc] T [K] P [GPa] P/(n_i kT)', C%C * ' F/(N_i kT) U/(N_i kT) C_V/(N_i k) chi_T chi_r', C%C * ' S/(N_i k)') C%C 111 format(1P,4E12.4,11E11.3) C%C end subroutine H2OFIT(RHO,T,PnkT,FNkT,UNkT,CV,CHIT,CHIR,PMbar,USPEC) * Version 06.07.15 * Fit of free energy to reproduce pressure and internal energy of H2O * (Wagner & Pruss + French et al. + Licari) * Input: RHO - mass density [g/cc] * T - temperature [K] * Output: PnkT - pressure / n_i kT, where n_i=2*n_H+n_O=3*n_{H2O} * FNkT - free energy / N_i kT (up to an additive constant) * UNkT - internal energy / N_i kT * CV - heat capacity /kN_i * CHIT - logarithmic derivative of pressure over temperature * CHIR - logarithmic derivative of pressure over density * PMbar - total pressure (Mbar) * USPEC - specific internal energy [erg/g] implicit double precision (A-H), double precision (O-Z) save parameter(PI=3.141592653d0,UN_T6=.3157746,C13=1.d0/3.d0) parameter(AUM=1822.88848d0) ! a.m.u./m_e parameter(Zmean=10.d0/3.d0,CMImean=18.d0/3.d0) parameter(TwoPI=2.d0*PI) parameter(DENSCONV=11.20587*CMImean, ! converts RHO --> n_i [au] * TnkCONV=8.31447d13/CMImean) ! TnkCONV*RHO*T6=n_i kT [erg/cc] parameter(aW=2.357,bW=340.8) !tabular van der Waals constants [au] parameter(TCRIT=.00205) ! critical temperature 647.15 K in Ha parameter(P1=2.35,P3=5.9,P4=3.78,P5=17.,P7=1.5,P8=.09, * QW=.00123797,PW=2.384,PQ=1.5d0,Q4=4.,Q1=.4,Q2=90., * PC1=.0069,PC2=.0031,PC3=.00558,PC4=.019) parameter(SBASE=4.9d0) T6=T/1.d6 TEMP=T6/UN_T6 ! T [au] DENSI=RHO/DENSCONV ! number density of all atomic nuclei [au] DENSMOL=DENSI/3.d0 ! number density of "molecules" (H+H+O) [au] PRESSI=DENSI*TEMP ! ideal-ions total pressure (normalization) RS=(.75d0/PI/DENSI/Zmean)**C13 ! r_s - electron density parameter GAME=1.d0/RS/TEMP ! electron Coulomb parameter Gamma_e GAMEsq=dsqrt(GAME) * 1. Superionic/plasma part of the EOS: ZNA=1.+P8/RS/GAMEsq ZNA1RS=-P8/RS/GAMEsq ! d ZNA / d ln RS ZNA1G=.5d0*ZNA1RS ! d ZNA / d ln GAME ZNA2RS=-ZNA1RS ! d ZNA1RS / d ln RS ZNA2RSG=.5d0*ZNA2RS ! d^2 ZNA / d ln RS d ln GAME ZNA2G=.5d0*ZNA2RSG ! d^2 ZNA / d ln GAME^2 ZNB=P1*RS/ZNA ZNB1RS=ZNB*(1.d0-ZNA1RS/ZNA) ! d ZNB / d ln RS ZNB1G=-ZNB*ZNA1G/ZNA ! d ZNB / d ln GAME ZNB2RS=ZNB1RS*(1.d0-ZNA1RS/ZNA)-ZNB*ZNA2RS/ZNA+ZNB*(ZNA1RS/ZNA)**2 ZNB2G=-ZNB1G*ZNA1G/ZNA-ZNB*ZNA2G/ZNA+ZNB*(ZNA1G/ZNA)**2 ZNB2RSG=-ZNB1RS*ZNA1G/ZNA-ZNB*ZNA2RSG/ZNA+ZNB*ZNA1G*ZNA1RS/ZNA**2 ZNC=1.d0+P5/GAME RS4=RS**P4 ZNC4=ZNC*dsqrt(ZNC) ! ZNC**P7 ZNE=P3*RS4/ZNC4 ZNE1RS=P4*ZNE ! d ZNE /d ln RS ZNE1G=ZNE*P7/ZNC*P5/GAME ! d ZNE /d ln GAME ZNE2RS=P4*ZNE1RS ZNE2G=ZNE1G*(ZNE1G/ZNE+P5/GAME/ZNC-1.d0) ZNE2RSG=P4*ZNE1G ZN=1.d0+ZNB+ZNE ZN1RS=ZNB1RS+ZNE1RS ! DZN/DlnRS ZN1G=ZNB1G+ZNE1G ! DZN/DlnGAME ZN2RS=ZNB2RS+ZNE2RS ! DZN/DlnRS ZN2G=ZNB2G+ZNE2G ! DZN/DlnGAME ZN2RSG=ZNB2RSG+ZNE2RSG ! DZN/DlnRS ZN1R=(ZN1G-ZN1RS)/3.d0 ! D ZN/Dln\rho ZN1T=-ZN1G ZN2R=(ZN2RS-2.d0*ZN2RSG+ZN2G)/9.d0 ! D2ZN/Dln\rho^2 ZN2T=ZN2G ZN2RT=-(ZN2G-ZN2RSG)/3.d0 ! D2ZN/Dln\rho DlnGAME ZEF=Zmean/ZN ! eff. ZDR=-ZN1R/ZN ! d ln ZEF / d ln rho ZDT=-ZN1T/ZN ! d ln ZEF / d ln\Gamma_e = - d ln ZEF / d ln T ZDRR=-ZN2R/ZN+(ZN1R/ZN)**2 ! d2 ln ZEF / d ln rho^2 ZDTT=-ZN2T/ZN+(ZN1T/ZN)**2 ! d2 ln ZEF / d ln T^2 ZDRT=ZN1R*ZN1T/ZN**2-ZN2RT/ZN ! d2 ln ZEF / d ln\rho d ln T DENSEF=DENSI*ZEF ! eff. n_e [a.u.] call ELECNR(DENSEF,TEMP,CHI, * FE,PE,UE,SE,CVE,CHITE,CHIRE) FNkTsi=FE*ZEF FEDR=PE*(1.d0+ZDR) FEDT=-UE+PE*ZDT FEDRR=PE*(CHIRE-1.d0)*(1.+ZDR)**2+PE*ZDRR FEDRT=(PE*(CHITE-1.d0)+PE*(CHIRE-1.d0)*ZDT)*(1.+ZDR)+PE*ZDRT FEDTT=UE-CVE+PE*(CHITE-1.d0)*ZDT+ + PE*(CHITE-1.d0)*ZDT+PE*(CHIRE-1.d0)*ZDT**2+PE*ZDTT FDR=FEDR*ZEF+FE*ZEF*ZDR FDT=FEDT*ZEF+FE*ZEF*ZDT FsiDRR=FEDRR*ZEF+FEDR*ZEF*ZDR+FEDR*ZEF*ZDR+FE*ZEF*ZDR**2+ + FE*ZEF*ZDRR FsiDTT=FEDTT*ZEF+2*FEDT*ZEF*ZDT+FE*ZEF*ZDT**2+FE*ZEF*ZDTT FsiDRT=FEDRT*ZEF+FEDR*ZEF*ZDT+FEDT*ZEF*ZDR+FE*ZEF*ZDR*ZDT+ + FE*ZEF*ZDRT PnkTsi=FDR UNkTsi=-FDT * 2. Nonideal molecular part of the EOS: cW=1.d0+(QW/TEMP)**PW cW1T=-PW*(QW/TEMP)**PW cW2T=-PW*cW1T bWPQ=bW*DENSMOL*dsqrt(bW*DENSMOL) ! (bW*DENSMOL)**PQ FNkTmol=(-aW*DENSMOL/TEMP+bW*DENSMOL+bWPQ*cW/PQ)/3.d0 PnkTmol=(-aW*DENSMOL/TEMP+bW*DENSMOL+bWPQ*cW)/3.d0 UNkTmol=-(aW*DENSMOL/TEMP+bWPQ*cW1T/PQ)/3.d0 FmDRR=(-aW*DENSMOL/TEMP+bW*DENSMOL+bWPQ*cW*PQ)/3.d0 FmDTT=-(aW*DENSMOL/TEMP-bWPQ*cW2T/PQ)/3.d0 FmDRT=(aW*DENSMOL/TEMP+bWPQ*cW1T)/3.d0 * 3. Nonideal part = combination of the superionic and molecular parts: X=Q4*dlog(Q1*RHO+Q2*TEMP) X1R=Q4*Q1*RHO/(Q1*RHO+Q2*TEMP) ! dx / d ln(RHO) X1T=Q4*Q2*TEMP/(Q1*RHO+Q2*TEMP) ! dx / d ln(T) X2R=Q4*Q1*Q2*RHO*TEMP/(Q1*RHO+Q2*TEMP)**2 ! d2x / d ln(RHO) ^2 X2T=X2R X2RT=-X2R YL=FERMIF(X) YH=1.d0-YL YH1X=YH*YL ! d YH / dx YH2X=YH1X*(YL-YH) ! d2 YH / dx2 YH1R=YH1X*X1R ! d YH / d ln(RHO) YH1T=YH1X*X1T ! d YH / d ln(T) YL1R=-YH1R ! d YL / d ln(RHO) YL1T=-YH1T ! d YL / d ln(T) YH2X=YH1X*(YL-YH) YL2X=-YH2X YH2R=YH2X*X1R**2+YH1X*X2R YH2T=YH2X*X1T**2+YH1X*X2T YH2RT=YH2X*X1R*X1T+YH1X*X2RT FNkTni=FNkTmol*YL+FNkTsi*YH PnkTni=PnkTmol*YL+PnkTsi*YH+(FNkTsi-FNkTmol)*YH1R UNkTni=UNkTmol*YL+UNkTsi*YH-(FNkTsi-FNkTmol)*YH1T FDRR=FmDRR*YL+FsiDRR*YH+2.*(PnkTsi-PnkTmol)*YH1R+ + (FNkTsi-FNkTmol)*YH2R FDTT=FmDTT*YL+FsiDTT*YH-2.*(UNkTsi-UNkTmol)*YH1T+ + (FNkTsi-FNkTmol)*YH2T FDRT=FmDRT*YL+FsiDRT*YH+ + (PnkTsi-PnkTmol)*YH1T-(UNkTsi-UNkTmol)*YH1R+ + (FNkTsi-FNkTmol)*YH2RT * 4. Ideal part of the EOS: THLmol=dsqrt(TwoPI/(18.d0*AUM*TEMP)) ! Thermal de Broglie wavelen. FNkTid=(dlog(DENSMOL*THLmol**3)-1.d0)/3.d0 ! id.-gas PnkTid=C13 UNkTid=.5d0 * 5. The sum of the nonideal and ideal EOSs: FNkT=FNkTni+FNkTid PnkT=PnkTni+PnkTid UNkT=UNkTni+UNkTid CV=UNkT-FDTT CHIR=FDRR/PnkT+1.d0 CHIT=FDRT/PnkT+1.d0 * 6. Thermal corrections: TTC=TEMP/TCRIT TL2=PC4*TTC ! PL2*TEMP ULB=TL2**2*dsqrt(TL2) ! (PL2*TEMP)**PL3, PL3=2.5 ULB1=1.d0+ULB FL=dlog(ULB1/ULB) UL=2.5d0/ULB1 ! - d FL3 / d ln T CVL=UL*(1.d0-1.5d0*ULB)/ULB1 ! -1.5d0=1.d0-PL3 * 7. Second (low-T) thermal correction: TTC2=TTC**2 ULC1=1.d0+TTC2 FC=-(PC1*dlog(ULC1/TTC2)+PC2*datan(TTC))/TCRIT-PC3/TEMP UC=((2.d0*PC1*TTC-PC2*TTC2)/ULC1-PC3)/TEMP CVC=2.d0/TCRIT*(PC1*(1.d0-TTC2)-PC2*TTC)/ULC1**2 * Total: FNkT=FNkT+FL+FC-SBASE ! the last constant fits entropy of Soubiran UNkT=UNkT+UL+UC CV=CV+CVL+CVC * Auxiliary: Tnk=TnkCONV*RHO*T6 ! n_i kT [erg/cc] PMbar=PnkT*Tnk/1.d12 ! P [Mbar] USPEC=UNkT*Tnk/RHO ! U [erg/g] return end subroutine ELECNR(DENSE,TEMP, * CHI,FEid,PEid,UEid,SEid,CVE,CHITE,CHIRE) * Version 16.02.14 * NON-RELATIVISTIC version of the ideal electron-gas EOS. * Input: DENSE - electron number density n_e [a.u.], TEMP - T [a.u.] * Output: CHI=\mu/kT, * FEid - free energy / N_e kT, UEid - internal energy / N_e kT, * PEid - pressure (P_e) / n_e kT, SEid - entropy / N_e k, * CVE - heat capacity / N_e k, * CHITE=(d ln P_e/d ln T)_V, CHIRE=(d ln P_e/d ln n_e)_T implicit double precision (A-H), double precision (O-Z) save parameter (PI=3.14159265d0) parameter (TwoPI=2.d0*PI) data KRUN/0/ if (KRUN.eq.0) then KRUN=1 SQPI=dsqrt(PI) endif CLE=dsqrt(TwoPI/TEMP) ! \loambda_e [a.u.] TPI=4.d0/(SQPI*CLE**3) ! prefactor FDENS=SQPI*CLE**3*DENSE/4.d0 ! \sqrt{\pi}\lambda_e^3 n_e/4 call FINVER(FDENS,1,CHI,XDF,XDFF) call FERINT(CHI,1,F12) call FERINT(CHI,2,F32) UEid=F32/FDENS PEid=UEid/1.5d0 FEid=CHI-PEid SEid=UEid-FEid XDF32=1.5d0*FDENS*XDF CHITE=2.5d0-XDF32/PEid CHIRE=XDF32*F12/F32 CVE=1.5d0*PEid*CHITE return end function FERMIF(X) ! Fermi function 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 FERMIF=F return end subroutine FERINT(X,N,F) ! Fermi integrals * F_q(x) : H.M.Antia 93 ApJS 84 101; relative error 0.01% * q=N-1/2=-1/2,1/2,3/2,5/2 (N=0,1,2,3) Version 16.02.14 * Stems from function FERINT v.07.03.95 implicit double precision (A-H), double precision (O-Z) save dimension A(0:7,0:3),B(0:7,0:3),C(0:11,0:3),D(0:11,0:3), * LA(0:3),LB(0:3),LC(0:3),LD(0:3) data A/1.71446374704454d7,3.88148302324068d7,3.16743385304962d7, * 1.14587609192151d7,1.83696370756153d6,1.14980998186874d5, * 1.98276889924768d3,1.d0, ! F_{-1/2} * 5.75834152995465d6,1.30964880355883d7,1.07608632249013d7, * 3.93536421893014d6,6.42493233715640d5,4.16031909245777d4, * 7.77238678539648d2,1.d0, ! F_{1/2} * 4.32326386604283d4,8.55472308218786d4,5.95275291210962d4, * 1.77294861572005d4,2.21876607796460d3,9.90562948053293d1, * 1.d0,0., ! F_{3/2} * 6.61606300631656d4,1.20132462801652d5,7.6725995316812d4, * 2.10427138842443d4,2.44325236813275d3,1.02589947781696d2, * 1.d0,0./, ! F_{5/2} * B/9.67282587452899d6,2.87386436731785d7,3.26070130734158d7, * 1.77657027846367d7,4.81648022267831d6,6.13709569333207d5, * 3.13595854332114d4,4.35061725080755d2, ! F_{-1/2} * 6.49759261942269d6,1.70750501625775d7,1.69288134856160d7, * 7.95192647756086d6,1.83167424554505d6,1.95155948326832d5, * 8.17922106644547d3,9.02129136642157d1, ! F_{1/2} * 3.25218725353467d4,7.01022511904373d4,5.50859144223638d4, * 1.95942074576400d4,3.20803912586318d3,2.20853967067789d2, * 5.05580641737527d0,1.99507945223266d-2, ! F_{3/2} * 1.99078071053871d4,3.79076097261066d4,2.60117136841197d4, * 7.97584657659364d3,1.10886130159658d3,6.35483623268093d1, * 1.16951072617142d0,3.31482978240026d-3/, ! F_{5/2} * C/-4.46620341924942d-15,-1.58654991146236d-12, * -4.44467627042232d-10,-6.84738791621745d-08, * -6.64932238528105d-06,-3.69976170193942d-04, * -1.12295393687006d-02,-1.60926102124442d-01, * -8.52408612877447d-01,-7.45519953763928d-01, * 2.98435207466372d+00,1.d0, ! F_{-1/2} * 4.85378381173415d-14,1.64429113030738d-11,3.76794942277806d-9, * 4.69233883900644d-7,3.40679845803144d-5,1.32212995937796d-3, * 2.60768398973913d-2,2.48653216266227d-1,1.08037861921488d0, * 1.91247528779676d+0,1.d0,0.d0, ! F_{1/2} * 2.80452693148553d-13,8.60096863656367d-11,1.62974620742993d-8, * 1.63598843752050d-6,9.12915407846722d-5,2.62988766922117d-3, * 3.85682997219346d-2,2.78383256609605d-1,9.02250179334496d-1, * 1.d0,0.d0,0.d0, ! F_{3/2} * 8.42667076131315d-12,2.31618876821567d-09,3.54323824923987d-7, * 2.77981736000034d-5,1.14008027400645d-3,2.32779790773633d-2, * 2.39564845938301d-1,1.24415366126179d00,3.18831203950106d0, * 3.42040216997894d0,1.d0,0.d0/, ! F_{5/2} * D/-2.23310170962369d-15,-7.94193282071464d-13, * -2.22564376956228d-10,-3.43299431079845d-08, * -3.33919612678907d-06,-1.86432212187088d-04, * -5.69764436880529d-03,-8.34904593067194d-02, * -4.78770844009440d-01,-4.99759250374148d-01, * 1.86795964993052d+00, 4.16485970495288d-01, ! F_{-1/2} * 7.28067571760518d-14,2.45745452167585d-11,5.62152894375277d-9, * 6.96888634549649d-7,5.02360015186394d-5,1.92040136756592d-3, * 3.66887808001874d-2,3.24095226486468d-1,1.16434871200131d0, * 1.34981244060549d0,2.01311836975930d-1,-2.14562434782759d-2, * 7.01131732871184d-13,2.10699282897576d-10,3.94452010378723d-8, * 3.84703231868724d-6,2.04569943213216d-4,5.31999109566385d-3, * 6.39899717779153d-2,3.14236143831882d-1,4.70252591891375d-1, * -2.15540156936373d-2,2.34829436438087d-3,0.d0, ! F_{3/2} * 2.94933476646033d-11,7.68215783076936d-09,1.12919616415947d-6, * 8.09451165406274d-5,2.81111224925648d-3,3.99937801931919d-2, * 2.27132567866839d-1,5.31886045222680d-1,3.70866321410385d-1, * 2.27326643192516d-2,0.d0,0.d0/, ! F_{5/2} * LA/7,7,6,6/,LB/7,7,7,7/,LC/11,10,9,1/,LD/11,11,10,9/ if (N.lt.0.or.N.gt.3) stop'FERINT: Invalid subscript' if (X.lt.2.d0) then T=dexp(X) UP=0. DOWN=0. do I=LA(N),0,-1 UP=UP*T+A(I,N) enddo do I=LB(N),0,-1 DOWN=DOWN*T+B(I,N) enddo F=T*UP/DOWN else T=1.d0/X**2 UP=0. DOWN=0. do I=LC(N),0,-1 UP=UP*T+C(I,N) enddo do I=LD(N),0,-1 DOWN=DOWN*T+D(I,N) enddo F=dsqrt(X)*X**N*UP/DOWN endif return end subroutine FINVER(F,N,X,XDF,XDFF) ! Inverse of Fermi integrals * Version 24.05.07 * X_q(f)=F^{-1}_q(f) : H.M.Antia 93 ApJS 84, 101 * q=N-1/2=-1/2,1/2,3/2,5/2 (N=0,1,2,3) * Input: F - argument, N=q+1/2 * Output: X=X_q, XDF=dX/df, XDFF=d^2 X / df^2 * Relative error: N = 0 1 2 3 * for X: 3.e-9, 4.2e-9, 2.3e-9, 6.2e-9 * jump at f=4: * for XDF: 6.e-7, 5.4e-7, 9.6e-8, 3.1e-7 * for XDFF: 4.7e-5, 4.8e-5, 2.3e-6, 1.5e-6 implicit double precision (A-H), double precision (O-Z) save dimension A(0:5,0:3),B(0:6,0:3),C(0:6,0:3),D(0:6,0:3), * LA(0:3),LB(0:3),LD(0:3) data A/-1.570044577033d4,1.001958278442d4,-2.805343454951d3, * 4.121170498099d2,-3.174780572961d1,1.d0, ! X_{-1/2} * 1.999266880833d4,5.702479099336d3,6.610132843877d2, * 3.818838129486d1,1.d0,0., ! X_{1/2} * 1.715627994191d2,1.125926232897d2,2.056296753055d1,1.d0,0.,0., * 2.138969250409d2,3.539903493971d1,1.d0,0.,0.,0./, ! X_{5/2} * B/-2.782831558471d4,2.886114034012d4,-1.274243093149d4, * 3.063252215963d3,-4.225615045074d2,3.168918168284d1, * -1.008561571363d0, ! X_{-1/2} * 1.771804140488d4,-2.014785161019d3,9.130355392717d1, * -1.670718177489d0,0.,0.,0., ! X_{1/2} * 2.280653583157d2,1.193456203021d2,1.16774311354d1, * -3.226808804038d-1,3.519268762788d-3,0.,0., ! X_{3/2} * 7.10854551271d2,9.873746988121d1,1.067755522895d0, * -1.182798726503d-2,0.,0.,0./, ! X_{5/2} * C/2.206779160034d-8,-1.437701234283d-6,6.103116850636d-5, * -1.169411057416d-3,1.814141021608d-2,-9.588603457639d-2,1.d0, * -1.277060388085d-2,7.187946804945d-2,-4.262314235106d-1, * 4.997559426872d-1,-1.285579118012d0,-3.930805454272d-1,1.d0, * -6.321828169799d-3,-2.183147266896d-2,-1.05756279932d-1, * -4.657944387545d-1,-5.951932864088d-1,3.6844711771d-1,1.d0, * -3.312041011227d-2,1.315763372315d-1,-4.820942898296d-1, * 5.099038074944d-1,5.49561349863d-1,-1.498867562255d0,1.d0/, * D/8.827116613576d-8,-5.750804196059d-6,2.429627688357d-4, * -4.601959491394d-3,6.932122275919d-2,-3.217372489776d-1, * 3.124344749296d0, ! X_{-1/2} * -9.745794806288d-3,5.485432756838d-2,-3.29946624326d-1, * 4.077841975923d-1,-1.145531476975d0,-6.067091689181d-2,0., * -4.381942605018d-3,-1.5132365041d-2,-7.850001283886d-2, * -3.407561772612d-1,-5.074812565486d-1,-1.387107009074d-1,0., * -2.315515517515d-2,9.198776585252d-2,-3.835879295548d-1, * 5.415026856351d-1,-3.847241692193d-1,3.739781456585d-2, * -3.008504449098d-2/, ! X_{5/2} * LA/5,4,3,2/,LB/6,3,4,3/,LD/6,5,5,6/ if (N.lt.0.or.N.gt.3) stop'FINVER: Invalid subscript' if (F.le.0.) stop'FINVER: Non-positive argument' if (F.lt.4.) then T=F UP=0. UP1=0. UP2=0. DOWN=0. DOWN1=0. DOWN2=0. do I=LA(N),0,-1 UP=UP*T+A(I,N) if (I.ge.1) UP1=UP1*T+A(I,N)*I if (I.ge.2) UP2=UP2*T+A(I,N)*I*(I-1) enddo do I=LB(N),0,-1 DOWN=DOWN*T+B(I,N) if (I.ge.1) DOWN1=DOWN1*T+B(I,N)*I if (I.ge.2) DOWN2=DOWN2*T+B(I,N)*I*(I-1) enddo X=dlog(T*UP/DOWN) XDF=1.d0/T+UP1/UP-DOWN1/DOWN XDFF=-1.d0/T**2+UP2/UP-(UP1/UP)**2-DOWN2/DOWN+(DOWN1/DOWN)**2 else P=-1./(.5+N) ! = -1/(1+\nu) = power index T=F**P ! t - argument of the rational fraction T1=P*T/F ! dt/df T2=P*(P-1.)*T/F**2 ! d^2 t / df^2 UP=0. UP1=0. UP2=0. DOWN=0. DOWN1=0. DOWN2=0. do I=6,0,-1 UP=UP*T+C(I,N) if (I.ge.1) UP1=UP1*T+C(I,N)*I if (I.ge.2) UP2=UP2*T+C(I,N)*I*(I-1) enddo do I=LD(N),0,-1 DOWN=DOWN*T+D(I,N) if (I.ge.1) DOWN1=DOWN1*T+D(I,N)*I if (I.ge.2) DOWN2=DOWN2*T+D(I,N)*I*(I-1) enddo R=UP/DOWN R1=(UP1-UP*DOWN1/DOWN)/DOWN ! dR/dt R2=(UP2-(2.*UP1*DOWN1+UP*DOWN2)/DOWN+2.*UP*(DOWN1/DOWN)**2)/ / DOWN X=R/T RT=(R1-R/T)/T XDF=T1*RT XDFF=T2*RT+T1**2*(R2-2.*RT)/T endif return end