C*********************************************************************** C Hydrogen bound opacity. Designed to give the total neutral H opacity C (bound-bound and bound free) at a given wavelength for a given list of C lines at specified plasma conditions. Occupation probability formalism C is used. Appropriate for applications like model stellar atmosphere C calculations. C C Also contains the code for normalised line profiles HLINOP which can C be used separately. C C Paul Barklem, Nik Piskunov, Uppsala August 2003. C C A number of parts from or based on code by Deane Peterson and Bob C Kurucz. C C Some contributions and significant input from Kjell Eriksson. C C We would very much appreciate bug reports to: barklem@astro.uu.se C C Table of Contents: C ------------------ C C REAL HBOP(WAVE,N,NLO,NUP,WAVEH,NH,NHE,NE,T,DOP,NPOP,NL,TOTAL,CONTIN) C REAL FUNCTION HLINOP(WAVE,NBLO,NBUP,WAVEH,T,XNE,H1FRC,HE1FRC,DOPPLE) C REAL FUNCTION HPROFL(N,M,WAVE,WAVEH,T,XNE,H1FRC,HE1FRC,DOPPH) C REAL FUNCTION STARK1(N,M,WAVE,WAVEH,T,XNE) C REAL FUNCTION HFNM(N,M) C REAL FUNCTION VCSE1F(X) C REAL FUNCTION SOFBET(B,P,N,M) C REAL*8 FUNCTION AIRVAC(W) C REAL*8 FUNCTION VACAIR(W) C REAL FUNCTION WCALC(NH, NE, NHW, NS, TEMP) C REAL FUNCTION HBF(WAVE,NLO) C SUBROUTINE LOGINT(X,Y,N,XI,YI) C SUBROUTINE POLINT(XA,YA,N,X,Y,DY) C C Variables to be set / aware of before use: C ------------------------------------------ C C HBOP : NLEVELS, NBF C HLINOP: CUTNEXT, REFLECT C HPROFL: RAYLCUT C C*********************************************************************** SUBROUTINE HBOP(WAVE, N, NLO, NUP, WAVEH, NH, NHE, NE, T, DOP, * NPOP, NL, TOTAL, CONTIN) C C Returns the total absorption coefficient due to H bound levels at WAVE, C for a given line list employing the occupation probability formalism C of Daeppen, Anderson & Milhalas ApJ 1987 (DAM). C C The extension of the occupation probability formalism to non-LTE by C Hubeny, Hummer and Lanz A&A 1994 (HHL) is implemented (though see C comments below). C C To use the code in LTE: C set NL = 0 C NPOP is then not used, population computed internally C To use the code in non-LTE: C set NL = number of states to treat in non-LTE C set NPOP = array of number densities for each state C states above NL will be treated in LTE if needed C C Note that Hubeny, Hummer and Lanz 1994 have extended the occupation C probability formalism to non-LTE, though in a somewhat unclear manner C in my opinion. The main point to be aware of is that w(HHL) does not C equal w(DAM). The first is a conditional probability and the second C is unconditional. The problem is that the definition in sect 2.1 of HHL C and derivation for w in HHL Appendix A is for the unconditional C probability w(DAM). In the rest of the paper it's conditional. The paper C doesn't really make this difference clear IMO. C C Thus w(j)(HHL) = P(j|i) = P(i & j)/P(i) = w(j)/w(i) for j>i C = 1 for j 3 p-states play a very small C role, and van der Waals might even dominate, there is however no C available theory for this. The lower state resonance broadening is C used as a guess. C IF (N.NE.1) THEN RESONT = HFNM(1,N)/(1.-1./ XN2) ELSE RESONT = HFNM(1,M)/(1.-1./ XM2) ENDIF RESONT = RESONT * 2.07E-24/GNM VDW = 4.45E-26/GNM*(XM2*(7.*XM2+5.))**0.4 STARK = 1.6678E-18*FREQNM*XKNM C C Approximately include fine structure. Exact pattern for alpha lines, C M infinite pattern used for others. C IF (N.GT.4 .OR. M.GT.10) THEN IFINS = 1 FINEST(1) = 0. FINSWT(1) = 1. ELSE IF (MMN.GT.1) THEN IFINS = LNCOMP(N) DO 1 I = 1,IFINS FINEST(I) = STCOMP(I, N)*1.D7 FINSWT(I) = STCPWT(I, N)/XN2 1 CONTINUE ELSE C C eg: Ly alpha IFINS=2, IPOS=1, FINEST=-7.3E9,3.7E9, FINSWT=1/3, 2/3 C IFINS = LNGHAL(N) IPOS = ISTAL(N) DO 2 I = 1,IFINS K = IPOS-1+I FINEST(I) = STALPH(K)*1.D7 FINSWT(I) = STWTAL(K)/XN2/3.D0 2 CONTINUE END IF END IF C C Now compute the profile for the given physical conditions. C C Firstly half-widths: DOPPH, HWSTK, HWLOR, are dnu/nu. DOP is the C doppler width dnu. C HWSTK = STARK*FO DOP = DOPPH*FREQNM C C For low Balmer lines use Barklem, Piskunov & O'Mara (2000, A&A 363, C 1091) instead of Ali & Griem for self-broadening. 1E6 factor is due C to SI -> CGS. C IF ((N.EQ.2) .AND. (M.LE.5)) THEN VBARH = SQRT(42008.*T) RESONT = BALFAC(M-N) * 1.E4*SIGMA(M-N) * * (VBARH/1.E4)**(1.-ALF(M-N)) RESONT = RESONT/FREQNM*1.E6/2./PI ENDIF C HWLOR = RESONT*H1FRC + VDW*T3NHE + RADAMP C C *Approximate* convolution is achieved by two cases: C (a) close to the core -> take the dominant mechanism. C (b) outside the core -> add all mechanisms. C C Mechanisms are divided into three: C (1) Doppler core with fine structure. C (2) (Lorentz) damping parts - self, natural and helium. C (3) Stark broadening. C C C Determine the largest half-width C NWID=1, Doppler C =2, Lorentz C =3, Stark C IF (DOPPH.GE.HWLOR.AND.DOPPH.GE.HWSTK) THEN NWID = 1 HFWID = DOPPH*FREQNM ELSE IF(HWLOR.GE.DOPPH.AND.HWLOR.GE.HWSTK) THEN NWID = 2 HFWID = HWLOR*FREQNM ELSE NWID = 3 HFWID = HWSTK*FREQNM END IF C ATTEMPT TO DO THIS BETTER - ASSUME FOR CORE STARK IS APPROX LORENTZIAN C AND USE VOIGT FUNCTION NWID = 1 HFWID = (DOPPH + HWLOR + HWSTK/5.) * FREQNM ! defines where we switch from Voigt to adding VA = (HWLOR + HWSTK/5.) / DOPPH ! VOIGT "a" C C HTOTAL is the profile. C HTOTAL=0.D0 C C Case 1: Doppler core including fine structure C IF ((ABS(DEL).GT.HFWID).OR.(NWID.EQ.1)) THEN DO 3 I = 1,IFINS D = DABS(FREQ-FREQNM-FINEST(I))/DOP c IF(D.LE.7.) HTOTAL = HTOTAL + EXP(-D*D)*FINSWT(I)/(SQRTPI*DOP) C C IF INSIDE HFWID WE USE VOIGT FUNCTION WITH DOP, STARK AND LOR C IF OUTSIDE JUST ADD THE DOPPLER PART - REST IS ADDED LATER IN DETAIL C IF (ABS(DEL).GT.HFWID) THEN IF(D.LE.7.) HTOTAL = HTOTAL + EXP(-D*D)*FINSWT(I)/(SQRTPI*DOP) ELSE CALL VOIGT(VA,D,XV) HTOTAL = HTOTAL + XV*FINSWT(I)/(SQRTPI*DOP) ENDIF 3 CONTINUE ENDIF C C Case 2: Self-broadening section - includes natural and vdW due to C helium for convenience. C IF ((ABS(DEL).GT.HFWID).OR.(NWID.EQ.2)) THEN C C All lines except Lyman alpha: C C The full classical expression for the line shape is used (eqn 4-73 in C Aller). This extends the profile to large detuning, but note the C impact approximation is not valid. C C The classical damping constant is frequency dependent. Assume similar C dependence of quantum mechanical damping (cf. eqn 4-116 Aller). C IF (.NOT.LYMANALF) THEN HHW = HWLOR*FREQNM HHW = HHW * FREQ*FREQ/(FREQNM*FREQNM) FREQSQ = FREQ*FREQ XSELF = 4.*FREQSQ*HHW/PI/ ; ((DEL*(FREQNM+FREQ))**2. + FREQSQ*4.*HHW**2.) C C Lyman alpha: C ELSE C C Natural broadening C C This is cut at some specified detuning in the red. C IF (WAVE.LT.RAYLCUT) THEN HHW = RADAMP*FREQNM HHW = HHW * FREQ*FREQ/(FREQNM*FREQNM) FREQSQ = FREQ*FREQ HTOTAL = HTOTAL + 4.*FREQSQ*HHW/PI/ ; ((DEL*(FREQNM+FREQ))**2. + FREQSQ*4.*HHW**2.) ENDIF C C Self broadening of Lyman alpha: C C Red wing of Lyman Alpha before the satellite region and blue wing C (ie. redward of ~1190 A). C C The impact approximation breaks down quickly here validity is ~ 1 A C at T~5000 K (Lortet & Roueff 1969, A&A 3, 462) but we use Ali & Griem C anyway (~ Barklem et al for 2p). C C A factor 1.13 is needed to match this region to the far red wing C (below) at -4000 cm^-1 (62 A) detuning. A linear transition is used. C IF (FREQ.GT.(82259.105-4000.)*CLIGHTCM) THEN XFAC = 1.+0.13*(DABS(DELW)/62.) IF (FREQ.GT.FREQNM) XFAC = 1.0 HHW = (XFAC*RESONT*H1FRC + VDW*T3NHE) * FREQNM XSELF = HHW/PI/(DEL*DEL+HHW*HHW) C C Far red wing of Lyman Alpha in satellite region: C C Self-broadening following Kurucz's implementation of Allard et al C (A&A 335, 1124, 1998), see Castelli & Kurucz (A&A 372, 260 2001) for C more details. C C Tables LYMANH2() are on a delta wavenumber grid spacing of 200 C waves/cm -22000, -21800, ...., -4000 waves/cm (91 points). The C profiles store log10(I(d omega)) for N(H) = 1e14 cm^-3. Assumed C insensitive to T, and linear scaling with N(H). Note I(d freq) = C I(d omega)/c C ELSE IF (FREQ.GT.20000.*CLIGHTCM) THEN SPACING = 200.*CLIGHTCM FREQ22000 = (82259.105-22000.)*CLIGHTCM C C If redward of last point (1660 -> 5000 Angstrom) extrapolate, C otherwise (1278 -> 1660 Angstrom) linear (in log10 profiles) C interpolation C IF (FREQ.LT.FREQ22000) THEN XLYMANH2 = (LYMANH2(2)-LYMANH2(1))/SPACING* * (FREQ-FREQ22000)+LYMANH2(1) ELSE ICUT = (FREQ-FREQ22000)/SPACING GRIDFREQ = ICUT*SPACING+FREQ22000 XLYMANH2 = (LYMANH2(ICUT+2)-LYMANH2(ICUT+1))/SPACING* * (FREQ-GRIDFREQ)+LYMANH2(ICUT+1) ENDIF XLYMANH2 = (10.**(XLYMANH2-14.))*H1FRC/CLIGHTCM XSELF = XLYMANH2 ELSE XSELF = 0. END IF C C End Lyman alpha section C ENDIF C HTOTAL = HTOTAL + XSELF C C End self-broadening section C ENDIF C C Case 3: Stark section C IF ((ABS(DEL).GT.HFWID).OR.(NWID.EQ.3)) THEN C C Stark wings due to ions and electrons. If Lyman alpha this is for C before (blueward of) the satellite region. C IF ((.NOT.(LYMANALF)).OR. ; (FREQ.GT.(82259.105-4000.)*CLIGHTCM)) THEN XSTARK = STARK1(N,M,WAVE,WAVEH,T,XNE) C C If Lyman alpha we match the static ion part to the Allard et al C value at -4000 cm^-1 detuning on red wing. C IF ((LYMANALF).AND.(FREQ.LT.(82259.105)*CLIGHTCM)) THEN WAVE4000 = 1.D8/(82259.105-4000.) XHOLT4000 = 0.5 * STARK1(N,M,WAVE4000,WAVEH,T,XNE) SAT4000 = (10.**(-11.07-14.))/CLIGHTCM*XNE XFAC = 1.+(SAT4000/XHOLT4000-1.)*DABS(DELW)/62. XSTARK = XSTARK * 0.5 * (1. + XFAC) ENDIF C C Stark wings due to protons in the satellite region: C C Following Kurucz's implementation of Allard et al (A&A 335, 1124, C 1998), see Castelli & Kurucz (A&A 372, 260 2001) for more details. C C Tables LYMANH2PLUS() are on a delta wavenumber grid spacing of 100 C waves/cm -15000, -14900, ...., -4000 waves/cm (111 points). C C The profiles store log10(I(d omega)) for N(H+) = 1e14 cm^-3. Assumed C insensitive to T, and linear scaling with N(H+). Note I(d freq) = C I(d omega)/c and we assume N(H+) = N(e-). C ELSE IF (FREQ.GT.20000.*CLIGHTCM) THEN SPACING=100.*CLIGHTCM FREQ15000=(82259.105-15000.)*CLIGHTCM C C If redward of last point (1487 -> 5000 Angstroem) extrapolate, C otherwise (1278 -> 1487 Angstroem) interpolation. C IF (FREQ.LT.FREQ15000) THEN XLYMANH2PLUS = (LYMANH2PLUS(2)-LYMANH2PLUS(1))/SPACING* ; (FREQ-FREQ15000)+LYMANH2PLUS(1) ELSE ICUT = (FREQ-FREQ15000)/SPACING GRIDFREQ = ICUT*SPACING+FREQ15000 XLYMANH2PLUS = (LYMANH2PLUS(ICUT+2)-LYMANH2PLUS(ICUT+1)) ; /SPACING*(FREQ-GRIDFREQ)+LYMANH2PLUS(ICUT+1) ENDIF XLYMANH2PLUS = (10.**(XLYMANH2PLUS-14.))*XNE/CLIGHTCM C C There may be a contribution due to electrons also at long range. C The static Holtsmark profile (1/2 the total STARK1) is used noting C that the correction for quantum effects suggested by Stehle (1994, C A&AS 104, 509 eqn 7) has been applied already. C XSTARK = XLYMANH2PLUS + 0.5 * STARK1(N,M,WAVE,WAVEH,T,XNE) ELSE XSTARK = 0. END IF HTOTAL = HTOTAL + XSTARK C C End Stark section C END IF C HPROFL = HTOTAL C RETURN END C*********************************************************************** FUNCTION STARK1(N,M,WAVE,WAVEH,T,XNE) C C Returns the Stark broadened line profile. The code follows Griem's C theories (mostly 1960 ApJ, 132, 883) with corrections to approximate C the Vidal, Cooper & Smith (1973, ApJS 25, 37) profiles. C C Area normalised to unity with frequency. C C by Deane Peterson & Bob Kurucz. C (adapted, corrected and comments added by PB) C REAL*8 WAVE,WAVEH,DELW,DEL,F,FO,CLIGHT,FREQ,FREQNM REAL*4 K DIMENSION Y1WTM(2,2),XKNMTB(4,3) LOGICAL LYMANALF SAVE C C Knm constants as defined by Griem (1960, ApJ 132, 883) for the long C range Holtsmark profile (due to ions only). Lyman and Balmer series C are from VCS, higher series from elsewhere. C DATA XKNMTB/0.0001716, 0.0090190, 0.1001000, 0.5820000, 1 0.0005235, 0.0177200, 0.1710000, 0.8660000, 2 0.0008912, 0.0250700, 0.2230000, 1.0200000/ C DATA Y1WTM/1.E18, 1.E17, 1.E16, 1.E14/ DATA N1/0/, M1/0/ C PARAMETER (CLIGHT = 2.9979258E18) PARAMETER (PI = 3.14159265359, SQRTPI = 1.77245385) PARAMETER (H = 6.62618E-27) !Planck in cgs PARAMETER (K = 1.38066E-16) !Boltzmann in cgs C C Variables depending on conditions C T4 = T/10000. T43 = T4**0.3 T3NHE = T43*HE1FRC XNE16 = XNE**0.1666667 PP = XNE16*0.08989/SQRT(T) ! the shielding parameter FO = XNE16**4*1.25E-9 ! Holtsmark normal field strength Y1B = 2./(1.+0.012/T*SQRT(XNE/T)) Y1S = T43/XNE16 C1D = FO*78940./ T C2D = FO**2/5.96E-23/XNE GCON1 = 0.2+0.09*SQRT(T4)/(1.+XNE/1.E13) GCON2 = 0.2/(1.+XNE/1.E15) C DELW = WAVE-WAVEH FREQNM = CLIGHT/WAVEH FREQ = CLIGHT/WAVE DEL = FREQ-FREQNM C C Variables dependent on line - compute first time only C IF((N.NE.N1).OR.(M.NE.M1)) THEN N1 = N M1 = M MMN = M-N XN = N XN2 = XN*XN XM = M XM2 = XM*XM XMN2 = XM2*XN2 XM2MN2 = XM2-XN2 GNM = XM2MN2/XMN2 C C Knm constants not tabulated from approximate asymptotic expression C (Griem 1960 eqn 33) where 1./(1.+.13/FLOAT(MMN)) appears to be a C correction factor to match to the tables. C IF ((MMN.LE.3).AND.(N.LE.4)) THEN XKNM = XKNMTB(N,MMN) ELSE XKNM = 5.5E-5/GNM*XMN2/(1.+.13/FLOAT(MMN)) END IF C C Some weighting factors which relate to y1, which is the velocity at C which the minimum impact parameter (where second order perturbation C theory breaks down) and the Lewis cutoff (the limit of validity of C the impact approximation) are the same. C IF(M.EQ.2) THEN Y1NUM = 550. ELSE IF (M.EQ.3) THEN Y1NUM = 380. ELSE Y1NUM = 320. END IF IF (MMN.LE.2 .AND. N.LE.2) THEN Y1WHT = Y1WTM(N,MMN) ELSE IF (MMN.LE.3) THEN Y1WHT = 1.E14 ELSE Y1WHT = 1.E13 END IF C C1CON = XKNM/WAVEH*GNM*XM2MN2 C2CON = (XKNM/WAVEH)**2 C ENDIF C C Compute line profile C C PRQS is the quasistatic ion contribution C FNS is the quasistatic electron contribution rel to PRQS C F is the impact electron contribution C C First compute the width of the impact electron profile roughly Griem C (1967, ApJ 147, 1092) eqn for w. C WTY1 = 1./(1.+XNE/Y1WHT) Y1SCAL = Y1NUM*Y1S*WTY1+Y1B*(1.-WTY1) C1 = C1D*C1CON*Y1SCAL C2 = C2D*C2CON G1 = 6.77*SQRT(C1) BETA = DABS(DELW)/FO/XKNM Y1 = C1*BETA Y2 = C2*BETA**2 IF ((Y2.LE.1.E-4).AND.(Y1.LE.1.E-5)) THEN GAM = G1*AMAX1(0.,0.2114+LOG(SQRT(C2)/C1))*(1.-GCON1-GCON2) ELSE GAM = G1*(0.5*EXP(-MIN(80.,Y1))+VCSE1F(Y1)-0.5*VCSE1F(Y2))* * (1.-GCON1/(1.+(90.*Y1)**3)-GCON2/(1.+2000.*Y1)) IF (GAM.LE.1.E-20) GAM = 0. END IF C C Compute individual quasistatic and impact profiles. C PRQS = SOFBET(BETA,PP,N,M) IF (GAM.GT.0.) THEN F = GAM/PI/(GAM*GAM+BETA*BETA) ELSE F = 0.D0 ENDIF C C Fraction of electrons which count as quasistatic. A fit to eqn 8 C (2nd term) of Griem (1967, ApJ 147, 1092). C P1 = (0.9*Y1)**2 FNS = (P1+0.03*SQRT(Y1))/(P1+1.) C C DBETA (=dBeta/dfreq) changes the area normalisation. C DSQRT(WAVE/WAVEH) corrects the long range part to dfreq**-5/2 C asymptote, (see Stehle and Hutcheon 1999, A&AS 140, 93). C DBETA = CLIGHT/FREQ/FREQ/XKNM/FO STARK1 = (PRQS*(1.+FNS)+F)*DBETA * DSQRT(WAVE/WAVEH) C C The red wing is multiplied by the Boltzmann factor to roughly account C for quantum effects (Stehle 1994, A&AS 104, 509 eqn 7). Assume C absorption case. If emission do for DEL.GT.0. C IF (DEL.LT.0.d0) STARK1 = STARK1 * DEXP(-DABS(H*DEL)/K/T) C END C*********************************************************************** FUNCTION HFNM(N,M) C C HFNM calculates hydrogen oscillator strengths C C From Kurucz codes. C SAVE DATA NSTR/0/,MSTR/0/ C HFNM=0. IF(M.LE.N) RETURN IF(N.NE.NSTR) THEN XN=N GINF=0.2027/XN**0.71 GCA=0.124/XN FKN=XN*1.9603 WTC=0.45-2.4/XN**3*(XN-1.) NSTR=N XM=M XMN=M-N FK=FKN*(XM/(XMN*(XM+XN)))**3 XMN12=XMN**1.2 WT=(XMN12-1.)/(XMN12+WTC) FNM=FK*(1.-WT*GINF-(0.222+GCA/XM)*(1.-WT)) MSTR=M ELSE IF(M.NE.MSTR) THEN XM=M XMN=M-N FK=FKN*(XM/(XMN*(XM+XN)))**3 XMN12=XMN**1.2 WT=(XMN12-1.)/(XMN12+WTC) FNM=FK*(1.-WT*GINF-(0.222+GCA/XM)*(1.-WT)) MSTR=M END IF HFNM=FNM C RETURN END C*********************************************************************** FUNCTION VCSE1F(X) C C E1 function calculator for VCS approximation. It's rough, but C arranged to be fast. X must be >=0. C C From Kurucz codes. C VCSE1F=0.0 IF(X.LE.0.0) RETURN IF(X.LE.0.01) THEN VCSE1F=-LOG(X)-0.577215+X ELSE IF(X.LE.1.0) THEN VCSE1F=-LOG(X)-0.57721566+X*(0.99999193+X*(-0.24991055+ + X*(0.05519968+X*(-0.00976004+ + X*0.00107857)))) ELSE IF(X.LE.30.) THEN VCSE1F=(X*(X+2.334733)+0.25062)/(X*(X+3.330657)+ + 1.681534)/X*EXP(-X) END IF C RETURN END C*********************************************************************** FUNCTION SOFBET(B,P,N,M) C C Calculates S(BETA,P) for Hydrogen lines ie. the Holtsmark profile for C quasistatic charged particles. The alpha and beta lines of the first C three series are explicitly included. All other cases use the H18 C profile. Profiles are normalised to full oscillator strength. Method C is based on Griem (1960, ApJ 132, 883). C C By Deane Peterson and Bob Kurucz. C C STORAGE FOR CORRECTIONS (P,BETA,IND),(P,IND),(P,IND) C DIMENSION PROPBM(5,15,7),C(5,7),D(5,7) DIMENSION PP(5),BETA(15) DIMENSION PROB1(75),PROB2(75),PROB3(75),PROB4(75),PROB5(75) DIMENSION PROB6(75),PROB7(75) DIMENSION C1(5),C2(5),C3(5),C4(5),C5(5),C6(5),C7(5) DIMENSION D1(5),D2(5),D3(5),D4(5),D5(5),D6(5),D7(5) EQUIVALENCE (PROPBM(1,1,1),PROB1(1)),(PROPBM(1,1,2),PROB2(1)) EQUIVALENCE (PROPBM(1,1,3),PROB3(1)),(PROPBM(1,1,4),PROB4(1)) EQUIVALENCE (PROPBM(1,1,5),PROB5(1)),(PROPBM(1,1,6),PROB6(1)) EQUIVALENCE (PROPBM(1,1,7),PROB7(1)) EQUIVALENCE (C(1,1),C1(1)),(C(1,2),C2(1)),(C(1,3),C3(1)) EQUIVALENCE (C(1,4),C4(1)),(C(1,5),C5(1)),(C(1,6),C6(1)) EQUIVALENCE (C(1,7),C7(1)) EQUIVALENCE (D(1,1),D1(1)),(D(1,2),D2(1)),(D(1,3),D3(1)) EQUIVALENCE (D(1,4),D4(1)),(D(1,5),D5(1)),(D(1,6),D6(1)) EQUIVALENCE (D(1,7),D7(1)) SAVE PROPBM,C,D,PP,BETA C C Lyman alpha C DATA PROB1/ 1-.980,-.967,-.948,-.918,-.873,-.968,-.949,-.921,-.879,-.821, 2-.950,-.922,-.883,-.830,-.764,-.922,-.881,-.830,-.770,-.706, 3-.877,-.823,-.763,-.706,-.660,-.806,-.741,-.682,-.640,-.625, 4-.691,-.628,-.588,-.577,-.599,-.511,-.482,-.484,-.514,-.568, 5-.265,-.318,-.382,-.455,-.531,-.013,-.167,-.292,-.394,-.478, 6 .166,-.056,-.216,-.332,-.415, .251, .035,-.122,-.237,-.320, 7 .221, .059,-.068,-.168,-.247, .160, .055,-.037,-.118,-.189, 8 .110, .043,-.022,-.085,-.147/ DATA C1 /-18.396, 84.674,-96.273, 3.927, 55.191/ DATA D1 / 11.801, 9.079, -0.651,-11.071,-26.545/ C C Lyman beta C DATA PROB2/ 1-.242, .060, .379, .671, .894, .022, .314, .569, .746, .818, 2 .273, .473, .605, .651, .607, .432, .484, .489, .442, .343, 3 .434, .366, .294, .204, .091, .304, .184, .079,-.025,-.135, 4 .167, .035,-.082,-.189,-.290, .085,-.061,-.183,-.287,-.374, 5 .032,-.127,-.249,-.344,-.418,-.024,-.167,-.275,-.357,-.420, 6-.061,-.170,-.257,-.327,-.384,-.047,-.124,-.192,-.252,-.306, 7-.043,-.092,-.142,-.190,-.238,-.038,-.070,-.107,-.146,-.187, 8-.030,-.049,-.075,-.106,-.140/ DATA C2 / 95.740, 18.489, 14.902, 24.466, 42.456/ DATA D2 / -6.665, -7.136,-10.605,-15.882,-23.632/ C C Balmer alpha C DATA PROB3/ 1-.484,-.336,-.206,-.111,-.058,-.364,-.264,-.192,-.154,-.144, 2-.299,-.268,-.250,-.244,-.246,-.319,-.333,-.337,-.336,-.337, 3-.397,-.414,-.415,-.413,-.420,-.456,-.455,-.451,-.456,-.478, 4-.446,-.441,-.446,-.469,-.512,-.358,-.381,-.415,-.463,-.522, 5-.214,-.288,-.360,-.432,-.503,-.063,-.196,-.304,-.394,-.468, 6 .063,-.108,-.237,-.334,-.409, .151,-.019,-.148,-.245,-.319, 7 .149, .016,-.091,-.177,-.246, .115, .023,-.056,-.126,-.189, 8 .078, .021,-.036,-.091,-.145/ DATA C3 /-25.088,145.882,-50.165, 7.902, 51.003/ DATA D3 / 7.872, 5.592, -2.716,-12.180,-25.661/ C C Balmer beta C DATA PROB4/ 1-.082, .163, .417, .649, .829, .096, .316, .515, .660, .729, 2 .242, .393, .505, .556, .534, .320, .373, .394, .369, .290, 3 .308, .274, .226, .152, .048, .232, .141, .052,-.046,-.154, 4 .148, .020,-.094,-.200,-.299, .083,-.070,-.195,-.299,-.385, 5 .031,-.130,-.253,-.348,-.422,-.023,-.167,-.276,-.359,-.423, 6-.053,-.165,-.254,-.326,-.384,-.038,-.119,-.190,-.251,-.306, 7-.034,-.088,-.140,-.190,-.239,-.032,-.066,-.103,-.144,-.186, 8-.027,-.048,-.075,-.106,-.142/ DATA C4 / 93.783, 10.066, 9.224, 20.685, 40.136/ DATA D4 / -5.918, -6.501,-10.130,-15.588,-23.570/ C C Paschen alpha C DATA PROB5/ 1-.819,-.759,-.689,-.612,-.529,-.770,-.707,-.638,-.567,-.498, 2-.721,-.659,-.595,-.537,-.488,-.671,-.617,-.566,-.524,-.497, 3-.622,-.582,-.547,-.523,-.516,-.570,-.545,-.526,-.521,-.537, 4-.503,-.495,-.496,-.514,-.551,-.397,-.418,-.448,-.492,-.547, 5-.246,-.315,-.384,-.453,-.522,-.080,-.210,-.316,-.406,-.481, 6 .068,-.107,-.239,-.340,-.418, .177,-.006,-.143,-.246,-.324, 7 .184, .035,-.082,-.174,-.249, .146, .042,-.046,-.123,-.190, 8 .103, .036,-.027,-.088,-.146/ DATA C5 /-19.819, 94.981,-79.606, 3.159, 52.106/ DATA D5 / 10.938, 8.028, -1.267,-11.375,-26.047/ C C Paschen beta C DATA PROB6/ 1-.073, .169, .415, .636, .809, .102, .311, .499, .639, .710, 2 .232, .372, .479, .531, .514, .294, .349, .374, .354, .279, 3 .278, .253, .212, .142, .040, .215, .130, .044,-.051,-.158, 4 .141, .015,-.097,-.202,-.300, .080,-.072,-.196,-.299,-.385, 5 .029,-.130,-.252,-.347,-.421,-.022,-.166,-.275,-.359,-.423, 6-.050,-.164,-.253,-.325,-.384,-.035,-.118,-.189,-.252,-.306, 7-.032,-.087,-.139,-.190,-.240,-.029,-.064,-.102,-.143,-.185, 8-.025,-.046,-.074,-.106,-.142/ DATA C6 /111.107, 11.910, 9.857, 21.371, 41.006/ DATA D6 / -5.899, -6.381,-10.044,-15.574,-23.644/ C C Balmer 18 C DATA PROB7/ 1 .005, .128, .260, .389, .504, .004, .109, .220, .318, .389, 2-.007, .079, .162, .222, .244,-.018, .041, .089, .106, .080, 3-.026,-.003, .003,-.023,-.086,-.025,-.048,-.087,-.148,-.234, 4-.008,-.085,-.165,-.251,-.343, .018,-.111,-.223,-.321,-.407, 5 .032,-.130,-.255,-.354,-.431, .014,-.148,-.269,-.359,-.427, 6-.005,-.140,-.243,-.323,-.386, .005,-.095,-.178,-.248,-.307, 7-.002,-.068,-.129,-.187,-.241,-.007,-.049,-.094,-.139,-.186, 8-.010,-.036,-.067,-.103,-.143/ DATA C7 /511.318, 1.532, 4.044, 19.266, 41.812/ DATA D7 / -6.070, -4.528, -8.759,-14.984,-23.956/ DATA PP/0.,.2,.4,.6,.8/ DATA BETA/1.,1.259,1.585,1.995,2.512,3.162,3.981,5.012,6.310, 1 7.943,10.,12.59,15.85,19.95,25.12/ C IF(B.GT.500.) THEN C C Very large B C B2=B*B SOFBET=(1.5/SQRT(B)+27./B2)/B2 RETURN END IF C C Other cases C CORR=1. B2=B*B SB=SQRT(B) INDX=7 MMN=M-N IF(N.LE.3.AND.MMN.LE.2) INDX=2*(N-1)+MMN C C Determine relevant Debye range C IM=MIN(INT(5.*P)+1,4) IP=IM+1 WTPP=5.*(P-PP(IM)) WTPM=1.-WTPP IF(B.LE.25.12) THEN C JP=2 1 IF(B.GT.BETA(JP).AND.JP.LT.15) THEN JP=JP+1 GO TO 1 END IF JM=JP-1 C WTBP=(B-BETA(JM))/(BETA(JP)-BETA(JM)) WTBM=1.-WTBP CBP=PROPBM(IP,JP,INDX)*WTPP+PROPBM(IM,JP,INDX)*WTPM CBM=PROPBM(IP,JM,INDX)*WTPP+PROPBM(IM,JM,INDX)*WTPM CORR=1.+CBP*WTBP+CBM*WTBM C C Get approximate profile for the inner part C PR1=0. PR2=0. WT=AMAX1(MIN(0.5*(10.-B),1.),0.) IF(B.LE.10.) PR1=8./(83.+(2.+0.95*B2)*B) IF(B.GE.8.) PR2=(1.5/SB+27./B2)/B2 SOFBET=(PR1*WT+PR2*(1.-WT))*CORR ELSE C C Asymptotic part for medium B's C CC=C(IP,INDX)*WTPP+C(IM,INDX)*WTPM DD=D(IP,INDX)*WTPP+D(IM,INDX)*WTPM CORR=1.+DD/(CC+B*SB) SOFBET=(1.5/SB+27./B2)/B2*CORR END IF C RETURN END C*********************************************************************** REAL*8 FUNCTION AIRVAC(W) REAL*8 W, WAVEN, WNEW C C W is air wavelength in Angstroms. WAVEN is air wavenumber which is C usually good enough, must iterate for exact solution C WAVEN = 1.D8 / W WNEW = W * (1.0000834213D0 + 2406030.D0 / (1.30D10 - WAVEN**2) + + 15997.D0 / (3.89D9 - WAVEN**2)) WAVEN = 1.D8 / WNEW WNEW = W * (1.0000834213D0 + 2406030.D0 / (1.30D10 - WAVEN**2) + + 15997.D0 / (3.89D9 - WAVEN**2)) WAVEN = 1.D8 / WNEW AIRVAC = W * (1.0000834213D0 + 2406030.D0 / (1.30D10 - WAVEN**2) + + 15997.D0 / (3.89D9 - WAVEN**2)) C RETURN END C*********************************************************************** REAL*8 FUNCTION VACAIR(W) REAL*8 W, WAVEN C C W is vacuum wavelength in Angstroms C WAVEN = 1.D8 / W VACAIR = W / (1.0000834213D0 + 2406030.D0 / (1.30D10 - WAVEN**2) + + 15997.D0 / (3.89D9 - WAVEN**2)) C RETURN END C*********************************************************************** REAL*8 FUNCTION WCALC(NH, NE, NHE, NS, TEMP) C C Computes the occupation probability for a H atom in C state with effective principle quantum number NS in a plasma C enviroment with NH, NE, NHE number densities of H, ions, C and He atoms respectively. This code assumes the perturbing C neutral H and He are in the ground state, (noting the hard C sphere approximation used is quite crude anyway) and that ions C are predominantly singly ionized (ie. N(Z=1 ions) = Ne). C C See eqn 4.71 Hummer & Milhalas (1988) C Ions are now treated via Appendix A Hubeny, Hummer & Lanz (1994) C which is a fit to detailed calculation including correlations, C thus the temperature dependence C C Sizes of neutral atoms adopted are sqrt() C C Coded by Paul Barklem and Kjell Eriksson, Aug 2003 C IMPLICIT NONE REAL NH, NE, NHE, NS, TEMP REAL*8 IONH, A0, E, PI, K, CHI, RIH, NEUTR, ION, NS2, NS4 REAL*8 F, A, BETAC, X, WNEUTR, WION, X1, X2 PARAMETER (IONH=2.17991E-11) PARAMETER (A0=5.29177E-9) PARAMETER (E=4.803207E-10) PARAMETER (PI=3.1415927) C NS2 = NS*NS NS4 = NS2*NS2 C C Neutral perturbers C CHI=IONH/(NS2) RIH=SQRT(2.5*NS4 + 0.5*NS2)*A0 X1=RIH + 1.73*A0 X2=RIH + 1.02*A0 NEUTR=NH*X1*X1*X1 + NHE*X2*X2*X2 WNEUTR = DEXP(-4.18879 * NEUTR) ! 4.18879 is 4*pi/3 C C Charged perturbers C K=1. IF(NS .gt. 3.) THEN C K=5.33333333*(NS/(NS+1.))**2 * (NS + 7./6.)/(NS2+NS+0.5) K=5.33333333*NS/(NS+1.)/(NS+1.) ENDIF c ION = NE*16.*(E*E/CHI/DSQRT(K))**3 c WION = DEXP(-4.*PI/3. * (ION)) IF ((NE.GT.10.).AND.(TEMP.GT.10.)) THEN ! just in case!!! A = 0.09 * DEXP(0.16667d0*LOG(NE)) / SQRT(TEMP) X = DEXP(3.15*DLOG(1.+A)) BETAC = 8.3E14 * EXP(-0.66667*LOG(NE)) * K / NS4 F = 0.1402d0*X*BETAC*BETAC*BETAC / / (1.+0.1285*X*BETAC*DSQRT(BETAC)) WION = F/(1.d0+F) ELSE WION = 1.0d0 ENDIF C WCALC = WION * WNEUTR RETURN END C*********************************************************************** REAL FUNCTION HBF(WAVE, NLO) C C Returns the bound-free cross section due to H in level NLO C at wavelength WAVE (in Angstroms). C C It gives an extrapolation redward of usual bound free jumps C so that we can use the occupation probability formalism C for "broadening" the jump due to nearby particles C C Data (Gaunt factors - Karzas & Latter) from MARCS code C INTEGER N(15), NLO REAL L1(12), G1(12), L2(12), G2(12), L3(14), G3(14) REAL L4(14), G4(14), L5(14), G5(14), L6(14), G6(14) REAL L7(15), G7(15), L8(15), G8(15), L9(6), G9(6) REAL L10(16), G10(16), L11(6), G11(6), L12(6), G12(6) REAL L13(6), G13(6), L14(6), G14(6), L15(6), G15(6) REAL G, W REAL*8 WAVE C DATA N/ 12, 12, 14, 14, 14, 14, 15, 15, 6, 16, 6, 6, 6, 6, 6/ C DATA L1/ 182., 241., 300., 356., 408., 456., * 538., 631., 729., 821., 877., 912./ DATA G1/ .9939, .9948, .9856, .9721, .9571, .9423, * .9157, .8850, .8531, .8246, .8076, .7973/ DATA L2/ 215., 301., 398., 503., 614., 729., * 965., 1313., 1824., 2525., 3144., 3647./ DATA G2/ 1.043, 1.058, 1.063, 1.061, 1.056, 1.049, * 1.033, 1.009, .9755, .9342, .9011, .8761/ DATA L3/ 222., 316., 424., 545., 677., 821., * 1132., 1641., 2525., 4103., 6034., 6933., * 7529., 8207./ DATA G3/ 1.053, 1.072, 1.081, 1.083, 1.081, 1.078, * 1.068, 1.051, 1.024, .9854, .9458, .9293, * .9189, .9075/ DATA L4/ 224., 321., 434., 561., 703., 859., * 1205., 1799., 2918., 5252., 8895., 10998., * 12576., 14588./ DATA G4/ 1.056, 1.077, 1.087, 1.091, 1.091, 1.089, * 1.082, 1.069, 1.047, 1.014, .9736, .9542, * .9408, .9247/ DATA L5/ 226., 324., 438., 569., 715., 877., * 1241., 1882., 3144., 6034., 11397., 15095., * 18235., 22794./ DATA G5/ 1.058, 1.080, 1.090, 1.095, 1.095, 1.094, * 1.089, 1.078, 1.060, 1.030, .9925, .9719, * .9563, .9358/ DATA L6/ 226., 325., 441., 573., 722., 887., * 1262., 1931., 3282., 6564., 13448., 18916., * 24120., 32915./ DATA G6/ 1.059, 1.081, 1.092, 1.097, 1.098, 1.097, * 1.092, 1.083, 1.067, 1.041, 1.006, .9851, * .9682, .9436/ DATA L7/ 227., 326., 442., 576., 726., 894., * 1275., 1961., 3372., 6933., 15095., 22347., * 29992., 36616., 44694./ DATA G7/ 1.060, 1.082, 1.093, 1.098, 1.099, 1.099, * 1.095, 1.086, 1.071, 1.047, 1.015, .9952, * .9776, .9641, .9494/ DATA L8/ 227., 326., 443., 578., 729., 897., * 1284., 1982., 3433., 7196., 16398., 25326., * 35615., 45361., 58446./ DATA G8/ 1.060, 1.082, 1.094, 1.099, 1.100, 1.100, * 1.096, 1.088, 1.074, 1.052, 1.022, 1.003, * .9853, .9708, .9540/ DATA L9/ 227., 901., 7383., 40703., 67537., 74126./ DATA G9/ 1.060, 1.101, 1.056, .9917, .9632, .9576/ DATA L10/ 227., 327., 445., 580., 732., 903., * 1294., 2006., 3507., 7529., 18235., 29992., * 45588., 63316., 72940., 91175./ DATA G10/ 1.060, 1.083, 1.095, 1.100, 1.102, 1.101, * 1.098, 1.091, 1.078, 1.058, 1.032, 1.014, * .9970, .9814, .9737, .9606/ DATA L11/ 82., 905., 7636., 49822., 96995.,111189./ DATA G11/ .9379, 1.102, 1.060, 1.002, .9703, .9632/ DATA L12/ 82., 905., 7720., 53950.,112562.,132138./ DATA G12/ .9379, 1.102, 1.062, 1.005, .9732, .9653/ DATA L13/ 82., 906., 7793., 57343.,130250.,154534./ DATA G13/ .9380, 1.102, 1.063, 1.008, .9758, .9672/ DATA L14/ 82., 907., 7846., 60381.,147056.,178775./ DATA G14/ .9380, 1.103, 1.064, 1.011, .9780, .9689/ DATA L15/ 82., 908., 7887., 63316.,162813.,207216./ DATA G15/ .9380, 1.103, 1.064, 1.014, .9803, .9703/ C W = WAVE G = 0. IF (NLO.EQ.1) CALL LOGINT(L1,G1,N(1),W,G) IF (NLO.EQ.2) CALL LOGINT(L2,G2,N(2),W,G) IF (NLO.EQ.3) CALL LOGINT(L3,G3,N(3),W,G) IF (NLO.EQ.4) CALL LOGINT(L4,G4,N(4),W,G) IF (NLO.EQ.5) CALL LOGINT(L5,G5,N(5),W,G) IF (NLO.EQ.6) CALL LOGINT(L6,G6,N(6),W,G) IF (NLO.EQ.7) CALL LOGINT(L7,G7,N(7),W,G) IF (NLO.EQ.8) CALL LOGINT(L8,G8,N(8),W,G) IF (NLO.EQ.9) CALL LOGINT(L9,G9,N(9),W,G) IF (NLO.EQ.10) CALL LOGINT(L10,G10,N(10),W,G) IF (NLO.EQ.11) CALL LOGINT(L11,G11,N(11),W,G) IF (NLO.EQ.12) CALL LOGINT(L12,G12,N(12),W,G) IF (NLO.EQ.13) CALL LOGINT(L13,G13,N(13),W,G) IF (NLO.EQ.14) CALL LOGINT(L14,G14,N(14),W,G) IF (NLO.EQ.15) CALL LOGINT(L15,G15,N(15),W,G) C HBF = 1.044E-26 * G * WAVE**3.d0 / (NLO**5.d0) RETURN END C*********************************************************************** SUBROUTINE LOGINT(X, Y, N, XI, YI) C C Does log interpolation and extrapolation of arrays X and Y C of length N. Returns value YI at XI. C IMPLICIT NONE INTEGER N, I, J REAL X(*), Y(*), XI, YI REAL XL(N), YL(N), XIL, YIL, XLS(3), YLS(3), DY C J = 0 XIL = LOG10(XI) DO 10 I = 1, N XL(I) = LOG10(X(I)) YL(I) = LOG10(Y(I)) IF ((XL(I).LT.XIL) .AND. J.EQ.0) J = I 10 CONTINUE IF (J.EQ.0) J = 1 IF (J.GT.N-3) J = N-3 C C polynomial 3 point interpolation or extrapolation C XLS(1) = XL(J) XLS(2) = XL(J+1) XLS(3) = XL(J+2) YLS(1) = YL(J) YLS(2) = YL(J+1) YLS(3) = YL(J+2) CALL POLINT(XLS, YLS, 3, XIL, YIL, DY) YI = 10.**YIL C RETURN END C*********************************************************************** C Numerical recipes polynomial interpolation routine SUBROUTINE POLINT(XA,YA,N,X,Y,DY) INTEGER I,N,NMAX,M,NS REAL XA(*),YA(*),X,Y,DY,DIF,DIFT,HO,HP,W,DEN PARAMETER (NMAX=10) DIMENSION C(NMAX),D(NMAX) NS=1 DIF=ABS(X-XA(1)) DO 11 I=1,N DIFT=ABS(X-XA(I)) IF (DIFT.LT.DIF) THEN NS=I DIF=DIFT ENDIF C(I)=YA(I) D(I)=YA(I) 11 CONTINUE Y=YA(NS) NS=NS-1 DO 13 M=1,N-1 DO 12 I=1,N-M HO=XA(I)-X HP=XA(I+M)-X W=C(I+1)-D(I) DEN=HO-HP IF(DEN.EQ.0.d0)PAUSE DEN=W/DEN D(I)=HP*DEN C(I)=HO*DEN 12 CONTINUE IF (2*NS.LT.N-M)THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY 13 CONTINUE RETURN END C***********************************************************************