* Date: Sun, 31 May 1998 15:10:50 +0200 (MET DST) * From: Bengt Edvardsson * To: Plez@Ferrum.fysik.lu.se FUNCTION HLINOP(WAVE,NBLO,NBUP,WAVEH,T,XNE, * HNORM,H1FRC,HE1FRC,DOPPLE) C C Hydrogen line opacity C * version received Aug 97 from Martin Asplund * hlinop = resulting opacity in cm2/gram stellar matter * wave = opacity interrogation wavelength in A real*8 * nblo = lower level main quantum number * nbup = upper level main quantum number * waveh = hydrogen line center wavelength in A real*8 * T = kinetic temperature in K * xne = electron number density in cm-3 * hnorm = sqrt(pi) e^2/(me c) * gf * 1/ntot * nHI/Q(HI) * exp(-elow/kT) * * Pg/kT * 1/rho * (1. -exp-(h nu/kT)) * h1frc = number density of H I in cm-3 * he1frc = number density of He I in cm-3 C DOPPLE is Doppler width in delta_lambda/lambda * REAL KAP0,KAP0RED,KAPRED,KAP0BLUE,KAPBLUE REAL*8 WAVEH,WAVE,WCON,EMERGE,WMERGE,WSHIFT,WTAIL REAL*8 BLUECUT,WLBLUE,REDCUT,WLRED,EMERGEH REAL*8 EHYD(100),CONTH(15) REAL*8 VACAIR LOGICAL FIRST SAVE DATA CONTH/ 1 109678.764,27419.659,12186.462,6854.871,4387.113, 2 3046.604,2238.320,1713.711,1354.044,1096.776, 3 906.426,761.650,648.980,559.579,487.456/ DATA FIRST/.TRUE./ C C Lyman series is treated in vacuum. The rest is converted to air C IF(NBLO.EQ.1) THEN IFVAC=1 ELSE IFVAC=0 END IF IF(FIRST) THEN cc IFVAC=1 EHYD(1)= 0.000D0 EHYD(2)= 82259.105D0 EHYD(3)= 97492.302D0 EHYD(4)=102823.893D0 EHYD(5)=105291.651D0 EHYD(6)=106632.160D0 EHYD(7)=107440.444D0 EHYD(8)=107965.051D0 DO 1 I=9,100 1 EHYD(I)=109678.764D0-109677.576D0/I**2 FIRST=.FALSE. END IF INGLIS=1194./XNE**.125 NMERGE=INGLIS-1.5 EMERGE= 109737.312D0/NMERGE**2 EMERGEH=109677.576D0/NMERGE**2 C C C Cutoff when the energy difference reaches ionization level. C This fraction of opacity was taken care of in hydrogen b-f opacity C calculations. C IF(NBUP.GT.NBLO+2) THEN C C Only blue lines can reach ionization energy C WSHIFT=1.D8/(CONTH(NBLO)-109677.576D0/81.D0**2) WMERGE=1.D8/(CONTH(NBLO)-EMERGEH) IF(WMERGE.LT.0.) WMERGE=WSHIFT+WSHIFT WCON=MAX(WSHIFT,WMERGE) WTAIL=1.D8/(1.D8/WCON-500.D0) WCON=MIN(WSHIFT+WSHIFT,WCON) IF(WTAIL.LT.0.) WTAIL=WCON+WCON WTAIL=MIN(WCON+WCON,WTAIL) IF(IFVAC.EQ.0) THEN WCON=VACAIR(WCON) WTAIL=VACAIR(WTAIL) ENDIF IF(WAVE.GE.WAVEH) THEN C C Red wing C REDCUT=WAVEH IF(NBUP.LE.10) THEN REDCUT=1.D8/(109678.764D0-109677.576D0/(NBUP-0.8D0)**2- - EHYD(NBLO)) IF(IFVAC.EQ.0) REDCUT=VACAIR(REDCUT) ENDIF WLRED=1.D8/(EHYD(NBUP-2)-EHYD(NBLO)) IF(IFVAC.EQ.0) WLRED=VACAIR(WLRED) HLINOP=0.0 IF(WAVE.LT.WCON) RETURN HLINOP=HPROFL(NBLO,NBUP,WAVEH,WAVE-WAVEH, * T,XNE,HNORM,H1FRC,HE1FRC,DOPPLE,KAP0) HLINOP=HLINOP*KAP0 IF(WAVE.LT.WTAIL) HLINOP=HLINOP*(WAVE-WCON)/(WTAIL-WCON) KAP0RED=KAP0*HFNM(NBLO,NBUP-2)/HFNM(NBLO,NBUP)/ / (EHYD(NBUP-2)-EHYD(NBLO))*(EHYD(NBUP)-EHYD(NBLO)) IF(WAVE.GT.REDCUT) THEN KAPRED=HPROFL(NBLO,NBUP-2,WAVEH,WAVE-WLRED, * T,XNE,HNORM,H1FRC,HE1FRC,DOPPLE,DUMMY) KAPRED=KAPRED*KAP0RED IF(WAVE.LT.WTAIL) * KAPRED=KAPRED*(WAVE-WCON)/(WTAIL-WCON) IF(HLINOP.LT.KAPRED) HLINOP=0.0 END IF ELSE C C Blue wing C BLUECUT=WAVEH IF(NBUP.LE.10) THEN BLUECUT=1.D8/(109678.764D0-109677.576D0/(NBUP+0.8D0)**2- - EHYD(NBLO)) IF(IFVAC.EQ.0) BLUECUT=VACAIR(BLUECUT) ENDIF WLBLUE=1.D8/(EHYD(NBUP+2)-EHYD(NBLO)) IF(IFVAC.EQ.0) WLBLUE=VACAIR(WLBLUE) HLINOP=0.0 IF(WAVE.LT.WCON.OR.WAVE.LT.WLBLUE) RETURN HLINOP=HPROFL(NBLO,NBUP,WAVEH,WAVE-WAVEH, * T,XNE,HNORM,H1FRC,HE1FRC,DOPPLE,KAP0) HLINOP=HLINOP*KAP0 IF(WAVE.LT.WTAIL) HLINOP=HLINOP*(WAVE-WCON)/(WTAIL-WCON) KAP0BLUE=KAP0*HFNM(NBLO,NBUP+2)/HFNM(NBLO,NBUP)/ / (EHYD(NBUP+2)-EHYD(NBLO))*(EHYD(NBUP)-EHYD(NBLO)) IF(WAVE.LT.BLUECUT) THEN KAPBLUE=HPROFL(NBLO,NBUP+2,WAVEH,WAVE-WLBLUE, * T,XNE,HNORM,H1FRC,HE1FRC,DOPPLE,DUMMY) KAPBLUE=KAP0BLUE*KAPBLUE IF(WAVE.LT.WTAIL) * KAPBLUE=KAPBLUE*(WAVE-WCON)/(WTAIL-WCON) IF(HLINOP.LT.KAPBLUE) HLINOP=0.0 ENDIF END IF ELSE C C Well isolated lines in the red C HLINOP=HPROFL(NBLO,NBUP,WAVEH,WAVE-WAVEH, * T,XNE,HNORM,H1FRC,HE1FRC,DOPPLE,KAP0) HLINOP=HLINOP*KAP0 END IF c correct normalization for higher lines according to N. Piskunov c multiply with pi^1/2=1.77245385 c MA970819 c Removed 15/07-1998 BPz. C Put back 19/09-2001 BPz if (nbup.ge.(nblo+2)) hlinop=hlinop*1.77245385 c c IF(NBUP.GT.NBLO+1) HLINOP=HLINOP*1.77245385 c IF(NBUP.GT.NBLO+1) HLINOP=HLINOP*2 C RETURN END FUNCTION HPROFL(N,M,WAVEH,DELW,T,XNE, * HNORM,H1FRC,HE1FRC,DOPPLE,KAPPA0) C C VERSION FINE STRUCTURE LIKE GENERAL BUT APPROXIMATELY INCLUDES FINE C STRUCTURE IN THE DOPPLER CORES. EXACT PATTERN C IS USED FOR ALPHA LINES, M INFINITE PATTERN C IS USED FOR ALL OTHER LINES. C from Deane Peterson C requires VCSE1F and SOFBET C REAL*8 DELW,WAVE,WAVEH,DOP,D,FREQ,FREQNM,FINEST(14),FINSWT(14),XN2 REAL KAPPA0 DIMENSION Y1WTM(2,2),XKNMTB(4,3) DIMENSION STCOMP(5,4),STALPH(34),ISTAL(4),LNGHAL(4),STWTAL(34), 1 STCPWT(5,4),LNCOMP(4) SAVE DATA XKNMTB/0.0001716,0.0090190,0.1001000,0.5820000,0.0005235, 1 0.0177200,0.1710000,0.8660000,0.0008912,0.0250700, 2 0.2230000,1.0200000/ DATA Y1WTM/1.E18,1.E17,1.E16,1.E14/ DATA N1/0/,M1/0/,RYDH/3.2880515E15/ C C Fine structure components for alpha lines in FREQ*10**-7 C DATA STALPH/-730.,370.,188.,515.,327.,619.,-772.,-473.,-369.,120., 1 256.,162.,285.,-161.,-38.3,6.82,-174.,-147.,-101.,-77.5,55.,126., 2 75.,139.,-60.,3.7,27.,-69.,-42.,-18.,-5.5,-9.1,-33.,-24./ C C Alpha component weights C DATA STWTAL/1.,2.,1.,2.,1.,2.,1.,2.,3.,1.,2.,1.,2.,1.,4.,6.,1.,2., 1 3.,4.,1.,2.,1.,2.,1.,4.,6.,1.,7.,6.,4.,4.,4.,5./ DATA ISTAL/1,3,10,21/, LNGHAL/2,7,11,14/ C C Fine structure for M.EQ.INFINITY IN FREQ*10**-7 C DATA STCOMP/0.,0.,0.,0.,0.,468.,576.,-522.,0.,0.,260.,290.,-33., 1 -140.,0.0,140.,150.,18.,-27.,-51./ C C Weights C DATA STCPWT/1.,0.,0.,0.,0.,1.,1.,2.,0.,0.,1.,1.,4.,3.,0.,1.,1., 1 4.,6.,4./ DATA LNCOMP/1,3,4,5/ PARAMETER (SQRTPI=1.77245385,CLIGHT=2.997925E18) C C Set up some variables depending on current physical conditions C TK=1.38054E-16*T TKEV=8.6171E-5*T XNE16=XNE**0.1666667 PP=XNE16*0.08989/SQRT(T) FO=XNE16**4*1.25E-9 Y1B=2./(1.+0.012/T*SQRT(XNE/T)) T4=T/10000. T43=T4**0.3 Y1S=T43/XNE16 T3NHE=T43*HE1FRC 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 C DOPPLE is Doppler width in delta_lambda/lambda C C DOP is Doppler width in delta_nu C DOP=DOPPLE*CLIGHT/WAVEH C C DOPPH Doppler width in delta_lambda/lambda C DOPPH=DOPPLE C C Normalization factor for opacity: C C HNORM=Z*ABUND*GF*XNFPH*BOLTH*XNATOM*XSTIM/RHO C Z=sqrt(pi)*e^2/(m*c) is in cm^2/s, C ABUND is fractional abundance relative to the total number of atoms, C GF is the oscillator strength, C XNFPH is the fraction of neutral hydrogen divided by partition function, C BOLTH=exp(-Elow/kT) where Elow is in ergs, C XNATOM is the total number of atoms in cm^-3 and C XSTIM=1-exp(-h*nu/kT), stimulated emission C RHO is the density in g*cm^-3. C Doppler width is in frequency units. C The units of KAPPA0 are cm^2/g. C KAPPA0=HNORM/DOP C C Set up for this line 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 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 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 FREQNM=RYDH*GNM DBETA=CLIGHT/FREQNM**2/XKNM WAVE=CLIGHT/FREQNM C1CON=XKNM/WAVE*GNM*XM2MN2 C2CON=(XKNM/WAVE)**2 RADAMP=1.389E9/XM**4.53/(1.+5./XM2/XM) IF(N.NE.1) RADAMP=RADAMP+1.389E9/XN**4.53/(1.+5./XN2/XN) RADAMP=RADAMP/FREQNM RESONT=HFNM(1,M)/XM/(1.-1./XM2) IF(N.NE.1) RESONT=RESONT+HFNM(1,N)/XN/(1.-1./XN2) C C Fudge to Baschek*2: C 2 IS FOR CONVERTING XNFPH TO XNFH: C RESONT=RESONT*5.593E-24/GNM*2. C RESONT=RESONT*5.593E-24/GNM C error in constant corrected 26nov95 C RESONT=RESONT*3.92E-24/GNM VDW=4.45E-26/GNM*(XM2*(7.*XM2+5.))**.4 STARK=1.6678E-18*FREQNM*XKNM C C Fine structure components C IF(N.GT.4.OR.M.GT.10) THEN IFINS=1 FINEST(1)=0. FINSWT(1)=1. ELSE IF(MMN.GT.1) THEN C C Use M.EQ.INF structure C 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 For alpha lines 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 given P,Pe,T conditions C DEL=-DELW/WAVE*FREQNM FREQ=FREQNM+DEL C C These half-widths: DOPPH, HWSTK, HWLOR, are in fact DNU/NU C HWSTK=STARK*FO HWLOR=RESONT*H1FRC+VDW*T3NHE+RADAMP C C Specify the largest half-width in case of core calculations 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 C Set flag if in line core C HPROFL=0. DOP=FREQNM*DOPPH C C Two big cases: first, we are close to the core C IF(ABS(DEL).LE.HFWID) THEN IF(NWID.EQ.1) THEN C C Compute Doppler core including fine structure C DO 3 I=1,IFINS D=ABS(FREQ-FREQNM-FINEST(I))/DOP C C Same normalization as the Voigt function C IF(D.LE.7.) HPROFL=HPROFL+EXP(-D*D)*FINSWT(I) 3 CONTINUE ELSE IF(NWID.EQ.2) THEN C C Compute Lorentz wings C IF(N.EQ.1.and.m.eq.2.AND.FREQ.LE.1.8474E15) THEN C C Insert Lyman alpha red cutoff a la Sando and Wormhoudt (1973 PHYS REV A 7,1889) C HWLOR=VDW*T3NHE+RADAMP FF=1.-FREQ/1.8474E15 * following statement never executed because of above if (suppressed by * BPz 3/06-1998, following MA version) cc IF(FREQ.GT.1.6E16) HWLOR=RESONT*H1FRC/(1.+EXP(-4492./T)* cc * SQRT(FF)*EXP(1398./T**0.333333*FF))+HWLOR HWLOR=RESONT*H1FRC/(1.+EXP(-4492./T)* * SQRT(FF)*EXP(1398./T**0.333333*FF))+HWLOR END IF HHW=FREQNM*HWLOR HPROFL=HPROFL+HHW/SQRTPI/(DEL**2+HHW**2)*DOP ELSE IF(NWID.EQ.3) THEN C C Compute Stark wings 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) GNOT=G1*MAX(0.,0.2114+LOG(SQRT(C2)/C1))*(1.-GCON1-GCON2) BETA=ABS(DEL)/FO*DBETA Y1=C1*BETA Y2=C2*BETA**2 IF(Y2.LE.1.E-4.AND.Y1.LE.1.E-5) THEN GAM=GNOT 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 PRQS=SOFBET(BETA,PP,N,M) F=0. IF(GAM.GT.0.) F=GAM/3.14159/(GAM**2+BETA**2) P1=(0.9*Y1)**2 FNS=(P1+.03*SQRT(Y1))/(P1+1.) C C Same normalization as all opacities C HPROFL=HPROFL+(PRQS*(1.+FNS)+F)/FO*DBETA*SQRTPI*DOP END IF ELSE C C Far off the core just add up all the components C C Compute Doppler core including fine structure C DO 4 I=1,IFINS D=ABS(FREQ-FREQNM-FINEST(I))/DOP C C Same normalization as all the opacities C IF(D.LE.7.) HPROFL=HPROFL+EXP(-D*D)*FINSWT(I) 4 CONTINUE C C Compute Lorentz wings C IF(N.EQ.1.and.m.eq.2.AND.FREQ.LE.1.8474E15) THEN C C Insert Lyman alpha red cutoff a la Sando and Wormhoudt (1973 PHYS REV A 7,1889) C HWLOR=VDW*T3NHE+RADAMP FF=1.-FREQ/1.8474E15 * following statement never executed because of above if (suppressed by * BPz 3/06-1998, following MA version) cc IF(FREQ.GT.1.6E16) HWLOR=RESONT*H1FRC/(1.+EXP(-4492./T)* cc * SQRT(FF)*EXP(1398./T**0.333333*FF))+HWLOR HWLOR=RESONT*H1FRC/(1.+EXP(-4492./T)* * SQRT(FF)*EXP(1398./T**0.333333*FF))+HWLOR END IF HHW=FREQNM*HWLOR HPROFL=HPROFL+HHW/SQRTPI/(DEL**2+HHW**2)*DOP C C Compute Stark wings 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) GNOT=G1*MAX(0.,0.2114+LOG(SQRT(C2)/C1))*(1.-GCON1-GCON2) BETA=ABS(DEL)/FO*DBETA Y1=C1*BETA Y2=C2*BETA**2 IF(Y2.LE.1.E-4.AND.Y1.LE.1.E-5) THEN GAM=GNOT 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 PRQS=SOFBET(BETA,PP,N,M) F=0. IF(GAM.GT.0.) F=GAM/3.14159/(GAM**2+BETA**2) P1=(0.9*Y1)**2 FNS=(P1+.03*SQRT(Y1))/(P1+1.) C C Same normalization as all opacities C HPROFL=HPROFL+(PRQS*(1.+FNS)+F)/FO*DBETA*SQRTPI*DOP END IF C c write(*,*) DELW,HPROFL RETURN END FUNCTION HFNM(N,M) C C HFNM calculates hydrogen oscillator strengths C SAVE DATA NSTR/0/,MSTR/0/ C HFNM=0. IF(M.LE.N) RETURN IF(N.NE.NSTR) THEN XN=N GINF=.2027/XN**.71 GCA=.124/XN FKN=XN*1.9603 WTC=.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-(.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-(.222+GCA/XM)*(1.-WT)) MSTR=M END IF HFNM=FNM C RETURN END FUNCTION VCSE1F(X) C C E1 function calculator for VCS approximation. C It's rough, but but arranged to be fast. X must be >=0. 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 FUNCTION SOFBET(B,P,N,M) C C GENERATES S(BETA,P) FOR HYDROGEN LINES. THE ALPHA AND BETA LINES C OF THE FIRST THREE SERIES ARE EXPLICITLY INCLUDED AND THE H18 C PROFILE IS USED FOR THE REST. C C THESE PROFILES HAVE BEEN RENORMALIZED TO FULL OSCILLATOR STRENGTH 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 PROB1,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, -.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 CORR=1. B2=B*B SB=SQRT(B) IF(B.GT.500.) THEN C C Very large B C SOFBET=(1.5/SB+27./B2)/B2*CORR RETURN END IF 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=MIN0(INT(5.*P)+1,4) IP=IM+1 WTPP=5.*(P-PP(IM)) WTPM=1.-WTPP IF(B.LE.25.12) THEN c c NP: my version of the loop commented out below 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 c This piece of junk form the original code can be used in books c c DO 1 J=2,15 c IF(B.LE.BETA(J)) GO TO 2 c 1 CONTINUE c 2 JM=J-1 c JP=J 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 aproximate profile for the inner part C PR1=0. PR2=0. WT=MAX(MIN(0.5*(10.-B),1.),0.) IF(B.LE.10.) PR1=8./(83.+(2.+.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 FUNCTION EXPI(N,X) C C LOW PRECISION VERSION 1.E-5 C EXPONENTIAL INTEGRAL FOR POSITIVE ARGUMENTS AFTER CODY AND C THACHER, MATH. OF COMP.,22,641(1968) C SAVE X1,EX,EX1,A0,A1,A2,B0,B1,C0,C1,C2,D1,D2,E0,E1,F1 DATA X1/-1.E20/ DATA A0,A1,A2,B0,B1/-4.43668255,4.42054938,3.16274620,7.68641124, 1 5.65655216/ DATA C0,C1,C2,D1,D2/.0012102205,.98147989,.75339742,1.6198645, 1 .29135151/ DATA E0,E1,F1/-.9969698,-.4257849,2.318261/ C IF(X.EQ.X1) THEN EX=EXP(-X) X1=X IF(X.GT.4.0) THEN EX1=(EX+EX*(E0+E1/X)/(X+F1))/X ELSE IF(X.GT.1.0) THEN EX1=EX*(C2+(C1+C0*X)*X)/(D2+(D1+X)*X) ELSE IF(X.GT.0.0) THEN EX1=(A0+(A1+A2*X)*X)/(B0+(B1+X)*X)-LOG(X) ELSE EX1=0.0 END IF END IF C EXPI=EX1 IF(N.GT.1) THEN N1=N-1 DO 1 I=1,N1 1 EXPI=(EX-X*EXPI)/FLOAT(I) END IF C RETURN END REAL*8 FUNCTION AIRVAC(W) REAL*8 W,WAVEN,WNEW C C W IS AIR WAVELENGTH IN ANGSTROMS C WAVEN IS AIR WAVENUMBER WHICH IS USUALLY GOOD ENOUGH C 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 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