C This program computes second overtone CIA spectrum of c H2-H2 pairs in the spectral range from 11,000 to 13,800 cm(-1) c at 10 cm(-1) intervals. C C ============================================================================ C Copyright (C) 1999 by Aleksandra Borysow, Jacek Borysow and Yi Fu C ============================================================================ 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 one of the authors: e-mail aborysow@.nbi.dk C Final request: If you publish your work which benefited from this C program, please acknowledge using this program C by referring to: C A. Borysow, J. Borysow and Y. Fu, C "Semiempirical Model of Collision Induced Infrared Absorption C Spectra of H_2-H_2 Complexes in the Second Overtone C Band of Hydrogen at Temperatures from 50 to 500 K", Icarus, 1999, submitted, C C ============================================================================ c C INPUT....## program reads file ``par.beta'' c ## program requests ratio of para to total H2 c ## 0.25 corresponds to ``normal hydrogen'' c ## by default computations are done at c ## thermal equilibrium c ## program requests temp (range is from 50 to 500 K) c c ========================================================================== c c The program will calculate WRONG absorption coefficient c for temperatures above 500 K because it assumes that c only ground vibrational states of H2 are populated c=========================================================================== c c OUTPUT is directed into file ``output.H2H2.03'' and contains absorption c coefficient divided by density squared multiplied by factor 10(8) (i.e. *1E8) c c Authors DO NOT guarantee that this program will produce c any useful results of any kind therefore everybody should c use it at her/his own risk. c PROGRAM H2H2model03_LT IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MARKER, labela, labelb dimension ibgama(6),LA1(6),LA2(6),LA(6),LL(6) dimension NVIB(10, 4), w(2), Npar(6) Dimension BM0(12,4), BTAU1(12,4), BTAU2(12,4) dimension alfatot(0:12, 600) common/ integ / rlo, rup, frac, Ninteg common/temp/temp COMMON /PRTSUM/ q1, w, NTYPE, JRANGE1, MARKER COMMON /SPECIT/ LIST common/ specit1/ FREQ(600), ABSCOEF(600) common/ read/ go common/om/omega_0 COMMON/VIBIF/ NVIB1, NVIB1p, NVIB2, NVIB2p COMMON /bpar/ Bpar(9, 4), NB DATA NTYPE, LIKE/ 0, 1 / DATA W(1), W(2)/1.0, 3.0/ ! rotational weights; data nf/ 600 / data NVIB/0, 0, 0, 1, 1, 1, 2, 2, 2, 3, ! == V1 & 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, ! == V1' & 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, ! == V2 & 3, 2, 3, 3, 2, 3, 3, 2, 3, 3/ ! == V2' C IBGAMA = 0 uses K0, IBGAMA = 1 uses BC, BGAMA = 2 uses K1 data ibgama/0,1,1,2,0,1/ ! IBGAMA = 0 uses K0, IBGAMA = 1 uses BC data LA1 /0,2,0,2,0,2/ data LA2 /0,0,2,0,2,2/ data LA /0,2,2,2,2,3/ data LL /1,3,3,1,1,3/ data Npar /7,4,4,7,7,4/ c c DATA BM0/-64.9252, -64.8355, -64.7565, -70.1186, -71.3592,0.d0, & 0.0d0,-64.5756,-64.6270,-65.8634,-65.1641,-65.4075, !B0 & -3.79051, -3.18043, -2.90725, 1.67553, 7.4384, 0.d0, & 0.0d0,-2.93532,-2.97008,-3.26495,-4.97074,-2.84093, !B1 & 1.69722, 1.22891, 1.1752, -1.51323, -4.11896, 0.d0, & 0.0d0, 1.13392, 1.15348, 1.30463, 2.15157,1.09193, !B2 & -0.237117, -0.120397, -0.156739, 0.37531, 0.769848, 0.d0, & 0.0d0,-0.123712,-0.127921,-0.145022,-0.244661,-0.122333/ !B3 DATA BTAU1/-11.3308, -12.9510, -10.6268, -14.4108, -13.2326,0.d0, & 0.0d0,-11.9383,-11.7162,-12.0872,-13.5173,-11.7599, !B0 & -1.19800, -0.469693, -2.77937, 3.35085, 1.09208,0.d0, & 0.0d0,-0.595612,-0.741825,-0.251777,1.39961,-0.712328, !B1 & 0.604604, -0.0484486, 1.12265, -2.03575, -0.809052,0.d0, & 0.0d0,0.0174677,0.115793,-0.140516,-1.01606,0.0745054, !B2 & -0.140778, 0.00534555, -0.16634, 0.327846, 0.141261,0.d0, & 0.0d0,-0.00576595,-0.0310663,0.0232436,0.183536,-0.0117238/ !B3 DATA BTAU2/-20.4644,-11.87746,-12.5763,0.0d0,-22.6701,0.d0, & 0.0d0,-12.61095,-13.0725,0.0d0,-9.8152,-13.1469, !B0 & 12.8185, -0.1798,0.297043, 0.0d0, 17.8713,0.d0, & 0.d0,-0.179975,0.0680707,0.d0,-4.75287,-0.159049, !B1 & -7.74865,-0.134118,-0.267489,0.d0,-9.97786,0.d0, & 0.d0,-0.129425,-0.291291,0.d0,2.21046,-0.129466, !B2 & 1.47426,0.0170518,0.00264082, 0.d0,1.72716,0.d0, & 0.d0,0.0152609,0.0527647,0.d0,-0.376787,0.0104684/ !B3 open(unit=10,file='output.H2H2.03', status='unknown') !ouput open(unit=8, file='par.beta', status='old') !input par for Beta_(R) | Bpar open(unit=13,file='param.out',status='unknown') write(*,*) 'The program will compute absorption coefficient' write(*,*) 'of RV CIA of pure hydrogen, at temperatures ' write(*,*) 'between 50 and 500K, in the second overtone band' write(*,*) 'at any given ortho-to-para H2 ratio' write(*,*) write (*,*) 'Final results: in the file: output.H2H2.03 ' write(*,*) write(*,*) 'Temperature? ' read (*,*) temp c write(*,*) 'input minimum, maximum, step of frequency' c read(*,*) FREQLO,FREQMAX,DELTASTEP FREQLO=11000.0 FREQMAX=13800.0 DELTASTEP=10.0 list = INT (FREQMAX - FREQLO)/DELTASTEP DO 10 I=1,list 10 FREQ(I)=FREQLO+DFLOAT(I-1)*DELTASTEP C+++++++++++++ Compute the partition function Q_total: +++++++++++++++++++ CALL PARTFCT0 (TEMP, Q_tot, Q_even, JRANGE1, MARKER, W ) c Q_ratio=Q_even/Q_tot c Q11=Q_tot Q2=Q_ratio Q3=Q_even c print *,'Equilibrium Hydrogen' print *,'Temp=',temp print *,'Partition function ratio, Q_para(J=even)/Q_total', & Q_even/Q_tot print *,'Write new ratio in decimal format' print *,'write "0" for no change (computations for H2 at ', & 'equilibrium)' print *,'hint: 0.25 corresponds to "normal hydrogen"' c c 111 continue print *,'note: Ratio must be smaller than 1.0' read *,Q_ratio if (q_ratio.ge.1.0) go to 111 if (Q_ratio.eq.0.0) go to 123 Q2=Q_ratio c c c CALL weightFCT (TEMP, Q11, Q2, Q3, MARKER, W ) c c Q11 is a new nonequlibrium partition function c c comes here for equilibrium 123 continue c Q1=Q11 write(10 , 1212) temp, freqlo, freqmax, deltastep 1212 format('## 2nd overtone band RV CIA spectra of H2-H2 ',/, & '## Temperature = ', f9.2, ' K', /, & '## Frequency: Min, Max, Step:',3f9.1, ' cm-1',/ & '## freq(cm-1) alfa*1E+8 (cm-1amagat-2)',/) write (10, 1213) Q_ratio 1213 format ('## Ratio of para to total hydrogen = ',f9.2,/) c c ALFA =1./(0.69519*TEMP) hbar = 1.054588757D-27 boltzk = 1.38054d-16 tau0 = hbar/(2.*boltzk*temp) C***************************************************************** TLG = DLOG10(TEMP) do 11 N_VIB = 1, 2 ! LOOP for Vibrational Transitions NVIB1 = NVIB(N_VIB, 1) NVIB1p = NVIB(N_VIB, 2) NVIB2 = NVIB(N_VIB, 3) NVIB2p = NVIB(N_VIB, 4) do 11 N_term = 1, 6 ! LOOP for La1La2LaL terms NB = Npar(N_term) read(8,*) labelb do N_B = 1, 4 ! LOOP for reading B_par of J_Dependance read(8,*) (Bpar(i, N_B), i = 1, 9) end do if ( Bpar(1, 1) .eq. 0.0 ) goto 11 IF(N_VIB.EQ.2 .AND. N_term.EQ.1) GOTO 11 GM0 = 0.0d0 GTau_1 = 0.0d0 GTau_2 = 0.0d0 Np_term = N_term+ (N_VIB-1)*6 DO J=1,4 GM0 = GM0 + BM0(Np_term,J)* TLG **(J-1) GTau_1 = GTau_1 + BTAU1(Np_term,J)*TLG **(J-1) GTau_2 = GTau_2 + BTAU2(Np_term,J)*TLG **(J-1) ENDDO XM0 = 10**GM0 Tau1 = 10**GTau_1 Tau2 = 10**GTau_2 write(13,*) gm0,gtau_1,gtau_2 write(*, *) 'N_VIB =, N_term = ', N_VIB, N_term C write(*, *)'M0', XM0, Tau1, Tau2 c * * * * ** * * * * * * * * * * ** * * * * * * * * * * ** * * * * * CALL ADDSPEC(TEMP,LIKE, LA1(N_term), LA2(N_term), & LA(N_term), LL(N_term), freqcut, & ibgama(N_term),XM0,tau1,tau2) do 198 l=1,list alfatot(0, l) = alfatot(0, l) + ABSCOEF(l) * 2 ! total alfatot((N_VIB-1)*6 + N_term, l) = ABSCOEF(l) * 2 ! separate terms 198 continue 9997 format(1x,f10.2, 2x, e12.4) 11 continue write (10,9997) (freq(l), alfatot(0,l)*1E+8, l=1, list) END SUBROUTINE ADDSPEC(TEMP, LIKE,LAMBDA1,LAMBDA2,LAMBDA,LVALUE, & CUT, ibgama, G0, tau1,tau2) C C This subroutine generates a detailed listing of CIA ALFA(OMEGA) c IBGAMA=0 means uses K0, IBGAMA=1 (uses BC), IBGAMA=2 (uses K1) C LIKE=1 for like systems (as H2-H2) C IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MESSAGE, MARKER COMMON /PRTSUM/ q1, w1(2), NTYPE, JRANGE1, MARKER COMMON /SPECIT/ LIST common/specit1/ FREQ(600),ABSCOEF(600) common/ integ / rlo, rup, frac, ninteg COMMON/VIBIF/ NVIB1, NVIB1p, NVIB2, NVIB2p COMMON/ROTIF/ J1, JP1, J2, JP2 data nf /600 / DATA SCALEF/1.D80/, BOLTZWN/ 0.6950304D0/ DATA HBAR,PI,CLIGHT/1.054588757D-27, 3.141592653589797D0, 1 2.9979250D10/ DATA Boltzk / 1.380658d-16 / DO 10 I=1, nf 10 ABSCOEF(I)=0.D0 TWOPIC=2.*PI*CLIGHT MESSAGE='NON-COMUL.' CALIB=TWOPIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*(2.68676D19)**2 CALIB=CALIB/DFLOAT(1+LIKE) FREQCUT = CUT ! do not alter the outside parameter CUT BETA1 = 1.D0/(BOLTZWN*TEMP) JPLUSL=JRANGE1+MAX0(LAMBDA1,LAMBDA2) DO 160 I1=1,JRANGE1 J1=I1-1 DO 160 IP1=1,JPLUSL JP1=IP1-1 CG1S=CLEBSQR(J1,LAMBDA1,JP1) IF (CG1S) 160,160,110 110 P1=H2POPL(NVIB1,J1,TEMP,W1)/Q1 OMEGA1=H2ELEV(NVIB1p,JP1)-H2ELEV(NVIB1,J1) DO 150 I2=1,JRANGE1 J2=I2-1 DO 150 IP2=1,JPLUSL JP2=IP2-1 CG2S=CLEBSQR(J2,LAMBDA2,JP2) IF (CG2S) 150,150,120 120 P2=H2POPL(NVIB2,J2,TEMP,W1)/Q1 OMEGA2=H2ELEV(NVIB2p,JP2)-H2ELEV(NVIB2,J2) FAC=CALIB*P1*P2*CG1S*CG2S FAC=FAC*scalej() ! DO 140 I=1,LIST FRQ=FREQ(I)-OMEGA1-OMEGA2 130 WKF=FREQ(I) * (1.-DEXP(- BETA1 * FREQ(I))) * FAC if (ibgama.eq.0) xx = g0 * bgama0(frq,tau1,tau2,temp) if (ibgama.eq.1) xx = g0 * bgama1(frq,tau1,tau2,temp) if (ibgama.eq.2) xx = g0 * bgama_H(frq,tau1,temp) ABSCOEF(I) = ABSCOEF(I) + xx * WKF 140 CONTINUE 150 CONTINUE 160 CONTINUE RETURN END SUBROUTINE PARTFCT0 (TEMP, Q, Q_even, JRANGE, MARKER, W) C c this function calculates total partition function Q for c hydrogen and partition function for even rotational levels c only. Thermal equilibrium is assumed. c IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MARKER DIMENSION W(2) DATA BOLTZWN /.6950304D0/ Q = 0.0 Qeven = 0.0 J = -1 70 J=J+1 DDD = H2POPL(0,J,TEMP,W) + h2popl(1,j,temp,w) + 1 h2popl(2,j,temp,w) + h2popl(3,j,temp,w) Q = Q + DDD IF (MOD(J,2).eq.0) Qeven = Qeven + DDD IF (DDD.GE.0.002*Q) GO TO 70 JRANGE=J+2 MARKER='HYDROGEN ' 200 continue Q_even=Qeven RETURN END SUBROUTINE weightFCT (TEMP, Q11, Q2, Q3, MARKER, W) C c this function caclulates rotational weights for c for new ratio of Q_even/Q_total (para/total) hydrogen c + new total c partition function Q_tot C IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MARKER DIMENSION W(2) DATA BOLTZWN /.6950304D0/ C Q_tot=Q11 Q_ratio=Q2 Q_even=Q3 W(1)=(Q_tot-Q_even)/Q_even*(Q_ratio/(1.0D0-Q_ratio)) Q_tot=W(1)*Q_even +(Q_tot-Q_even) Q11=Q_tot 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 This function is for low temperatures. C C All D, B, H and wavenumber are from the reference by C Uwe Fink T.A. Wiggins and D.H. Rank, JMS 18, 384-395 (1965) C Checked by Yi Fu, Jan 29, 1998 IMPLICIT double precision (A-H,O-Z) EH2(I,B, D, H)=B*DFLOAT(I) -D*DFLOAT(I)**2 + H *DFLOAT(I)**3 if(nvib.eq.0) h2elev=eh2(jrot*(jrot+1),59.3362D0, & 4.583D-2, 4.87D-5) if(nvib.eq.1) h2elev=eh2(jrot*(jrot+1),56.3722D0, & 4.3880D-2, 3.97D-5) + 4161.1782d0 C 4161.1782cm^-1 is the wavenumber of the vibrational C transition of v= 0-> 1 if(nvib.eq.2) h2elev=eh2(jrot*(jrot+1),53.4823D0, & 4.272D-2,3.97D-5) + 8087.d0 C 8087.000d0 cm^-1 is the wavenumber of the vibrational C transition of v= 0 -> 2 if(nvib.eq.3) h2elev=eh2(jrot*(jrot+1),50.6245D0, & 4.08d-2, 3.68d-5) + 11782.355d0 C 11782.355d0cm^-1 is the wavenumber of the vibrational C transition of v= 0-> 3 if(nvib.gt.3) stop 888 C 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=dABS(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 FUNCTION BGAMA0(FNU,TAU1,TAU2,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/ TAU0 = HBOK /(2.*TEMP) TAU4=DSQRT(TAU1*TAU1+(TAU0)**2) OMEGA=TWOPIC*FNU XNOM=1.d0/(TAU2*TAU2)+OMEGA*OMEGA X=TAU4*dSQRT(XNOM) TAU12=TAU1/TAU2 TAU12=DMIN1(TAU12,430.d0) BGAMA0=(TAU1/PI)*DEXP(TAU12+FNU/(2.d0*BKW*TEMP))* 1 XK0(X) end FUNCTION CLEBSQR (L,LAMBDA,LP) 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 REAL*8 (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 END FUNCTION FCTL (N) C Simple FACTORIALS program IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) P(Z)=((((-2.294720936D-4)/Z-(2.681327160D-3))/Z+(3.472222222D-3))/ 1Z+(8.3333333333D-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 FUNCTION BGAMA1(FNU,TAU1,TAU2,TEMP) C BIRNBAUM - COHEN (BC) 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 BGAMA_H(FNU,tau ,TEMP) c "K1" one parameter model lineshape (14-Jul-1989 ) c with Egelstaff time, by A. Borysow and J. Borysow, 1989 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,.6950304256, 1 7.638280918D-12, 3.141592654/ TAU3=DSQRT( 3./2.*TAU**2 + (HBOK/(TEMP*2.))**2 ) OMEGA=TWOPIC*FNU ARG = DABS(OMEGA)*TAU3 IF(ARG.EQ.0.D0) GO TO 5 BGAMA_H =1./PI * DEXP(FNU/(2.*BKW*TEMP))* XK1(ARG) 1 * (TAU/TAU3)**2 * TAU * (3./2.)**(3./2.) RETURN 5 BGAMA_H =1./PI * DEXP(FNU/(2.*BKW*TEMP)) 1 * (TAU/TAU3)**2 * TAU * (3./2.)**(3./2.) C LIMIT K1(X) (X-->0) = 1/X; SO C LIMIT X*K1(X) (X-->0) = X/X = 1. EXACT c G1 = (beta *hbar)/(tau**2) c G2 = (5./3.)*beta*beta*hbar*hbar/(tau**4) + 2./(tau**2) c where beta = 1./(boltzk * temp) in cgs units RETURN END SUBROUTINE DERIVAS(FCT,X0,DX1, Y,Y1,Y2,Y3,Y4,IFLAG) C COMPUTES OTH..4TH DERIVATIVES OF Y=FCT(X) AT X=X0 USING FIVE C ABSCISSAE AND CENTRAL DIFFERENCES. FCT MUST BE DECLARED EXTERNAL C IN CALLING ROUTINE. NOTE THAT FCT(X) MUST BE DEFINED FOR C X0-2DX<=X<=X0+2DX C A TEST OF THE DEFINITION OF THE FIRST DERIVATIVE: IFLAG=1 OUTPUT C SEE ABRAMOWITZ AND STEGUN, P. 914 IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION F(7) 9 X=X0 H=DX1 HSQ=H*H IFLAG=0 DO 10 I=1,7 10 F(I)=FCT(X-DFLOAT(4-I)*H) Y=F(4) A=F(1)+F(7) A1=F(1)-F(7) B=F(2)+F(6) B1=F(2)-F(6) C=F(3)+F(5) C1=F(3)-F(5) C THESE ARE HIGHEST ORDER DERIVATIVES Y1 = (-A1+9.*B1-45.*C1)/(60.*H) Y2 = (A-13.5*B+135.*C-245.*Y)/(90.*HSQ) Y3 = (A1-8.*B1+13.*C1)/(8.*HSQ*H) Y4 = (-A+12.*B-39.*C+56.*Y)/(6.*HSQ*HSQ) RETURN END function scalej() IMPLICIT double precision (A-H,O-Z) external va common /temp/temp DATA ANGAU, FOURPI, HBAR/ 0.52917706, 12.566370614359, & 1.05458876D-34/ T= TEMP * 1.380662D-23 FMP=2.0*1.67265D-27 F2=HBAR*HBAR/(12.*FMP*T*T*ANGAU**2*1.D-20) h=0.05 g0=0.0 g1=0.0 do 100 i=30,160 r=h*dfloat(i) ! in A.U. dr=r/100.0 CALL BETAJ(R, B0, BJ) CALL DERIVAS(VA,R,DR,V,DV,DDV,Z3,Z4,JFLAG) G = DEXP(-V/T)*R*R Gp=G*(1.+F2*((DV/T-4./R)*DV-2.*DDV)) g0=g0+gp*B0**2 g1=g1+gp*BJ**2 100 continue g0=FOURPI*G0*h*(ANGAU*1.D-8)**3 g1=FOURPI*G1*h*(ANGAU*1.D-8)**3 scalej=g1/g0 return end SUBROUTINE BETAJ(R, BETA0, BETA) IMPLICIT double precision (A-H,O-Z) CHARACTER*10 NOTE,LGAS,LABEL dimension BJ(9) COMMON /bpar/ Bpar(9, 4), n COMMON /COMUNIC/ NOTE(2),LGAS,LABEL,FMP,EPSP(2),RMP(2),LVALUE COMMON/ROTIF/ J1, JP1, J2, JP2 DATA DEBYE / 2.541769713d-18 / x = R-6.0 do i = 1, 9 BJ(i) = Bpar(i, 1) / R**N + Bpar(i, 2) * & DEXP( Bpar(i, 3) * x + Bpar(i, 4) * x**2 ) end do Beta = BJ(1) + BJ(2) * J1*(J1+1) + BJ(3) *(J1*(J1+1))**2 + & BJ(4) * J2*(J2+1) + BJ(5) * (J2*(J2+1))**2 + & BJ(6) * JP1*(JP1+1) + BJ(7) * (JP1*(JP1+1))**2 + & BJ(8) * JP2*(JP2+1) + BJ(9) * (JP2*(JP2+1))**2 BETA = Beta * debye*1.0d-6 BETA0= BJ(1)* debye*1.0d-6 RETURN END FUNCTION VA(R) C C double potential H2(v=0) - H2(v=0) initial ---> C H2(v=1) - H2(v=1) as final C H2(v=0) - H2(v=0) initial as the initial and final for 2nd C overtone C ============================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) fACTOR = (4.803250D-10)**2/(0.52916607D-8) *(1.D-7) C POTENTIAL H2(v=0) - H2(v=0) CJB4.28 NOTE='00 - 00 H2' RM = 6.5*0.52917706D0 RMIN1= 6.5D0 E= 0.11163D-3 EPSP= 0.11163D-3 * FACTOR A= 0.37996D7 GAMA= 2.39588D0 ALPHA=14.84169D0 D = 1.02354D0 C6= -12.134D0 C8= -221.26D0 C10= -0.50358D4 C6S = C6 /( E * RMIN1**6 ) C8S = C8 /( E * RMIN1**8 ) C10S = C10 /( E * RMIN1**10 ) X = R/RM V_REP = A * X**GAMA * DEXP(- ALPHA*X) F = 1.D0 IF(X .LE. D ) F = DEXP(- ( D/X-1.)**2) V_DISP = F * (C6S/X**6 + C8S/X**8 + C10S/X**10) FF = V_REP + V_DISP VA= FF * EPSP RETURN END