* interpolation of the thermodynamic functions and Rosseland mean * opacities from tables "hmag*.dat" * Last change: 19.06.03 * Insignificant modifications: "pause" changed ro "write" 07.11.09 * output to a file is added to MAIN 25.03.10 * Output is in the same format as the tables being interpolated; * see the description in the file "hmagnet.txt" * or at http://www.ioffe.ru/astro/NSG/Hmagnet/ * REFERENCES: A.Y.Potekhin & G.Chabrier, 2003, Astrophys. J. 585, 955; * * A.Y.Potekhin & G.Chabrier, 2004, Astrophys. J. 600, 317 * * Please send your questions or comments to Alexander Potekhin, * E-mail: palex@astro.ioffe.ru * This MAIN program block (may be commented-out) serves only to * demonstrate the use of HMINTER, * the basic interpolation subroutine placed below. C%C implicit double precision (A-H), double precision (O-Z) C%C character KEY C%C parameter(BLGmin=11.9d0,BLGmax=15.d0, C%C * TLGmin=5.3d0,TLGmax=7.6d0,RLGmin=-7.4d0,RLGmax=3.6d0) C%C dimension AOUT(13) C%C data AOUT/13*0./ C%C open(2,file='hminter.out') C%C write(*,110) C%C write(2,110) C%C 1 write(*,'('' INPUT: lg B(G): ''$)') C%C read*,BLG C%C if (BLG.lt.BLGmin.or.BLG.gt.BLGmax) then C%C print*,'choose Lg B between',BLGmin,' and',BLGmax C%C goto 1 C%C endif C%C 2 write(*,'('' Lg T(K): ''$)') C%C read*,TLG C%C if (TLG.lt.TLGmin.or.TLG.gt.TLGmax) then C%C print*,'choose Lg T between',TLGmin,' and',TLGmax C%C goto 2 C%C endif C%C RHLmin=(TLG-6.)*3.+RLGmin C%C RHLmax=(TLG-6.)*3.+RLGmax C%C write(*,'(2F10.4)') TLG,BLG C%C write(2,'(2F10.4)') TLG,BLG C%C 3 write(*,'('' Lg rho (g/cm^3): ''$)') C%C read*,RHL C%C if (RHL.lt.RHLmin.or.RHL.gt.RHLmax) then C%C print*,'choose Lg rho between',RHLmin,' and',RHLmax C%C goto 3 C%C endif C%C RLG=RHL-(TLG-6.)*3. C%C do Kout=1,13 C%C call HMINTER(BLG,TLG,RLG,Kout,OUTPUT,Kret) C%C AOUT(Kout)=OUTPUT C%C enddo C%C* ---------------------- OUTPUT: ----------------------------- * C%C write(*,114) RLG,(AOUT(I),I=1,13) C%C write(2,114) RLG,(AOUT(I),I=1,13) C%C write(*,'('' New density? (Y/N) ''$)') C%C read(*,'(A)') KEY C%C if (KEY.ne.'n'.and.KEY.ne.'N') goto 3 C%C write(*,'('' New temperature? (Y/N) ''$)') C%C read(*,'(A)') KEY C%C if (KEY.ne.'n'.and.KEY.ne.'N') goto 2 C%C write(*,'('' New B? (Y/N) ''$)') C%C read(*,'(A)') KEY C%C if (KEY.ne.'n'.and.KEY.ne.'N') goto 1 C%C stop C%C 110 format(' INTERPOLATED EOS AND OPACITIES.'/ C%C * ' Format of the output:'/ C%C * 22X,'thermodynamic functions',14X,'|',8X, C%C * 'number fractions',12X,'lg(opacities)'/ C%C * ' lg(R) lg(P) PV/(NkT) U/(NkT) S/(Nk) Cv/(Nk) chit chir ', C%C * ' x(H) x(H0) x(H2) x(pert.) long. transv.'/ C%C * ' lg(T) lg(B)'/ C%C * ' (1) (2) (3) (4) (5) (6) (7) (8)', C%C * ' (9) (10) (11) (12) (13) (14)') C%C 114 format(F6.2,F9.4,2F8.3,4F7.3,1P,4E9.2,0P,3F7.3) C%C end ** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * This subroutine interpolates the thermodynamic functions * * and Rosseland mean opacities of magnetized hydrogen * * from data files of type "hmag1?_?.dat" or "hmn1?_?.dat * * REFERENCES: A.Y.Potekhin & G.Chabrier, 2003, Astrophys. J. 585, 955; * * A.Y.Potekhin & G.Chabrier, 2004, Astrophys. J. 600, 317 * * Web: http://www.ioffe.ru/astro/NSG/Hmagnet/ * * E-mail: palex@astro.ioffe.ru (Alexander Potekhin) * * ---------------------------------------------------- Version 24.05.02 subroutine HMINTER(BLG1,TLG1,RLG1,Kout,OUTPUT,Kret) * Input: BLG - lg(B[G}), TLG - lg(T[K]), RLG - lg(R), where R=rho/T_6^3 * Output is regulated by the key Kout as follows: * Kout OUTPUT * 1 PLG=lg(P[bar]), where P is the pressure * 2 PNkT=PV/(NkT), where N is number of protons (free and bound) * 3 UNkT=U/(NkT) * 4 ENTROPY=S/(Nk) * 5 CVREL=C_V/Nk * 6 CHIT=\chi_T = (d log P/d log T)_V * 7 CHIR=\chi_\rho = (d log P/d log rho)_T * 8 FHOPT - the "TRUE" ("optical") atomic fraction * 9 FHOPT0 - the true ground-state atomic fraction * 10 FH2OPT2 - number of H2 molecules divided by N * 11 FCLUST - the fraction of protons comprised in clusters * 12 ROSLONG=lg(K0), K0 is the effective Rosseland mean opacity * [cm**2/g] for radiative transport along B * 13 ROSTRAN=lg(K1) - the same but perpendicular to B * Kret=0 if the calculation is normal; Kret=-1 otherwise (then output=0) implicit double precision (A-H), double precision (O-Z) character*12 AFNAME,AFNAME3 save parameter (NTmax=34,NRmax=55,MAXB=17) parameter (NTmax3=19,NRmax3=54,MAXB3=16) parameter (MAXT=NTmax+1,MAXR=NRmax+1) parameter (MAXT3=NTmax3+1,MAXR3=NRmax3+1) parameter (MAXDAT=13) dimension AT(MAXT),AR(MAXR),AB(MAXB),ADAT(MAXR,MAXT,MAXB,MAXDAT), * AFNAME(MAXB) dimension AT3(MAXT3),AR3(MAXR3),AB3(MAXB3), * ADAT3(MAXR3,MAXT3,MAXB3,MAXDAT), * AFNAME3(MAXB3) data AFNAME/'hmag11_9.dat','hmag12_0.dat','hmag12_1.dat', * 'hmag12_2.dat','hmag12_3.dat','hmag12_4.dat','hmag12_5.dat', * 'hmag12_6.dat','hmag12_7.dat','hmag12_8.dat','hmag12_9.dat', * 'hmag13_0.dat','hmag13_1.dat','hmag13_2.dat','hmag13_3.dat', * 'hmag13_4.dat','hmag13_5.dat'/ data AB/11.9d0,12.0d0,12.1d0,12.2d0,12.3d0, * 12.4d0,12.5d0,12.6d0,12.7d0,12.8d0, * 12.9d0,13.0d0,13.1d0,13.2d0,13.3d0, * 13.4d0,13.5d0/ data AFNAME3/ * 'hmn13_5.dat','hmn13_6.dat','hmn13_7.dat','hmn13_8.dat', * 'hmn13_9.dat','hmn14_0.dat','hmn14_1.dat','hmn14_2.dat', * 'hmn14_3.dat','hmn14_4.dat','hmn14_5.dat','hmn14_6.dat', * 'hmn14_7.dat','hmn14_8.dat','hmn14_9.dat','hmn15_0.dat'/ data AB3/13.5d0,13.6d0,13.7d0,13.8d0, * 13.9d0,14.0d0,14.1d0,14.2d0, * 14.3d0,14.4d0,14.5d0,14.6d0, * 14.7d0,14.8d0,14.9d0,15.0d0/ if (BLG1.lt.13.5d0) then call HMINTER3(BLG1,TLG1,RLG1,Kout,OUTPUT,Kret, * MAXB,MAXT,MAXR,MAXDAT,AT,AR,AB,ADAT,AFNAME) else call HMINTER3(BLG1,TLG1,RLG1,Kout,OUTPUT,Kret, * MAXB3,MAXT3,MAXR3,MAXDAT,AT3,AR3,AB3,ADAT3,AFNAME3) endif return end subroutine HMINTER3(BLG1,TLG1,RLG1,Kout,OUTPUT,Kret, * MAXB,MAXT,MAXR,MAXDAT,AT,AR,AB,ADAT,AFNAME) * This subroutine must be used together with HMINTER. * --------------------------------------------------- Version 19.06.03 * Insignificant modification ("pause" changed ro "write") 07.11.09 implicit double precision (A-H), double precision (O-Z) character*12 FNAME,AFNAME save dimension AT(MAXT),AR(MAXR),AB(MAXB),ADAT(MAXR,MAXT,MAXB,MAXDAT), * AFNAME(MAXB) data KRUN/-1/,KWARN/-1/,KWB/-1/,KWR/-1/,KWT/-1/ Kret=0 if (Kout.lt.1.or.Kout.gt.13) stop'HMINTER3: Kout out of range' if ((BLG1.lt.13.5d0.and.KRUN.ne.2).or. . (BLG1.ge.13.5d0.and.KRUN.ne.3)) then ! Reading print*,'Reading EOS and opacity tables...' do IBLG=1,MAXB FNAME=AFNAME(IBLG) open(13,file=FNAME,status='OLD') write(*,'(I3,'': '',A12)') IBLG,FNAME do I=1,5 read(13,'(A)') ! skip lines enddo do ITLG=1,MAXT read(13,*) TLG,BLG if (BLG.ne.AB(IBLG)) stop'HMINTER3: B mismatch' if (IBLG.eq.1) then AT(ITLG)=TLG else if (AT(ITLG).ne.TLG) stop'HMINTER3: T mismatch' endif do IRLG=1,MAXR read(13,*) RLG,(ADAT(IRLG,ITLG,IBLG,IDAT),IDAT=1,MAXDAT) if (IBLG.eq.1.and.ITLG.eq.1) then AR(IRLG)=RLG else if (AR(IRLG).ne.RLG) stop'HMINTER3: R mismatch' endif do IDAT=8,11 ! take lg (number fractions) DAT=dlog10(dmax1(ADAT(IRLG,ITLG,IBLG,IDAT),1.d-21)) ADAT(IRLG,ITLG,IBLG,IDAT)=DAT enddo enddo ! next R enddo ! next T close(13) enddo ! next B if (BLG1.lt.13.5d0) KRUN=2 if (BLG1.ge.13.5d0) KRUN=3 print*,' done.' endif if (RLG1.gt.3.851+3.*(BLG1/2.-TLG1).and.KWARN.lt.0) then print*,CHAR(7) print*,'HMINTER3 Warning: B is weakly quantized at this R' write(*,'('' (lg B ='',F9.3,'', lg R ='',F9.3,2H).)') BLG1,RLG1 print*,'A danger of large errors due to oscillations exists.' print*,'It is not recommended to use interpolation ', * 'under such conditions. [ to continue]' read(*,'(A)') KWARN=1 endif call HUNT(AB,MAXB,BLG1,IB) if (IB.eq.0.or.IB.eq.MAXB) then if (KWB.lt.0) then write(*,'('' HMINTER3: B out of range; output := 0.'', * '' [ to continue]''$)') read(*,'(A)') endif KWB=1 OUTPUT=0. Kret=-1 return endif IBM=max0(1,IB-1) IBP=min0(MAXB,IB+2) call HUNT(AT,MAXT,TLG1,IT) if (IT.eq.0.or.IT.eq.MAXT) then if (KWT.lt.0) then write(*,'('' HMINTER3: T out of range; output := 0.'', * '' [ to continue]''$)') read(*,'(A)') endif KWT=1 OUTPUT=0. Kret=-1 return endif call HUNT(AR,MAXR,RLG1,IR) if (IR.eq.0.or.IR.eq.MAXR) then if (KWR.lt.0) then write(*,'('' HMINTER3: R out of range; output := 0.'', * '' [ to continue]''$)') read(*,'(A)') endif KWR=1 OUTPUT=0. Kret=-1 return endif * ITM=max0(1,IT-1) IT1=IT+1 ITP=min0(MAXT,IT+2) IRM=max0(1,IR-1) IR1=IR+1 IRP=min0(MAXR,IR+2) * BM,TM: call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,ITM,IBM,Kout),ADAT(IR,ITM,IBM,Kout), * ADAT(IR1,ITM,IBM,Kout),ADAT(IRP,ITM,IBM,Kout),BMTM) * BM,T0: if (IT.eq.ITM) then BMT0=BMTM else call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,IT,IBM,Kout),ADAT(IR,IT,IBM,Kout), * ADAT(IR1,IT,IBM,Kout),ADAT(IRP,IT,IBM,Kout),BMT0) endif * BM,T1: call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,IT1,IBM,Kout),ADAT(IR,IT1,IBM,Kout), * ADAT(IR1,IT1,IBM,Kout),ADAT(IRP,IT1,IBM,Kout),BMT1) * BM,TP: if (ITP.eq.IT1) then BMTP=BMT1 else call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,ITP,IBM,Kout),ADAT(IR,ITP,IBM,Kout), * ADAT(IR1,ITP,IBM,Kout),ADAT(IRP,ITP,IBM,Kout),BMTP) endif * BM: Cubic interpolation in TLG: call CINTERP3s(AT(ITM),AT(IT),AT(IT1),AT(ITP),TLG1,IT,MAXT, * BMTM,BMT0,BMT1,BMTP,VBM) * B0: if (IBM.eq.IB) then VB0=VBM goto 1 endif * B0,TM: call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,ITM,IB,Kout),ADAT(IR,ITM,IB,Kout), * ADAT(IR1,ITM,IB,Kout),ADAT(IRP,ITM,IB,Kout),B0TM) * B0,T0: if (IT.eq.ITM) then B0T0=B0TM else call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,IT,IB,Kout),ADAT(IR,IT,IB,Kout), * ADAT(IR1,IT,IB,Kout),ADAT(IRP,IT,IB,Kout),B0T0) endif * B0,T1: call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,IT1,IB,Kout),ADAT(IR,IT1,IB,Kout), * ADAT(IR1,IT1,IB,Kout),ADAT(IRP,IT1,IB,Kout),B0T1) * B0,TP: if (ITP.eq.IT1) then B0TP=B0T1 else call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,ITP,IB,Kout),ADAT(IR,ITP,IB,Kout), * ADAT(IR1,ITP,IB,Kout),ADAT(IRP,ITP,IB,Kout),B0TP) endif * B0: Cubic interpolation in TLG: call CINTERP3s(AT(ITM),AT(IT),AT(IT1),AT(ITP),TLG1,IT,MAXT, * B0TM,B0T0,B0T1,B0TP,VB0) 1 CONTINUE * B1,TM: IB1=IB+1 call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,ITM,IB1,Kout),ADAT(IR,ITM,IB1,Kout), * ADAT(IR1,ITM,IB1,Kout),ADAT(IRP,ITM,IB1,Kout),B1TM) * B1,T0: if (IT.eq.ITM) then B1T0=B1TM else call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,IT,IB1,Kout),ADAT(IR,IT,IB1,Kout), * ADAT(IR1,IT,IB1,Kout),ADAT(IRP,IT,IB1,Kout),B1T0) endif * B1,T1: call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,IT1,IB1,Kout),ADAT(IR,IT1,IB1,Kout), * ADAT(IR1,IT1,IB1,Kout),ADAT(IRP,IT1,IB1,Kout),B1T1) * B1,TP: if (ITP.eq.IT1) then B1TP=B1T1 else call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,ITP,IB1,Kout),ADAT(IR,ITP,IB1,Kout), * ADAT(IR1,ITP,IB1,Kout),ADAT(IRP,ITP,IB1,Kout),B1TP) endif * B1: Cubic interpolation in TLG: call CINTERP3s(AT(ITM),AT(IT),AT(IT1),AT(ITP),TLG1,IT,MAXT, * B1TM,B1T0,B1T1,B1TP,VB1) * BP: if (IBP.eq.IB1) then VBP=VB1 goto 2 endif * BP,TM: call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,ITM,IBP,Kout),ADAT(IR,ITM,IBP,Kout), * ADAT(IR1,ITM,IBP,Kout),ADAT(IRP,ITM,IBP,Kout),BPTM) * BP,T0: if (IT.eq.ITM) then BPT0=BPTM else call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,IT,IBP,Kout),ADAT(IR,IT,IBP,Kout), * ADAT(IR1,IT,IBP,Kout),ADAT(IRP,IT,IBP,Kout),BPT0) endif * BP,T1: call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,IT1,IBP,Kout),ADAT(IR,IT1,IBP,Kout), * ADAT(IR1,IT1,IBP,Kout),ADAT(IRP,IT1,IBP,Kout),BPT1) * BP,TP: if (ITP.eq.IT1) then BPTP=BPT1 else call CINTERP3s(AR(IRM),AR(IR),AR(IR1),AR(IRP), * RLG1,IR,MAXR, * ADAT(IRM,ITP,IBP,Kout),ADAT(IR,ITP,IBP,Kout), * ADAT(IR1,ITP,IBP,Kout),ADAT(IRP,ITP,IBP,Kout),BPTP) endif * BP: Cubic interpolation in TLG: call CINTERP3s(AT(ITM),AT(IT),AT(IT1),AT(ITP),TLG1,IT,MAXT, * BPTM,BPT0,BPT1,BPTP,VBP) 2 CONTINUE * Cubic interpolation in B: call CINTERP3s(AB(IBM),AB(IB),AB(IB1),AB(IBP),BLG1,IB,MAXB, * VBM,VB0,VB1,VBP,OUTPUT) if (Kout.ge.8.and.Kout.le.11) then if (OUTPUT.lt.-20.) then OUTPUT=0. else OUTPUT=10.d0**OUTPUT endif endif 50 return end subroutine CINTERP3s(ZM,Z0,Z1,ZP,Z,N0,MXNV,VM,V0,V1,VP,VF) * Given 4 values of Z and 4 values of V, find VF corresponding to 5th Z * Version 20.05.02 * This is a variant of CINTERP3 (v13.05.99) with reduced output * Output: VF - interpolated value of function implicit double precision (A-H), double precision (O-Z) save if (N0.le.0.or.N0.ge.MXNV) stop'CINTERP3s: N0 out of range' X=Z-Z0 H=Z1-Z0 ! basic interval XH=X/H if (N0.gt.1) then HM=Z0-ZM ! left adjoint interval V01=((V1-V0)/H**2+(V0-VM)/HM**2)/ / (1./H+1./HM) ! left derivative endif if (N0.lt.MXNV-1) then HP=ZP-Z1 ! right adjoint interval V11=((V1-V0)/H**2+(VP-V1)/HP**2)/ / (1./H+1./HP) ! right derivative endif if (N0.gt.1.and.N0.lt.MXNV-1) then ! Cubic interpolation C2=3.*(V1-V0)-H*(V11+2.*V01) C3=H*(V01+V11)-2.*(V1-V0) VF=V0+V01*X+C2*XH**2+C3*XH**3 goto 10 endif if (N0.eq.1) then ! Quadratic interpolation C2=V0-V1+V11*H VF=V1-V11*(H-X)+C2*(1.-XH)**2 else ! N0=MXNV-1 C2=V1-V0-V01*H VF=V0+V01*X+C2*XH**2 endif 10 return end subroutine HUNT(XX,N,X,JLO) * W.H.Press, B.P.Flannery, S.A.Teukolsky, W.T.Vetterling * Numerical Receipes(Cambridge Univ., 1986) * Given an array XX of length N, and given a value X, * returns a value JLO such that X is between XX(JLO) and XX(JLO+1). * XX must be monotonic, either increasing or decreasing. * JLO=0 or JLO=N is returned to indicate that X is out of range. * JLO on input is taken as the initial guess for JLO on output. implicit double precision (A-H), double precision (O-Z) dimension XX(*) logical ASCND ASCND=XX(N).gt.XX(1) ! true if ascending order, false otherwise if (JLO.le.0.or.JLO.gt.N) then ! Input guess not useful. JLO=0 JHI=N+1 goto 3 ! go immediately to bisection endif INC=1 ! set the hunting increment if (X.ge.XX(JLO).eqv.ASCND) then ! Hunt up: 1 JHI=JLO+INC if (JHI.gt.N) then ! Done hunting, since off end of table JHI=N+1 elseif (X.ge.XX(JHI).eqv.ASCND) then ! Not done hunting JLO=JHI INC=INC+INC goto 1 endif else ! Hunt down: JHI=JLO 2 JLO=JHI-INC if (JLO.lt.1) then ! Done hunting, since off end of table JLO=0 elseif (X.lt.XX(JLO).eqv.ASCND) then ! Not done hunting JHI=JLO INC=INC+INC ! so double the increment goto 2 ! and try again endif ! Done hunting, value bracketed endif * Hunt is done, so begin the final bisection phase: 3 if (JHI-JLO.eq.1) return JM=(JHI+JLO)/2. if (X.ge.XX(JM).eqv.ASCND) then JLO=JM else JHI=JM endif goto 3 end