BLOCK DATA BH2HE01 C ======================================= C Copyright (C) 1992 Aleksandra Borysow C ======================================= C Copyright Notice: C You may copy and distribute this program freely for non-commercial purposes, C provided that you keep the package together and this copyright notice intact. C You may not alter or modify the files; this helps to ensure that all C distributions of H2HE01 are the same. C If you make any modifications, then you must give the files new names, C other than the present. C ======================================= C C Direct all your inquires to the author: e-mail aborysow@nbi.dk C If you decide to use this program, please acknowledge it by c referring to the original paper: c Aleksandra Borysow, c "New model of collision-induced infrared absorption C spectra of H_2-He pairs in the 2-2.5 micrometer region C at temperatures from 20 to 300K: An update", C Icarus, vol. 96, p. 169 -- 175, (1992). c Version running on the SUN c *********************************************************** c Program written by: Aleksandra Borysow c Previously at: Physics Department c Michigan Technological University c Houghton, MI 49931 - 1295 c e-mail: aborysow@phy.mtu.edu, c Currently at: c Niels Bohr Institute for Astronomy, Physics and Geophysics, c University Observatory, c Juliane Maries Vej 30, c DK-2100 Copenhagen, c Denmark c e-mail: aborysow@nbi.dk c *********************************************************** c c Program to model H2-He CIA spectra at low temperatures c for planetary interest (20-300K) c accounts for either equilibrium, or normal H2 c ================================================================ c Program asks interactively for: c TEMPERATURE [K], choose from 20 - 300K c min. frequency (FREQLO) (in cm^[-1]) c max. frequency (FREQMAX); c frequency interval (in cm^{-1}) c c Program computes absorption spectrum at temperatures c between 20 and 300K H2-He fundamental band c J-dependence accounted for (see paper ....****************) c Program uses model lineshapes for J1=J1'=J2=J2':==0 c and rescales their magnitude by the j-dependent c M0 moments : c SEMICLASSICAL M0 moments are used for each J1,J1',J2,J2' c Semiclassical M0 ~~ quantum M0 at temp. >100K (within 1-2% ) c and within much better than 1% at 300K; c NOTE: the inaccuracy is in a term that corresponds by itself c by ~10%max so the approximation is fully justified. c IMPLICIT double precision (A-H,O-Z) COMMON/SPEKTRM/ LIKE, nvib_F,nvib_I, NTYPE, 1 FMAX, FREQLO, DFREQ common /wagi/ w(2) DATA LIKE, nvib_F, nvib_I / 0, 1, 0 / c nvib_F - upper state, nvib_I - lower state DATA NTYPE / 1 / ! choose whenever you want to fix W() DATA W(1), W(2) /1., 3./ ! rotational weights; equilibrium DATA FMAX / 4000.d0/ ! to be kept this way c C NTYPE=0 for equilibrium hydrogen; NTYPE=1 for normal H2 c NF - sets the maximum number of points where the spectrum can be c requested, set as a constant to 601 (fixed now) END PROGRAM H2HE01 c c Program H2HE01 is COPYRIGHT (C) by Aleksandra Borysow 1991 C Program generates a detailed listings of the roto-vibrational spectra c of H2-He in the fundamental band c (0-->1 vibrational transition) at temperatures less than 300K c c Program uses BB lineshapes, as described in the paper c by G. Birnbaum and A. Borysow, Molecular Physics, p.57-68, c vol. 73, 1991. c IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MESSAGE, MARKER CHARACTER*11 labeta common/integ/ rlo, rup, frac, Ninteg common/temp/temp COMMON /PRTSUM/ JRANGE1,MARKER common/prtsum1/ q1, w1(2) common /wagi/ w(2) common/ bb/ labeta, j1, j1p, j2, j2p common/ola/ xxm0 common/sumrul/ g0, g1 COMMON /SPECIT/ LIST common/ specit1/ FREQ(601), ABSCOEF(601), ALFATOT(601) COMMON/SPEKTRM/ LIKE, nvib_F, nvib_I, NTYPE, 1 FMAX, FREQLO, DFREQ common/ read/ go common/om/omega_0 common/ uncor/ xm0, tau1, tau2 common/corr/ omega1, omega2, tauh COMMON/YMAX/YMAXXR, YMAXXL Z(T,a,b,c,d) = A * DEXP( B * dlog(T) + C * (dlog(T))**2 + 1 D * (dlog(T))**3 ) open(unit=10, file='output.h2he', 1 status='unknown',form='formatted') NF = 601 ! with a new Fortran release, NF must be fixed write(*,*) ' Which temperature? ' read(*,*) temp if ( (temp. lt. 18.) .or. (temp.gt.300.)) stop 1299 c Temperatures should be within the range 20-300K C CALCULATE OMEGA_0 C DETERMINING THE PURELY VIBRATIONAL TRANSITION FREQUENCY:(CM^-1) JJ=0 OMEGA_0 = H2ELEV(NVIB_F,JJ) - H2ELEV(NVIB_I,JJ) write(*, 6998) NVIB_I,NVIB_F, OMEGA_0 6998 FORMAT(' VIBRATIONAL TRANSITION FREQUENCY FOR ',I1,' --> ', 1 I1,' IS EQUAL TO',/,F12.5,' CM-1') write(*, 7799) 7799 FORMAT(' The spectrum will be computed from FREQLO', 1 ' to FREQMAX',/,' FREQLO=?') read(*,*) FREQLO write(*, 1256) 1256 format(' FREQMAX=?') read(*,*) freqmax write(*, 1266) 1266 FORMAT(' Frequency STEP in cm^{-1} = ? ') read(*,*) DELTASTEP list = INT (FREQMAX - FREQLO)/DELTASTEP IF(list.GT.601) STOP 963 write(*, 1233) FREQLO, FREQMAX, DELTASTEP, list 1233 FORMAT (' THE SPECTRUM WILL BE COMPUTED FOR FREQUENCIES',/, 1 ' FROM', F12.3,' TO ', F12.3, ' EVERY', F10.1, ' IN ',I4, 1 ' STEPS ') DO 10 I=1,list FREQ(I)=FREQLO+DFLOAT(I-1)*DELTASTEP 10 ALFATOT(I)=0. write(10, 1212) list, freqlo, deltastep 1212 format(' Absorption will be given in ',i4,' points',/, 1 ' From ', f10.2, ' cm-1, every ', f10.2, ' cm-1',/) c =============================================================== c Compute the partition function Q1: CALL PARTFCT (TEMP, Q1, NTYPE, JRANGE1, MARKER, W ) w1(1) = w(1) w1(2) = w(2) ALFA =1./(0.69519*TEMP) hbar = 1.054588757D-27 boltzk = 1.38054d-16 tau0 = hbar/(2.*boltzk*temp) c ***************************************************************** c compute radial distribution function g(R) at temperature TEMP ldist=1 labeta =' call g(R)' rlo = 1.d0 rup = 12.d0 Ninteg = 400 frac = 100.d0 CALL MOMTS (ldist, labeta, RLO,RUP,FRAC,Ninteg,TEMP,XXM0) c ***************************************************************** labeta = '<0|B0001|1>' C leave "labeta" as is, data already for h2he... XM0 = Z(temp,0.26224d-63,-0.75953d0,0.23467d0,-0.94335d-2) XM1 = z(temp,0.49529d-50,-0.68931d0,0.22073d0,-0.93922d-2) g0 = xm0 g1 = xm1 XM2=Z(temp,0.26848d-36,-0.79461d0,0.19498d0,0.20814d-2) c calculate tau1 and tau2 of K0 from the three above moments DELT = (TAU0*XM0)**2-4.*(XM1*XM1/XM0+XM1/TAU0-XM2)*TAU0*TAU0*XM0 TAU1=(-DSQRT(DELT)-TAU0*xm1)/(2.*(XM1*XM1/XM0+XM1/TAU0-XM2)) TAU1=DSQRT(TAU1) TAU2=TAU0*TAU1*xm0/(xm1*TAU1*TAU1-TAU0*xm0) omega1 = Z(temp,0.29314d2,-0.80060d0,0.15804d0,-0.37519d-2) omega2 = Z(temp,0.16982d2,-0.25270d0,0.15023d-1,0.63805d-2) tauH = Z(temp,0.43044d-12,-0.49836d0,-0.48001d-2,0.11347d-2) C write(10,*) labeta,xm0,xm1,xm2,tau1,tau2,omega1,omega2, tauH ibgama = 0 ! means, X==K0 freqcut = 4000.d0 call range(ibgama) CALL ADDSPEC(TEMP,LIKE, 0, 0, 0, 1, freqcut,nvib_F,nvib_I, 2 ibgama, XM0, tau1,tau2, omega1, omega2, tauH ) DO 920 I=1,list 920 ALFATOT(I) = ABSCOEF(I) c ***************************************************** c 0445 Labeta='<0|B0445|1>' c data for H2-He OK here; c Moments fitted 00|0445|01 XM0 = Z(temp,0.88266d-66,-0.61877d0,0.18463d0,-0.39652d-2) XM1 = Z(temp,0.27524d-52,-0.68438d0,0.20160D0,-0.46843d-2) g0 = xm0 g1 = xm1 XM2 = Z(temp, 0.15703d-38,-0.67076d0,0.15878d0,0.67472d-2) c calculate tau1 and tau2 of K0 from the three above moments DELT = (TAU0*XM0)**2-4.*(XM1*XM1/XM0+XM1/TAU0-XM2)*TAU0*TAU0*XM0 TAU1=(-DSQRT(DELT)-TAU0*xm1)/(2.*(XM1*XM1/XM0+XM1/TAU0-XM2)) TAU1=DSQRT(TAU1) c parameters for correction term: omega1 = z(temp,0.30144d+02,-0.81980d0,0.16656d0,-0.42591d-2) omega2 = z(temp,0.17688d+2,-0.21369d0,0.44073d-2,0.73000d-2) tauH = z(temp,0.36923d-12,-0.51052d0,0.23942d-2,-0.12392d-3) C write(10,*)labeta,xm0,xm1,xm2,tau1,tau2,omega1,omega2, tauH ibgama = 1 ! means, ==BC call range(0) CALL ADDSPEC(TEMP,LIKE, 0, 4, 4, 5, FMAX, nvib_F,nvib_I, 2 ibgama, XM0, tau1,tau2, omega1, omega2, tauH ) DO 27 I=1,list 27 ALFATOT(I)= ALFATOT(I)+ ABSCOEF(I) c * * * * * * * * * * * * * ** * * * * * * * * * * * c 0221 Labeta='<0|B0221|1>' c data for H2-He OK here c Moments fitted 00|0221|01 XM0 = Z(temp,0.22683d-64,-0.64631d+00,0.19421d+00,-0.44912d-02) XM1 = Z(temp,0.50240d-51,-0.64046d0,0.19290d0,-0.43866d-2) g0 = xm0 g1 = xm1 XM2 = Z(temp,0.29112d-37,-0.68647d0,0.15610d0,0.73255d-2) c calculate tau1 and tau2 of K0 from the three above moments DELT = (TAU0*XM0)**2-4.*(XM1*XM1/XM0+XM1/TAU0-XM2)*TAU0*TAU0*XM0 TAU1=(-DSQRT(DELT)-TAU0*xm1)/(2.*(XM1*XM1/XM0+XM1/TAU0-XM2)) TAU1=DSQRT(TAU1) c parameters for correction term: omega1 = z(temp,0.31140d2,-0.77014d0,0.14912d0,-0.27394d-2) omega2 = z(temp,0.17672d2,-0.19539d0,-0.47614d-3,0.77027d-2) tauH = z(temp,0.40303d-12,-0.50558d0,-0.35074d-3,0.34315d-3) C write(10,*)labeta,xm0,xm1,xm2,tau1,tau2,omega1,omega2,tauH ibgama = 0 ! means, ==K0 call range(0) CALL ADDSPEC(TEMP,LIKE, 0, 2, 2, 1, FMAX, nvib_F,nvib_I, 2 ibgama, XM0, tau1,tau2, omega1, omega2, tauH ) DO 24 I=1,list 24 ALFATOT(I)= ALFATOT(I)+ ABSCOEF(I) c * * * * * * * * * * * * * ** * * * * * * * * * * * labeta='<0|B0223|1>' c data for H2-He OK here C Moments fitted 00|0223|01 XM0 = Z(temp,0.36021d-64,-0.61751d-1,0.34871d-1,0.43329d-2) XM1 = Z(temp,0.46654d-51,-0.26894d0,0.98138d-1,0.52341d-3) XM2=Z(temp,0.25990d-37,-0.67277d0,0.17143d0,0.37406d-2) g0 = xm0 g1 = xm1 c find tau1 and tau2 of BC... TTA=XM0*TAU0/XM1 XXX = (XM2*TTA-XM0*(1.d0 + TAU0**2/TTA))/(XM0*(TAU0/TTA)**2) TAU1=DSQRT(XXX) TAU2=TTA/TAU1 omega1 = Z(temp,0.16545d+2,-0.10967d1,0.26608d0,-0.11929d-1) omega2 = Z(temp,0.18147d2,-0.42255d0,0.54097d-1,0.40845d-2) tauH = Z(temp,0.44703d-12,-0.48130d0,-0.32231d-2,0.13402d-3) C write(10,*) labeta,xm0,xm1,xm2,tau1,tau2,omega1,omega2,tauH ibgama = 1 ! means, ==BC call range(1) CALL ADDSPEC(TEMP,LIKE, 0, 2, 2, 3, FMAX, nvib_F,nvib_I, 2 ibgama, XM0, tau1,tau2, omega1, omega2, tauH ) DO 324 I=1,list 324 ALFATOT(I)= ALFATOT(I)+ ABSCOEF(I) c * * * * * * * * * * * ** * * * * * * * * * * * MESSAGE = 'GRANDTOTAL' CALL ALFAMOM (TEMP,MESSAGE,10,LIKE,-1,-1,-1,-1) 11 continue stop end subroutine addspec( temp, $ like,lambda1,lambda2,lambda,lvalue, cut, nvib_F, nvib_I, 2 ibgama, G0, tau1,tau2, omega11, omega22, tauH ) C c Subroutine ADDSPEC designed for (molecule-atom) c collisions; will assume lambda==0, lambda2.ne.0 C This subroutine generates a detailed listing of CIA ALFA(OMEGA). C FREQCUT (or, CUT) serves two purposes: cut-off for moment integrals, and C to skip interpolation of G(OMEGA) for OMEGA > FREQCUT. c IBGAMA=0 means uses K0, IBGAMA=1 (uses BC) C LIKE=1 for like systems (as H2-H2) C IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MESSAGE, MARKER CHARACTER*11 labeta COMMON /PRTSUM/ JRANGE1,MARKER common/sumrul/ gg0, gg1 common/prtsum1/ q1, w1(2) COMMON /SPECIT/ LIST common/specit1/ FREQ(601),ABSCOEF(601), ALFATOT(601) common /BB/ labeta, J1, Jp1, J2, Jp2 common/ola/ XM00 common/gam/ gamma(3) common/om/omeg_0 common/ integ / rlo, rup, frac, ninteg COMMON/YMAX/YMAXXR, YMAXXL DATA SCALEF/1.D80/, BOLTZWN/ 0.6950304D0/ DATA HBAR,PI,CLIGHT/1.054588757D-27, 3.141592653589797D0, 1 2.9979250D10/ DATA Boltzk / 1.380658d-16 / if(lambda1.ne.0) stop 999 ! should be 0 nf = 601 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) WRITE (10,270) LAMBDA1,LAMBDA2,LAMBDA,LVALUE,JRANGE1,MARKER, 1 W1(1),W1(2) 270 FORMAT (//,' LAMBDA1,LAMBDA2, LAMBDA,LVALUE=',2I3,2X,2I3,8X, 1 'RANGE1=',I2,/, ' OUTPUT ADDSPEC',10X,A10,'HYDROGEN',2F10.2,/, 2 1X,45(1H=) ) C Molecule interacting with an atom: XM00 = G0 ! pure quantum result c Checking, M0 semiclassical, how different from quantum ??? ldist = 0 j1=0 jp1=0 j2=0 jp2=0 CALL MOMTS (ldist, labeta, RLO,RUP,FRAC,Ninteg,TEMP, XM00p) write(10, 7551) LABETA, XM00, XM00p 7551 format(/,' For ', a11, ' and all J=0, quantum M0=',e12.5,/, 1 12x, ' ... and semiclassical M0 =', e12.5 ,/) JPLUSL=JRANGE1+LAMBDA2 DO 150 I2=1,JRANGE1 J2=I2-1 DO 150 IP2=1,JPLUSL JP2=IP2-1 CG2S=CLEBSQR(J2,LAMBDA2,JP2) IF (CG2S.eq.0.) go to 150 P2=H2POPL(nvib_I,J2,TEMP,W1)/Q1 OMEGA2=H2ELEV(nvib_f,JP2)-H2ELEV(NVIB_i,J2) c H2 does undergoes roto-vibrational transition ldist = 0 CALL MOMTS (ldist, labeta, RLO,RUP,FRAC,Ninteg,TEMP, XM0j ) c Below a correction made for different M0(J2,J2'): c write(10,*) 'J2,Jp2,P2; M0j/M00:', J2,Jp2,P2, xm0j/xm00 FAC =CALIB *P2 *CG2S* (Xm0j/Xm00) DO 140 I=1,LIST FRQ=FREQ(I)-OMEGA2 c IF (DABS(FRQ)-FREQCUT) 130,140,140 130 if (frq.gt.ymaxxR) go to 140 if (frq.lt.ymaxxL ) go to 140 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) yy = g0/2. * (bgama_h(frq-(omega11+omega22),tauH,temp) - 1 bgama_h(frq-(omega22-omega11), tauH, temp)) gg = xx + yy ABSCOEF(I) = ABSCOEF(I) + gg * WKF 140 CONTINUE 150 CONTINUE CALL ALFAMOM (TEMP,MESSAGE,10,LIKE,LAMBDA1,LAMBDA2,LAMBDA,LVALUE) C c checking the sum rules: c I need M0 and M1 of the translational lineshape: G0 and G1 (unperturbed) vib = omeg_0 * 2.*pi*clight b_rot = (h2elev(0,1) - h2elev(0,0))/2.d0 beta = 1./(boltzk * temp) GG1 = GG1 + 1 GG0/2.*(omega11+omega22) - GG0/2.*(omega22-omega11) XW1 = GG1 + 2.*pi*clight*b_rot* 1 dfloat(lambda2*(lambda2+1)) * GG0 + GG0 * vib alpha1 = 2.*2.*pi**2/(3.*hbar*clight)*XW1 gama1 = 2.* pi*pi/(3.*clight*boltzk*temp)*GG0 write(10, 2222) alpha1, gama1 2222 format(' Moments from the sum rules:',/,' alfa1=', e12.5,/, 1 ' gama1=', e12.5,/) C RETURN END SUBROUTINE PARTFCT (TEMP, Q, NTYPE, JRANGE, MARKER, W) C C NTYPE = 0 for equilibrium hydrogen C NTYPE = 1 redefines the weights W for "normal" hydrogen with a C ratio of ortho:para H2 of 3:1. C IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MARKER DIMENSION W(2) DATA BOLTZWN /.6950304D0/ C 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) c 1 h2popl(4,j,temp,w) + h2popl(5,j,temp,w) Q = Q + DQ J = J + 1 if ( (w(2).eq.0.).and.(mod(j,2).eq.1)) go to 30 ! odd J: do not compare IF (DQ.GT.Q/990.) GO TO 30 J=-1 S=0. IF (NTYPE.EQ.1) GO TO 60 40 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) c 1 h2popl(4,j,temp,w) + h2popl(5,j,temp,w) S= S + DDD/Q IF (S.LE.0.999D0) GO TO 40 JRANGE=J+2 MARKER='EQUILIBRM ' GO TO 90 C C "NORMAL" HYDROGEN REQUIRES REDEFINITION OF W(2): C 60 SEV=0. 70 J=J+1 DS=H2POPL(0, J, TEMP, W)+ h2popl(1,j,temp,w) + 1 h2popl(2,j,temp,w) + h2popl(3,j,temp,w) c 1 h2popl(4,j,temp,w) + h2popl(5,j,temp,w) S=S+DS SEV=SEV+DS J=J+1 S=S+H2POPL(0, J, TEMP, W)+ h2popl(1,j,temp,w) + 1 h2popl(2,j,temp,w) + h2popl(3,j,temp,w) c 1 h2popl(4,j,temp,w) + h2popl(5,j,temp,w) IF (S.LE.0.999*Q) GO TO 70 JRANGE=J+2 SODD=S-SEV W(2)=W(2)*(3.*SEV)/SODD C C DEFINITION OF "NORMAL" HYDROGEN: 3*S(EVEN) = S(ODD) C Q=4.*SEV MARKER='NORMAL (!)' C 90 continue WRITE (10,310) marker, TEMP,Q,W(1),W(2) 310 FORMAT (1x,a10, ' HYDROGEN',/, 1 ' THE ROTATION PARTITION SUM OF H2 ',/, 1 ' AT TEMPERATURE =', F8.2, 8H K IS Q=, E12.4,/, 1 ' Weights J_even/odd:', 2F10.3,/ ) 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 SUBROUTINE ALFAMOM(TEMP,MESSAGE,K,LIKE, 1 LAMBDA1,LAMBDA2,LAMBDA,LVALUE) C C Writes out the ALFA (absorption) array C IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MESSAGE COMMON /SPECIT/ LIST common/specit1/ FREQ(601),Abscoef(601),Alfatot(601) DIMENSION XI(1), aa(601), abs2(601), absp(601), fx(1) DATA EPS, XI(1) /1.D-4, 3000.D0/, PI /3.141592653589797D0/ DATA CLIGHT,CLOSCHM /2.9979250D10, 2.68675484D19/ DATA HBAR,BK / 1.054588757D-27, 1.38066244D-16/ do 1555 i=1, list absp(i) = 0.d0 abs2(i) = 0.d0 1555 continue IF(MESSAGE.EQ. 'GRANDTOTAL' ) go to 11 if(message.eq. 'NON-COMUL.') go to 12 stop 987 11 do 22 i=1, list 22 aa(i)=alfatot(i) go to 13 12 do 23 i=1, list 23 aa(i)=abscoef(i) 13 continue TWOPIC= 2.*PI*CLIGHT CALIB = TWOPIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*CLOSCHM**2 1 /DFLOAT(1+LIKE) WRITE (K,330) FREQ(1),FREQ(LIST),FREQ(2)-FREQ(1),TEMP, 1 LAMBDA1,LAMBDA2,LAMBDA,LVALUE,MESSAGE 330 FORMAT (/,' ABSORPTION FROM',F10.1,' CM-1 TO',F7.1,' CM-1.', 1 ' every',F8.2, 10x, /, ' Temperature:', f8.2, /, 1 2X,4I2,5X,' IN UNITS OF CM-1 AMAGAT-2.',4X,A10/) WRITE (K,340) ( AA(I),I=1,LIST) 340 FORMAT ( 90(10E13.5,/)) IF(MESSAGE.EQ.'GRANDTOTAL' ) 1 open(unit=11,file='plot',status='unknown', 1 form='formatted') IF(MESSAGE.EQ. 'GRANDTOTAL' ) 1 write(11,2999) (freq(i), aa(i), i=1, list) 2999 format( f10.2, e12.5) IF(MESSAGE.EQ. 'GRANDTOTAL' ) close(11) c testing and computing the integrals (moments) call spline(list, 1,2,eps,freq,aa,xi,fx,alfa1,nr,abs2) alfa1 = alfa1 * twopic/closchm**2 do 250 i=2, list e2 = dexp(hbar*twopic*freq(i)/(bk*temp)) div = freq(i)*(e2-1.)/(e2+1.) absp(i) = aa(i)/div 250 continue fr = (freq(2)/freq(3))**2 absp(1) = ( absp(2) - fr * absp(3))/(1.-fr) call spline(list,1,2,eps,freq,absp,xi,fx,gama1,nr,abs2) gama1 = gama1 * hbar/(2.*bk*temp*closchm**2) c here I have alfa1 and gama1 write(10, 1111) message,freq(1),freq(list),alfa1,gama1 1111 format(1x,a10,3x, 1 ' Moments alfa1 and gama1 from the integration',/, 1 ' from ', f10.2, 'to', f10.2,/, ' alfa1:', e12.5,/, 1 ' gama1:',e12.5,//) 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 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 SUBROUTINE MOMTS (IFIRST, labeta, RLO,RUP,FRAC,N,TEMP, G0) c c If IFIRST = 1 then do only preliminary set-up; set ifirst=0 c when calling moments, c G0 is the zeroth semiclassical moment C C COMPUTES 0-TH, SEMICLASSICAL TRANSLATIONAL SPECTRAL MOMENT C WORKS WITH "DOUBLE" POTENTIAL FUNCTIONS IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*10 LGAS, LABEL, NOTE(2) COMMON /COMUNIC/ NOTE(2),LGAS,LABEL,FFMP,EPSP(2),RMN(2), 1 LVALUE COMMON/EPOT/EPS(2),RM(2),FM,RMIN,RMAX,HHH DATA ANGAU / 0.529917706 / C RLO,RUP - LOWER AND UPPER LIMIT OF INTEGRATION IN ANG C N - NUMBER OF STEPS, FRAC REQUIRED BY DERIVAS, USUALLY SET=100. C FMP - REDUCED MASS IN HYDROGEN MASS, LVALUE = L C LGAS - 10 CHARACTER DESCRIPTION OF SYSTEM if (ifirst.eq.1) go to 222 go to 333 222 BBB=BETACOM(0.D0) ! call BETA_JJ.for VVV=VACOMUN(0.D0) ! call VA(R,k) RM(1)=RMN(1)/ANGAU EPS(1)=EPSP(1) RM(2)=RMN(2)/ANGAU EPS(2)=EPSP(2) 10 FORMAT(3F10.5,I5) 11 FORMAT(F10.5) RLO1=RLO/ANGAU RUP1=RUP/ANGAU H=(RUP1-RLO1)/DFLOAT(N) T= TEMP * 1.380662D-23 FMP=FFMP/2.d0 FMPI=2.*FMP*1.67265D-27 return 333 CALL MOMENTS(RLO1,H,N,T,LVALUE,FRAC,FMPI,Gamma, gammap) G0=GAMMA RETURN END SUBROUTINE MOMENTS(RLO,H,N,T,LVALUE,FRAC,FMP,G0, g0p) C COMPUTES ZEROTH MOMENT (G0) FOR ANY L VALUE c G0p is the confirmation of the accuracy; lower order diff. IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON/ORDER/Y1P,Y2P,Y3P,Y4P EXTERNAL VTILDE,BETA DATA ANGAU, FOURPI, HBAR/ 0.52917706, 12.566370614359, 1 1.05458876D-34/ SP=0. SQ = 0. R=RLO-H F2=HBAR*HBAR/(12.*FMP*T*T*ANGAU**2*1.D-20) DO 10 I=1,N R=R+H DR=R/FRAC B = beta(R) CALL DERIVAS(VTILDE,R,DR,V,DV,DDV,Z3,Z4,JFLAG) DVP=Y1P DDVP=Y2P G = DEXP(-V/T)*R*R G2=G*(1.+F2*((DV/T-4./R)*DV-2.*DDV)) G2P=G*(1.+F2*((DVP/T-4./R)*DVP-2.*DDVP)) SP=SP+B*B*G2 SQ=SQ+B*B*G2P 10 CONTINUE G0 = FOURPI*SP*H*(ANGAU*1.D-8)**3 G0p = FOURPI*SQ*H*(ANGAU*1.D-8)**3 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) COMMON/ORDER/Y1P,Y2P,Y3P,Y4P 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 FOLLOWING ARE LOWER ORDER DERIVATIVES FOR TESTING Y1P=(B1-8.*C1)/(12.*H) Y2P=(16.*C-B-30.*Y)/(12.*HSQ) Y3P=(2.*C1-B1)/(2.*H*HSQ) Y4P=(B-4.*C+6.*Y)/(HSQ*HSQ) 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 VTILDE(R) IMPLICIT DOUBLE PRECISION(A-H,O-Z) K=0 VTILDE = VA(R,K) RETURN END FUNCTION BETA(R) IMPLICIT REAL*8 (A-H,O-Z) character*11 labeta CHARACTER*10 NOTE(2),LGAS,LABEL COMMON /BB/ LABETA, J1, J1p, J2, J2p common/ ola/ XM0 COMMON /COMUNIC/ NOTE(2),LGAS,LABEL,FMP,EPSP(2),RMP(2),LVALUE DATA DEBYE / 2.541769713d-18 / DATA E0, A0 / 4.803250D-10, 0.52917706D-8 / RP = R -5.70D0 J = J2 JP = J2p if (labeta .eq. '<0|B0001|1>') go to 10 if (labeta .eq. '<0|B0223|1>') go to 30 if (labeta .eq. '<0|B0221|1>') go to 60 if (labeta .eq. '<0|B0445|1>') go to 90 write(*, 555) labeta 555 format(' Unpredicted LABEL in BETA:', a11) stop 999 ! then, what term is that??? 10 continue ! 0001 here LAMDA = 0 LVALUE = 1 LABEL = '00|0001|01' C =(24.7 +0.216 *dfloat(J *(J+1)) -0.247 *dfloat(JP *(JP+1))) D0 = -0.000751 D1 = -0.00000868 *dfloat(J *(J+1)) D2 = 0.00000798 * dfloat(JP *(JP+1)) N = 7 REST = D0 *DEXP((-1.538 -0.051 *RP) *RP) 1 +D1 *DEXP((-1.675 -0.056 *RP) *RP) 2 +D2 *DEXP((-1.683 -0.049 *RP) *RP) BETA = C/R**N + REST go to 100 30 continue ! 0223 LAMDA = 2 LVALUE =3 LABEL = '00|0223|01' N =4 C = (-.205 -0.00393 * dfloat(J *(J+1)) 1 +0.00385 * dfloat(JP *(JP+1))) D0 = -0.000119 D1 = -0.00000101 * dfloat(J *(J+1)) D2 = 0.000000666 * dfloat(JP *(JP+1)) REST = D0 *DEXP((-1.583 -0.053 *RP) *RP) 1 +D1 *DEXP((-1.615 -0.057 *RP) *RP) 2 +D2 *DEXP((-1.729 -0.084 *RP) *RP) BETA = C/R**N + REST go to 100 60 continue ! 0221 here LAMDA = 2 LVALUE = 1 LABEL = '00|0221|01' N = 7 C = (-2.86 -0.00494 * dfloat(J *(J+1)) 1 -0.0548 * dfloat(JP *(JP+1))) D0 = -0.000170 D1 = -0.00000145 * dfloat(J *(J+1)) D2 = 0.000000806 * dfloat(JP *(JP+1)) rest = D0 *DEXP((-1.640 -0.018 *RP) *RP) 1 +D1 *DEXP((-1.702 -0.022 *RP) *RP) 2 +D2 *DEXP((-1.918 -0.089 *RP) *RP) BETA = C/R**N + REST go to 100 90 continue ! 0445 here LAMDA = 4 LVALUE = 5 LABEL = '00|0445|01' N = 6 C = (-0.511 -0.00540 * dfloat(J *(J+1)) 1 + 0.00502 * dfloat(JP *(JP+1))) D0 = -0.0000218 D1 = -0.000000132 * dfloat(J *(J+1)) D2 = 0.0000000501 * dfloat(JP *(JP+1)) rest = D0 *DEXP((-1.873 -0.097 *RP) *RP) 1 +D1 *DEXP((-2.023 -0.138 *RP) *RP) 2 +D2 *DEXP((-2.309 -0.223 *RP) *RP) BETA = C/R**N + REST 100 continue beta = beta * debye return ENTRY BETACOM ! same for all terms BETA = 0.D0 RETURN END FUNCTION VA(R,K) c C H2-He ab initio potential by Meyer, Hariharan and Kutzelnigg, C J.C.P.73(1980)1880, supplemented by Schaefer and Koehler, C Physica 129A(1985)469 (isotropic component, vib. averages) C DIFFERENT initial and final state potentials: C Different vibrational averages, like for v=0,1) c Analytical representation: HFD model fitted to numerical data c (27-Nov-89) C ==================================================== c c v=0 (initial) --> v'=1 (final) c ===================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*10 NOTE(2), LGAS, LABEL COMMON /EPOT/ EPS(2),RV(2),FM,RMIN,RMAX,HQ COMMON /COMUNIC/ NOTE(2),LGAS,LABEL,FMP,EPSP(2),RM(2),LVALUE COMMON/VVP/IVVP DIMENSION A(2), GAMA(2), ALPHA(2), D(2), RMIN1(2), E(2), 1 C6(2), C8(2), C6S(2), C8S(2) IVVP=K K1=K+1 X = R/RV(K1) V_REP = A(K1) * X**GAMA(K1) * DEXP(- ALPHA(K1)*X) F = 1.D0 IF(X .LE. D(K1) ) F = DEXP(- ( D(K1)/X-1.)**2) V_DISP = F * (C6S(K1)/X**6 + C8S(K1)/X**8 ) FF = V_REP + V_DISP VA = FF * EPS(K1) RETURN ENTRY VACOMUN C Start-up procedure: INITIAL STATE IS FIRST: c Potential H2(v=0)-He C FACTOR = (4.803250D-10)**2/(0.52916607D-8) *(1.D-7) NOTE(1)='MEYER - 00' LGAS='H2-He mix.' FMP=2.6603 C write(6, 777) lgas, fmp 777 format(' In VACOMUN: Gas=', a10, ' FMP=', f10.2, 'a.u.') RM(1) = 6.4 * 0.52917706 RMIN1(1) = 6.4d0 EPSP(1)= 0.42957d-4 * factor E(1) = 0.42957d-4 A(1)= 0.95289d6 GAMA(1)= 0.10029d1 ALPHA(1)= 0.13966d2 D(1) = 0.13724d1 C6(1)= -0.3587d1 C8(1)= -0.1038d3 C C FINAL STATE: C POTENTIAL H2(v=1) - He C NOTE(2)='MEYER - 11' RM(2) = 6.4 * 0.52917706 RMIN1(2) = 6.4d0 EPSP(2)= 0.41879d-4 * factor E(2) = 0.41879d-4 A(2) = 0.10923d7 GAMA(2)= 0.10759d1 ALPHA(2)= 0.14026d2 D(2) = 0.1462d1 C6(2)= -0.3656d1 C8(2)= -0.12426d3 DO 10 I=1,2 C6S(I) = C6(I) /( E(I) * RMIN1(I)**6 ) C8S(I) = C8(I) /( E(I) * RMIN1(I)**8 ) 10 continue c VA=0. RETURN END SUBROUTINE SPLINE (L,M,K,EPS,X,Y,T,SS,SI,NR,S2) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(L), Y(L), T(M), SS(M), S2(L) DELY(I)=(Y(I+1)-Y(I))/(X(I+1)-X(I)) B(I)=(X(I)-X(I-1))*0.5/(X(I+1)-X(I-1)) C(I)=3.*(DELY(I)-DELY(I-1))/(X(I+1)-X(I-1)) N=L N1=N-1 DO 10 I=2,N1 10 S2(I)=C(I)/1.5 OMEGA=1.0717968 IC=0 C IF (NR.EQ.-1) GOTO 20 S2(N)=0. S2(1)=0. 20 ETA=0. IC=IC+1 SM=DABS(S2(1)) DO 30 I=2,N IF (DABS(S2(I)).GT.SM) SM=DABS(S2(I)) 30 CONTINUE EPSI=EPS*SM DO 50 I=2,N1 W=(C(I)-B(I)*S2(I-1)-(0.5-B(I))*S2(I+1)-S2(I))*OMEGA IF (DABS(W)-ETA) 50,50,40 40 ETA=DABS(W) 50 S2(I)=S2(I)+W IF (ETA-EPSI) 60,20,20 ENTRY IXPOLAT C C THIS ENTRY USEFUL WHEN ITERATION PREVIOUSLY COMPLETED C N=L N1=N-1 IC=-1 60 IF (K.EQ.2) GO TO 260 NR=0 DO 250 J=1,M I=1 IF (T(J)-X(1)) 110,210,80 80 IF (T(J)-X(N)) 100,190,150 90 IF (T(J)-X(I)) 200,210,100 100 I=I+1 GO TO 90 110 NR=NR+1 HT1=T(J)-X(1) HT2=T(J)-X(2) YP1=DELY(1)+(X(1)-X(2))*(2.*S2(1)+S2(2))/6. IF (K) 140,130,120 120 SS(J)=YP1+HT1*S2(1) GO TO 250 130 SS(J)=Y(1)+YP1*HT1+S2(1)*HT1*HT1/2. GO TO 250 140 SS(J)=S2(I) GO TO 250 150 HT2=T(J)-X(N) HT1=T(J)-X(N1) NR=NR+1 YPN=DELY(N1)+(X(N)-X(N1))*(S2(N1)+2.*S2(N))/6. IF (K) 180,170,160 160 SS(J)=YPN+HT2*S2(N) GO TO 250 170 SS(J)=Y(N)+YPN*HT2+S2(N)*HT2*HT2/2. GO TO 250 180 SS(J)=S2(N) GO TO 250 190 I=N 200 I=I-1 210 HT1=T(J)-X(I) HT2=T(J)-X(I+1) PROD=HT1*HT2 S3=(S2(I+1)-S2(I))/(X(I+1)-X(I)) SS2=S2(I)+HT1*S3 DELSQS=(S2(I)+S2(I+1)+SS2)/6. IF (K) 240,220,230 220 SS(J)=Y(I)+HT1*DELY(I)+PROD*DELSQS GO TO 250 230 SS(J)=DELY(I)+(HT1+HT2)*DELSQS+PROD*S3/6. GO TO 250 240 SS(J)=SS2 250 CONTINUE 260 SI=0. DO 270 I=1,N1 H=X(I+1)-X(I) 270 SI=SI+0.5*H*(Y(I)+Y(I+1))-H**3*(S2(I)+S2(I+1))/24. IF (K.EQ.2) NR=IC RETURN C END FUNCTION CLEBSQR (L,LAMBDA,LP) C SQUARE OF CLEBSCH-GORDAN COEFFICIENT (L,LAMBDA,0,0;LP,0) C FOR INTEGER ARGUMENTS ONLY C IMPLICIT REAL*8 (A-H,O-Z) FC=DFLOAT(2*LP+1) GO TO 10 C ENTRY THREEJ2 C 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 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; unique for "H" function 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.883651568E11,.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 RANGE(IBGAMA) implicit double precision (a-h,o-z) COMMON/YMAX/ymaxxr, ymaxxl COMMON/UNCOR/ G0, TAU1, TAU2 COMMON/CORR/ omega1, omega2, tauH common/temp/temp c choosing YMAXX here: if(ibgama.eq.1) gmax =g0 * bgama1(0.,tau1,tau2,temp) if(ibgama.eq.0) gmax =g0 * bgama0(0.,tau1,tau2,temp) STEP = 100.D0 XS = 100.D0 omegapl = omega1 + omega2 omegamn = omega2 - omega1 4441 if(ibgama.eq.1) gx = g0 * bgama1(xs,tau1,tau2,temp) if(ibgama.eq.0) gx =g0 * bgama0(xs,tau1,tau2,temp) gy = g0/2. * ( bgama_h(xs- omegapl,tauH,temp) - 1 bgama_h(xs - omegamn, tauH, temp) ) GLOWR = GX + GY XS = XS + STEP if(xs.gt.5000. ) go to 6565 IF(GLOWR. GT. GMAX/5000.D0) GO TO 4441 6565 YMAXXR = XS-STEP if(ymaxxr.lt.1000.d0) ymaxxr = 2000. XS = -100.D0 if(ibgama.eq.1) gx =g0 * bgama1(0.,tau1,tau2,temp) if(ibgama.eq.0) gx =g0 * bgama0(0.,tau1,tau2,temp) 6441 XS = XS - STEP if(xs.le.-5000.) go to 7766 if(ibgama.eq.1) gx =g0 * bgama1(xs,tau1,tau2,temp) if(ibgama.eq.0) gx =g0 * bgama0(xs,tau1,tau2,temp) gy = g0/2. * ( bgama_h(xs- omegapl,tauH,temp) - 1 bgama_h(xs - omegamn, tauH, temp) ) if((gx+gy).le.0.d0) go to 7766 GLOWL = GX + GY c write(6, 122) xs, gx, gy, glowl IF(GLOWL. GT. GMAX/5000.D0) GO TO 6441 7766 YMAXXL = XS + STEP if(ymaxxL.gt.-1000.d0) ymaxxL = -2000.d0 C write (6, 133) ymaxxR, ymaxxL C 133 format(/,' Selected ymaxxR, ymaxxL', 2f12.3,/) return end