c This program computes RT CIA for H2_H2 pairs c at temperatures between 600 and 7400K c Note: this program does not handle non-equilibrium H2. c C ============================================================================ C Copyright (C) 1993 Chunguang Zheng and Aleksandra Borysow C ============================================================================ C Copyright Notice: C You may use this program for your scientific applications, C but please do not distribute it yourself. C Direct all your inquires to aborysow@phy.mtu.edu C Final request: If you publish your work which benefited from this C program, please ACKNOWLEDGE using this program and QUOTE C the original paper describing the procedure: c c author = " C. Zheng and A. Borysow", c title = "Rototranslational {CIA} spectra of {H$_2$-H$_2$} c at temperatures between 600 and 7,000~{K}", c journal = "Astrophys.\ J.", c pages = " 960-965 ", year = " 1995", volume = " 441" c C ============================================================================ C BLOCK DATA INPUT implicit double precision (a-h,o-z) dimension fit(2,3,5),P(0:3,0:2),Q(2,0:3,0:2) COMMON /BLOCKIN/ fit COMMON /PQ/P,Q C*********************************************************************** C Q is quadrupole & P is polarizability for H2. C j is rotational quantum number. vib is vibrational quantum number. C Q(dj,vib,i) are fitted parameters for Q=Q0+Q1*J(J+1)+Q2*(J(J+1))^2 C dj=1,2 for delta(j)=0,2. C P(vib,i) are fitted parameters for P=P0+P1*J(J+1)+P2*(J(J+1))^2 C delta(j)=0. C*********************************************************************** data (Q(1,0,i),i=0,2)/0.96683,0.00113673,-3.14771e-07/ data (Q(1,1,i),i=0,2)/1.06931,0.00118117,-5.49e-07/ data (Q(1,2,i),i=0,2)/1.17193,0.00113273,-8.2635e-07/ data (Q(1,3,i),i=0,2)/1.26542,0.000770076,9.36233e-08/ data (P(0,i),i=0,2)/5.41315,0.00474039,-6.10846e-07/ data (P(1,i),i=0,2)/5.88304,0.0050382, -7.66904e-07/ data (P(2,i),i=0,2)/6.37242,0.00531034,-9.71942e-07/ data (P(3,i),i=0,2)/6.8769,0.00554729,-1.23633e-06/ data (Q(2,0,i),i=0,2)/0.974011,0.0012223, -5.25832e-07/ data (Q(2,1,i),i=0,2)/1.07511,0.00109792,-7.24724e-07/ data (Q(2,2,i),i=0,2)/1.17736,0.000789709,-7.19323e-07/ data (Q(2,3,i),i=0,2)/1.26615,0.000223587,4.02125e-07/ C*********************************************************************** C fit(m,n,i) are fitted parameters for M0, Tau1, Tau2. C m=1 for 2021 term and m=2 for 2023 term C n=1, 2, 3 for M0, Tau1, Tau2 respectively C polynomial fits were used for M0-Temp and log10(Tau1,Tau2)-log10(Temp) C*********************************************************************** data (fit(1,1,i),i=1,5)/-1.15381e-63, 5.09166e-66, + 1.43838e-69,-2.56086e-74,0.0/ data (fit(1,2,i),i=1,5)/-13.1952,-0.0080046, + 0.00373508,-0.020582,0.0/ data (fit(1,3,i),i=1,5)/42.0169,-66.6006, + 30.184,-6.16637,0.478008/ data (fit(2,1,i),i=1,5)/1.53992e-62,2.59111e-65, + -2.79125e-70,4.41603e-74,0.0/ data (fit(2,2,i),i=1,5)/-7.2311,-4.91122, + 1.31344,-0.133142,0.0/ data (fit(2,3,i),i=1,5)/-15.9277,2.72415, + -0.911435,0.0789137,0.0/ C********************************************************************** END PROGRAM H2_RTCIA implicit double precision (a-h,o-z) dimension alfa(601),freq(601),abs1(601),fit(2,3,5) COMMON /BLOCKIN/ fit COMMON /PRTSUM/ Q,W(2),JRANGE1 data w/1.0,3.0/ write(*,*)' Input MINIMUM, MAXIMUM & INTERVAL FREQUENCY in cm-1' read (*,*) FNUMIN,FNUMAX,DNU NF=INT((FNUMAX-FNUMIN)/DNU+0.5)+1 IF (NF.GT.601) then NF=601 FNUMAX=FNUMIN+DFLOAT(NF-1)*DNU write(*,*)' Number of intervals is larger than 600 and is changed +to 600.','MAXIMUM FREQUENCY is now ',FNUMAX end if write (*,*) ' Input temperature 600--7400K' read (*,*) temp if(temp.lt.600. .or. temp.gt.7400. ) then write(*,*) ' Temperature out of range!!!!! ' stop 999 else end if write(*,*) ' Your results in file: out ' open(unit=3, file='out', status='unknown') CALL PARTFCT (TEMP, Q, 0, JRANGE1, W ) DO 10 I=1,NF FREQ(I)=FNUMIN+DFLOAT(I-1)*DNU 10 alfa(i)=0.d0 x=dlog10(temp) s1=fit(1,1,1)+fit(1,1,2)*temp+fit(1,1,3)*temp**2 + +fit(1,1,4)*temp**3+fit(1,1,5)*temp**4 t11=10.0**(fit(1,2,1)+fit(1,2,2)*x+fit(1,2,3)*x*x + +fit(1,2,4)*(x**3)+fit(1,2,5)*(x**4)) t21=10.0**(fit(1,3,1)+fit(1,3,2)*x+fit(1,3,3)*x*x + +fit(1,3,4)*(x**3)+fit(1,3,5)*(x**4)) s2=fit(2,1,1)+fit(2,1,2)*temp+fit(2,1,3)*temp**2 + +fit(2,1,4)*temp**3+fit(2,1,5)*temp**4 t12=10.0**(fit(2,2,1)+fit(2,2,2)*x+fit(2,2,3)*x*x + +fit(2,2,4)*(x**3)+fit(2,2,5)*(x**4)) t22=10.0**(fit(2,3,1)+fit(2,3,2)*x+fit(2,3,3)*x*x + +fit(2,3,4)*(x**3)+fit(2,3,5)*(x**4)) do 500 NV1=0,3 do 500 NV2=0,3 CALL ADDSPEC(S1,T11,T21,TEMP,2,0,2,1,NV1,NV2,NF,FREQ,ABS1,1) do 296 l=1,nf 296 alfa(l)=alfa(l)+2.*abs1(l) CALL ADDSPEC(S2,T12,T22,TEMP,2,0,2,3,NV1,NV2,NF,FREQ,ABS1,2) do 196 l=1,nf 196 alfa(l)=alfa(l)+2.*abs1(l) 500 CONTINUE write(3,140) temp 140 format('## RT-CIA SPECTRA alfa(omega) of H2-H2 ',/, + 2h##,5x,'temperature =',f8.1,' K',/,'##',2x, + 'freq(cm-1)',5x,'alfa(cm-1amagat-2)') do 195 l=1,nF 195 write(3, 160) freq(l), alfa(l) 160 format ( f10.1,8x,e14.4) close (3) END SUBROUTINE ADDSPEC (G0,TAU1,TAU2,TEMP,LAMBDA1,LAMBDA2, & LAMBDA,LVALUE,NVIB1,NVIB2,NF,FREQ,ABSCOEF,jp) C C THIS PROGRAM GENERATES A LISTING OF THE CIA TR ALFA(OMEGA) C IF BOTH LAMBDA1 AND LAMBDA2 ARE NEGATIVE: SINGLE TRANSITIONS; C DOUBLE TRANSITIONS ARE ASSUMED OTHERWISE. C LIKE=1 FOR LIKE SYSTEMS (AS H2-H2), SET LIKE=0 ELSE. C implicit double precision (a-h,o-z) dimension FREQ(601),ABSCOEF(601),P(0:3,0:2),Q(2,0:3,0:2) COMMON /PRTSUM/ QQ,WH2(2),JRANGE1 COMMON /PQ/P,Q DATA B0,D0/59.3392,0.04599/ DATA CLOSCHM,BOLTZWN/2.68675484E19,.6950304/ DATA HBAR,PI,CLIGHT/1.054588757d-27,3.1415926535898,2.9979250E10/ do 777 i=1, nf 777 abscoef(i) = 0.d0 PIC=PI*CLIGHT CALIB=PIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*CLOSCHM**2 BETA=1./(BOLTZWN*TEMP) LIST=NF C******ROTATIONAL SPECTRUM FOR THE DETAILED LISTING******************* JPLUSL=JRANGE1+MAX0(LAMBDA1,LAMBDA2) DO 50 I1=1,JRANGE1 J1=I1-1 DO 50 IP1=1,JPLUSL JP1=IP1-1 CG1S=CLEBSQR(J1,LAMBDA1,JP1) IF (CG1S) 50,50,10 10 P1 =H2POPL(nvib1,J1,TEMP,WH2)/QQ IF (P1.LT.0.001) GO TO 50 OMEGA1=H2ELEV(nvib1,JP1)-H2ELEV(nvib1,J1) DO 40 I2=1,JRANGE1 J2=I2-1 DO 40 IP2=1,JPLUSL JP2=IP2-1 CG2S=CLEBSQR(J2,LAMBDA2,JP2) IF (CG2S) 40,40,20 20 P2 =H2POPL(nvib2,J2,TEMP,WH2)/QQ IF (P2.LT.0.001) GO TO 40 OMEGA2=H2ELEV(nvib2,JP2)-H2ELEV(nvib2,J2) FAC=CALIB*P1*P2*CG1S*CG2S DO 30 I=1,LIST FRQ=FREQ(I)-OMEGA1-OMEGA2 WKI=FREQ(I)*(1.-DEXP(-BETA*FREQ(I))) WKF=WKI*FAC if (jp.eq.1) then XBG=G0*BGAMA0(FRQ,TAU1,TAU2,TEMP) else XBG=G0*BGAMA1(FRQ,TAU1,TAU2,TEMP) end if IF (jp.ne.2) then CP = 1.0 CQ = 1.0 GOTO 170 end if if (jp1-j1) 110,120,130 110 a=dfloat(jp1) ndj=2 goto 100 120 a=dfloat(j1) ndj=1 goto 100 130 a=dfloat(j1) ndj=2 100 cq=(q(ndj,nvib1,0)+q(ndj,nvib1,1)*a*(a+1) + +q(ndj,nvib1,2)*(a*(a+1))**2)/q(2,0,0) a=dfloat(j2) cp=(p(nvib2,0)+p(nvib2,1)*a*(a+1) + +p(nvib2,2)*(a*(a+1))**2)/p(0,0) 170 ABSCOEF(I)=ABSCOEF(I)+XBG*WKF*(cp*cq)**2 30 CONTINUE 40 CONTINUE 50 CONTINUE RETURN END C*************************************************************** FUNCTION CLEBSQR (L,LAMBDA,LP) C C SQUARE OF CLEBSCH-GORDAN COEFFICIENT (L,LAMBDA,0,0;LP,0) C FOR INTEGER ARGUMENTS ONLY C NOTE THAT LAMBDA SHOULD BE SMALL, MAYBE @10 OR SO. C implicit double precision (a-h,o-z) FC=DFLOAT(2*LP+1) GO TO 10 C ENTRY THREEJ2 C C THIS ENTRY RETURNS THE SQUARED 3-J SYMBOL L LAMBDA LP C 0 0 0 C INSTEAD OF THE CLEBSCH-GORDAN COEFFICIENT C (LIMITATION TO INTEGER ARGUMENTS ONLY) C C NOTE THAT THE THREE-J SYMBOLS ARE COMPLETELY SYMMETRIC IN THE C ARGUMENTS. IT WOULD BE ADVANTAGEOUS TO REORDER THE INPUT ARGUMENT C LIST SO THAT LAMBDA BECOMES THE SMALLEST OF THE 3 ARGUMENTS. C FC=1. 10 CLEBSQR=0. IF (((L+LAMBDA).LT.LP).OR.((LAMBDA+LP).LT.L).OR.((L+LP).LT.LAMBDA) 1) RETURN IF (MOD(L+LP+LAMBDA,2).NE.0) RETURN IF ((L.LT.0).OR.(LP.LT.0).OR.(LAMBDA.LT.0)) RETURN F=1./DFLOAT(L+LP+1-LAMBDA) IF (LAMBDA.EQ.0) GO TO 30 I1=(L+LP+LAMBDA)/2 I0=(L+LP-LAMBDA)/2+1 DO 20 I=I0,I1 20 F=F*DFLOAT(I)/DFLOAT(2*(2*I+1)) 30 P=FC*F*FCTL(LAMBDA+L-LP)*FCTL(LAMBDA+LP-L) CLEBSQR=P/(FCTL((LAMBDA+L-LP)/2)*FCTL((LAMBDA+LP-L)/2))**2 RETURN C END FUNCTION FCTL (N) implicit double precision (a-h,o-z) P(Z)=((((-2.294720936d-4)/Z-(2.681327160d-3))/Z+(3.472222222d-3))/ 1Z+(8.333333333d-2))/Z+1. FCTL=1. IF (N.LE.1) RETURN IF (N.GT.15) GO TO 20 J=1 DO 10 I=2,N 10 J=J*I FCTL=DFLOAT(J) RETURN 20 Z=DFLOAT(N+1) FCTL=(DEXP(-Z)*(Z**(Z-0.5))*P(Z)*2.5066282746310) RETURN END CZ G0*bgama0(,,,) or G0*bgama1(...), G0 is the zeroth moment (=M0) FUNCTION BGAMA1(FNU,TAU1,TAU2,TEMP) C BIRNBAUM S CIA LINE SHAPE MODEL (K1) C NORMALIZATION SUCH THAT 0TH MOMENT EQUALS UNITY C FNU IS THE FREQUENCY IN CM-1; TEMP IN KELVIN. implicit double precision (a-h,o-z) DATA TWOPIC,BKW,HBOK,PI/1.883651568d11,.6950304256d0, 1 7.638280918d-12, 3.141592654d0/ TAU3=dSQRT(TAU2*TAU2+(HBOK/(TEMP*2.))**2) OMEGA=TWOPIC*FNU DENOM=1.d0+(OMEGA*TAU1)**2 X=(TAU3/TAU1)*dSQRT(DENOM) AAA=TAU2/TAU1 AAA=dMIN1(AAA,430.d0) BGAMA1=(TAU1/PI)*dEXP(AAA+FNU/(2.d0*BKW*TEMP))* 1 XK1(X)/DENOM RETURN END FUNCTION BGAMA0(FNU,TAU5,TAU6,TEMP) C K0 LINE SHAPE MODEL C NORMALIZATION SUCH THAT 0TH MOMENT EQUALS UNITY C FNU IS THE FREQUENCY IN CM-1; TEMP IN KELVIN. implicit double precision (a-h,o-z) DATA TWOPIC,BKW,HBOK,PI/1.883651568d11,.6950304256d0, 1 7.638280918d-12, 3.141592654d0/ TAU4=dSQRT(TAU5*TAU5+(HBOK/(TEMP*2.d0))**2) OMEGA=TWOPIC*FNU XNOM=1.d0/(TAU6*TAU6)+OMEGA*OMEGA X=TAU4*dSQRT(XNOM) TAU56=TAU5/TAU6 TAU56=dMIN1(TAU56,430.d0) BGAMA0=(TAU5/PI)*dEXP(TAU56+FNU/(2.d0*BKW*TEMP))* 1 XK0(X) RETURN END FUNCTION XK0(X) C MODIFIED BESSEL FUNCTION K0(X) C ABRAMOWITZ AND STEGUN P.379 implicit double precision (a-h,o-z) IF(X-2.) 10,10,20 10 T=(X/3.75d0)**2 FI0=(((((.0045813*T+.0360768)*T+.2659732)*T 1 +1.2067492)*T+3.0899424)*T+3.5156229)*T+1. T=(X/2.)**2 P=(((((.00000740*T+.00010750)*T+.00262698)*T 1 +.03488590)*T+.23069756)*T+.42278420)*T+ 2 (-.57721566) X=ABS(X) XK0=-dLOG(X/2.)*FI0+P RETURN 20 T=(2./X) P=(((((.00053208*T-.00251540)*T+.00587872)*T 1 -.01062446)*T+.02189568)*T-.07832358)*T+ 2 1.25331414 X=dMIN1(X,330.d0) XK0=dEXP(-X)*P/dSQRT(X) RETURN END FUNCTION XK1(X) C MODIFIED BESSEL FUNCTION K1(X) TIMES X C PRECISION IS BETTER THAN 2.2e-7 EVERYWHERE. C ABRAMOWITZ AND S,TEGUN, P.379; TABLES P.417. implicit double precision (a-h,o-z) IF(X-2.) 10,10,20 10 T=(X/3.75)**2 FI1=X*((((((.00032411*T+.00301532)*T+.02658733)*T+.15084934) 1 *T+.51498869)*T+.87890594)*T+.5) T=(X/2.)**2 P=(((((-.00004686*T-.00110404)*T-.01919402)*T-.18156897)*T- 1 .67278579)*T+.15443144)*T+1. XK1=X*dLOG(X/2)*FI1+P RETURN 20 T=2./X P=(((((-.00068245*T+.00325614)*T-.00780353)*T+.01504268)*T- 1 .03655620)*T+.23498619)*T+1.25331414 X=dMIN1(X,330.d0) XK1=dSQRT(X)*dEXP(-X)*P RETURN END SUBROUTINE PARTFCT (TEMP, Q, NTYPE, JRANGE, W) C C NTYPE = 0 for equilibrium hydrogen C IMPLICIT double precision (A-H,O-Z) DIMENSION W(2) Q=0. J=0 30 DQ=H2POPL(0,J,TEMP,W) + h2popl(1,j,temp,w) + 1 h2popl(2,j,temp,w) + h2popl(3,j,temp,w) CZ 1 h2popl(4,j,temp,w) + h2popl(5,j,temp,w) Q = Q + DQ J = J + 1 IF (DQ.GT.Q/500.) GO TO 30 JRANGE=J RETURN END FUNCTION H2POPL (NVIB, J, T, W) C UNNORMALIZED (!) Boltzmann factor of the roto- vibrational levels C of the H2 molecule. (To be normalized by Q = partition sum) C IMPLICIT double precision (A-H,O-Z) DIMENSION W(2) DATA BOLTZWN /.6950304D0/ C 20 DDD = H2ELEV(nvib,J) H2POPL = DFLOAT(2*J+1)*W(1+MOD(J,2))*DEXP(-DDD/(BOLTZWN*T)) RETURN END FUNCTION H2ELEV (NVIB, JROT) C Roto- vibrational energy levels of the H2 molecule C Needed for the computation of the partition sum Q of H2 ****** C IMPLICIT double precision (A-H,O-Z) C c The following results are from LEVELS (levels.lev file) of c Kolos et all. (HSINGLX.FOR function), fitted as function of J(J+1) c Fitting program FITLEVELS.FOR (.for), below only for v=0, 1, 2 C NOTE: ROTATIONAL CONSTANTS B AND D HAVE BEEN FIXED WHILE FITTING THE C ROTATIONAL LEVELS, AFTER S.L.BRAGG, J.W.BRAULT,W.H.SMITH, AP.J.,VOL.263, C P.999-1004, (1982). EH2(I, A, B, D, H, XL, G, P, R)= A + B*DFLOAT(I) - 1 D*1.D-2*DFLOAT(I)**2 + H * 1.D-5 * DFLOAT(I)**3 - 1 XL*1.D-8*DFLOAT(I)**4 + G * 1.D-11*DFLOAT(I)**5 - 1 P * 1.D-14 * DFLOAT(I)**6 + R*1.D-17*DFLOAT(I)**7 C A - THE ROTATIONAL ENERGY (NVIB, J=0) c |E(v=0,J=0) - E(v=1,J=0)| = 4162.14 cm-1 if(nvib.eq.0) h2elev=eh2(jrot*(jrot+1), -0.361132D5, 1 59.33451D0, 4.56510D0, 4.568777D0, 4.366951D0, 2.863721D0, 1 1.051977D0, 0.156903D0) + 36113.d0 c add 36113.d0 to each energy to prevent h2elev/(kT) being too c big at small temperatures c ====================================================== if(nvib.eq.1) h2elev=eh2(jrot*(jrot+1), -0.319511D5, 1 56.3742D0, 4.4050D0, 4.515341D0, 4.614924D0, 3.301549D0, 1 1.32896D0, 0.212922D0) + 36113.d0 if(nvib.eq.2) h2elev=eh2(jrot*(jrot+1),-0.280238D5, 53.482D0, 1 4.28D0, 4.552623D0, 4.957454D0, 3.695921D0, 1.469646D0, 1 0.199078D0) + 36113.d0 if(nvib.eq.3) h2elev=eh2(jrot*(jrot+1), -.243268d+05, 1 0.5050d+02, 0.3818d+01, 0.2393d+01, -.8313d+00, -.4652d+01, 1 -.4660d+01, -.1628d+01) + 36113.d0 if(nvib.eq.4) h2elev=eh2(jrot*(jrot+1), -.208566d+05, 1 0.4762d+02, 0.3605d+01, 0.1511d+01, -.3982d+01, -.1106d+02, 1 -.1108d+02, -.4167d+01) + 36113.d0 if(nvib.eq.5) h2elev=eh2(jrot*(jrot+1), -.176133d+05, 1 0.4481d+02, 0.3511d+01, 0.1099d+01, -.6643d+01, -.1913d+02, 1 -.2174d+02, -.9433d+01) + 36113.d0 if(nvib.gt.5) stop 888 C RETURN END