C N. R. BADNELL NRB v6.0.1 05/11/21 C PROGRAM TPSMX C C----------------------------------------------------------------------- C C TEST-DRIVER/WRAPPER FOR L-CHANGING COLLISIONS (DELTA-N=0~DELTA-E). C C FOR DIPOLE TRANSITIONS, COMPARES THE ORIGINAL PENGELLY & SEATON C RESULTS (PS64) MNRAS 127, 165 (1964) WITH THE MODIFIED PENGELLY & C SEATON RESULTS (PSM20) OF BADNELL ET AL MNRAS XXX, XXX (2021) WHICH C IMPROVES-UPON THE ORIGINAL MODIFIED PENGELLY SEATON (PSM17) METHOD C INTRODUCED BY GUZMAN ET AL MNRAS 464, 312 (2017). C C INTERACTIVE VERSION ASKS USER FOR TEMPERATURE, DENSITY, N-RANGE, C L-RANGE AND MULTIPOLE (=DELTA-L= 1 OR 0 FOR RESOLVED OR UNRESOLVED.) C C HARDWIRED IS THE SET-UP OF PROJECTILE & TARGET MASSES & CHARGES. C DEFAULT IS PROTONS ON HYDROGEN - OTHER EXAMPLES ARE GIVEN BELOW. C C OUTPUT IS MAXWELLIAN RATE COEFFICIENTS PSM20 & PS64 (DIPOLE ONLY). C THE PSM20 RESULTS ARE JUST LABELLED PSM (AS PSM17 IS ALSO POSSIBLE.) C ALL RESULTS ARE "ANALYTIC", HENCE INTERNAL "X" LABEL. C C BASED-ON TLCHNG v6.0 18/08/21 C C----------------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) C LOGICAL bout,bout7,bsum,brslv C PARAMETER (MXNVAL=10000) !MAX PRINCIPAL QUANTUM NO. ALLOW C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) PARAMETER (DTEN=10.0D0) PARAMETER (D1P10=1.D10) PARAMETER (D1M20=1.D-20) C PARAMETER (IZERO=0) PARAMETER (IONE=1) C PARAMETER (TKAY=1.5789D5) !KELVIN TO RYDBERGS PARAMETER (BOHRA=5.29177D-9) !CM TO A.U. PARAMETER (EMASS=5.485798D-4) !ELECTRON MASS A.M.U C !SAMPLE PROJECTILES PARAMETER (PROTON= 1.00727647) PARAMETER (DEUTERON=2.01355323) PARAMETER (ALPHAP= 4.00150618) C !SAMPLE TARGETS PARAMETER (HYDROGEN= 1.007825) PARAMETER (DEUTERIUM=2.014102) PARAMETER (TRITIUM= 3.016049) PARAMETER (HELIUM3=3.016029) PARAMETER (HELIUM4=4.002603) PARAMETER (HELIUM=HELIUM4) PARAMETER (CARBON=12.0) C PARAMETER (T2NEBIG=1.D+25) !MAX ALLOWED VALUE OF T**2/Ne PARAMETER (TOLT2NE=1.D-15) !MIN ALLOWED VALUE OF T**2/Ne c data ntest0/20/ !tolerance parameter for largest "safe" n for T/N C PI=ACOS(-DONE) !F2003 PARAMETER C C*********************************************************************** C write(*,*) write(*,*) & '*** L-CHANGING COLLISION RATE COEFFICIENTS (v6.0.1) ***' write(*,*) write(*,*) & '*** MODIFIED PENGELLY & SEATON (PSM20) IN USE (R/R_1) ***' write(*,*) C C I/O SET-UP C (SET UNIT NUMBER IOUT NEGATIVE TO SUPPRESS WRITES THERE) C iout= 0 !=0 rate coeffs to screen bout=iout.ge.0 if(iout.gt.0)open(iout,file='tpsmx.out') !=6 to outfile c if(iout.gt.0)write(iout, &"(/'# *** L-CHANGING COLLISION RATE COEFFICIENTS (v6.0) ***'/)") if(iout.gt.0)write(iout, &"('# *** MODIFIED PENGELLY & SEATON (PSM20) IN USE (R/R_1) ***'/)" &) c iout7= 7 !.lt. to suppress writes bout7=iout7.gt.0 if(bout7)open(iout7,file='tpsmx.sum') !summd/unresolvd rate coeffs c if(bout7)write(iout7, &"(/'# *** L-CHANGING COLLISION RATE COEFFICIENTS (v6.0) ***'/)") if(bout7)write(iout7, &"('# *** MODIFIED PENGELLY & SEATON (PSM20) IN USE (R/R_1) ***'/) &") C C*********************************************************************** C************************** BEGIN USER SET-UP ************************** C*********************************************************************** C C SET-UP PROJECTILE AND TARGET CHARGES (ZP,ZT) AND MASSES (PMASS,TMASS): C ZP=1 !PROJECTILE CHARGE PMASS=PROTON !PROJECTILE MASS C C FOR IONIC TARGETS, ATTEMPT TO ACCOUNT FOR MASS OF MISSING ELECTRONS C NZ0=1 !TARGET NUCLEAR CHARGE NELC=1 !NO. OF TARGET ELCTRNS (INC RYD) C ZT=NZ0-NELC+1 !TARGET CHARGE (AS SEEN BY RYD) C TMASS=HYDROGEN !TARGET MASS (INC. RYD) C TMASS=TMASS-(NZ0-NELC)*EMASS !APPROX FOR IONIC TARGETS C C*********************************************************************** C*************************** END USER SET-UP *************************** C*********************************************************************** C C----------------------------------------------------------------------- c write(*,"(/' TARGET MASS=',f14.6,' AMU' & ,5x,'TARGET CHARGE(RYD)=',f4.1)")tmass,zt write(*,"(/' PROJECTILE MASS=',f10.6,' AMU' & ,5x,'PROJECTILE CHARGE=',f5.1/)")pmass,zp if(iout.gt.0)write(iout,"(/'# TARGET MASS=',f14.6,' AMU' & ,5x,'TARGET CHARGE(RYD)=',f4.1)")tmass,zt if(iout.gt.0)write(iout,"(/'# PROJECTILE MASS=', f10.6,' AMU' & ,5x,'PROJECTILE CHARGE=',f5.1/)")pmass,zp if(bout7)write(iout7,"(/'# TARGET MASS=',f14.6,' AMU' & ,5x,'TARGET CHARGE(RYD)=',f4.1)")tmass,zt if(bout7)write(iout7,"(/'# PROJECTILE MASS=',f10.6,' AMU' & ,5x,'PROJECTILE CHARGE=',f5.1/)")pmass,zp C C----------------------------------------------------------------------- C C READ TEMPERATURE IN KELVIN AND DENSITY IN /CM^3 FOR ANY CUT-OFF C (DENE.LT.0 CAN FLAG A RANGE OF DENSITIES, SEE BELOW.) C 1 WRITE(0,*) WRITE(0,*)"Enter T(K) and N(/cm^3): (Negative to EXIT)" C READ(*,*)TEE,DENE C IF(TEE.LE.DZERO.OR.DENE.EQ.DZERO)STOP 'WE ARE DONE' C IF(DENE.LT.DZERO)THEN C WRITE(0,*)"Enter No. densities and Log increment (canbe .lt.0):" C READ(*,*)NDEN,DENINC C IF(NDEN.LE.0)STOP 'NEED NDEN.GT.0 FOR DENE.LT.0 FLAG' IF(DENINC.eq.DZERO)STOP 'NEED DENINC.NE.0 FOR DENE.LT.0 FLAG' C DENE=-DENE if(deninc.gt.0)then denmx=log10(dene)+(nden-1)*deninc denmx=dten**denmx else denmx=dene endif ELSE NDEN=1 DENINC=DZERO denmx=dene ENDIF c t2ne=tee*tee/denmx c if(t2ne.lt.tolt2ne)then write(*,"(/' *** Temperature/Density ratio too small', & ' - set T**2/Ne =',1p,e9.2,' .gt.',e9.2)")t2ne,tolt2ne go to 1 endif c if(t2ne.gt.t2nebig)then write(*,"(/' *** SR.PSMX: Temperature/Density ratio too large' & ,' - set T**2/Ne =',1p,e9.2,' .lt.',e9.2)")t2ne,t2nebig go to 1 endif c write(*,*) write(*,"(/' T=',1p,e9.2,' K ',' N=',e9.2,' /cm**3')") & tee,abs(dene) debye=tee/(tkay*8*pi*abs(dene)*bohra**3) !debye radius in a.u. debye=sqrt(debye) write(*,"(/' R_Debye= ',6x,1p,e8.2,' au',3x,1p,e8.2,' cm'/)") & debye,debye*bohra C C USER SPECIFIES RYDBERG L-RANGE AND N-RANGE TO LOOP-OVER AND DELTA-L C (THE RANGES CAN BE IN EITHER DIRECTION, DESPITE MIN,MAX LABELS) C NINC=1 LINC=1 C 10 WRITE(0,*) WRITE(0,*) & "ENTER lmin, lmax, delta-l; nmin, nmax: (5*0 to EXIT)" C READ(*,*)L1MIN,L1MAX,LD0,NMIN,NMAX !,LINC,NINC c write(*,*) C IF(NMIN+NMAX.EQ.0)GO TO 1 C C DIRECTION OF L-LOOP IS GOVERNED BY L1MIN AND L1MAX, SINCE LINC IS C NOT NORMALLY READ, BUT IN CASE READ, CHECK AND RESET IF NECESSARY. C LINC=ABS(LINC) LINC=MAX(IONE,LINC) C C RESTRICT MULTIPOLES FOR PRODUCTION CODE C LD=ABS(LD0) LD=MAX(LD,IONE) IF(LD.GT.1)THEN !FOR PENGELLY SEATON WRITE(0,*)'*** SET MULTIPOLE (LD) .LE.1 FOR PENGELLY SEATON' GO TO 10 ENDIF C C----------------------------------------------------------------------- c c set some (outout) flags c bsum=bout.and.abs(linc).eq.ld !remove bout? bout7=bout7.and.bsum brslv=.false. C C SET-UP FOR UNRESOLVED PS CASE. C LDX=LD0 LDM=LD0 C IF(LD0.EQ.0)THEN !SYNC RESOLVED QM WITH UNRESOLVED PS C write(*,*) write(*,*) & '*** MODIFIED PENGELLY & SEATON USES UNRESOLVED DIPOLE ***' if(iout.gt.0)write(iout, &"(/'# *** MODIFIED PENGELLY & SEATON USES UNRESOLVED DIPOLE ***')" & ) if(bout7)write(iout7, &"(/'# *** MODIFIED PENGELLY & SEATON USES UNRESOLVED DIPOLE ***'/) &") C IF(L1MIN.GT.L1MAX)THEN C C L=L1MIN !REVERSE TO SUPPRESS BACKWARDS RECURSION C L1MIN=L1MAX C L1MAX=L c c comment-out above and uncomment "crev" below to "test" c backwards recursion - a little less accurate than forwards c also need to uncomment "crev" flag l10=-l10-1 above call sr.psmx c *and* in sr.psmx uncomment l1=-l1-1 c crev l1min=l1min+1 ldm=-ld crev ENDIF C IF(L1MIN.LE.L1MAX)THEN L1MIN=L1MIN-1 LDM=LD ENDIF c c see if we can generate resolved ps from unresolved ps by recurrence c brslv=(ldm.gt.0.and.l1min.le.0 & .or.ldm.lt.0.and.l1min+1.ge.max(nmin,nmax)) if(brslv)then ldx=ldm write(*,*) write(*,*) & '*** RESOLVED MODIFIED PENGELLY & SEATON FROM RECURSION ***' write(*,*) if(iout.gt.0)write(iout,"(/ & '# *** RESOLVED MODIFIED PENGELLY & SEATON FROM RECURSION ***' & /)") endif C ENDIF C C ALLOW L TO RANGE IN EITHER DIRECTION (AND CHECK/ADJUST MIN VALUE) C LI=L1MIN LF=L1MAX IF(LI.LE.LF)THEN LI=MAX(-LDX,LI,IZERO) LF0=LF ELSE LF=MAX(-LDX,LF,IZERO) LI0=LI LINC=-LINC ENDIF C C CHECK MAX N FOR PRODUCTION CODE C NMN=MIN(NMIN,NMAX) IF(NMN.GT.MXNVAL)THEN !FOR PRODUCTION CODE WRITE(0,*)'*** NMIN.GT.MXNVAL TOO LARGE FOR PRODUCTION CODE' GO TO 10 ENDIF NMN=MAX(IZERO,NMN) C NMX=MAX(NMIN,NMAX) IF(NMX.GT.MXNVAL)THEN !FOR PRODUCTION CODE WRITE(0,*)'*** REDUCING NMAX TO MXNVAL FOR PRODUCTION CODE' NMX=MXNVAL ENDIF C NINC=MAX(IONE,NINC) C C ALLOW N TO RANGE IN EITHER DIRECTION C IF(NMIN.LE.NMAX)THEN NI=NMN NF=NMX ELSE NI=NMX NF=NMN NINC=-NINC ENDIF c c try warn and flag largest recommended n for the smallest t2ne, c where psm will agree with qm to within 20%. c this is for pure debye cut-off only. c tbd: accommodate lifetime/delta-e cut-off - more severe. c t2=sqrt(t2ne) ntest=5*ntest0*sqrt(t2) if(-nmx.gt.ntest)then write(*, & "(/'(N.B. Debye only) Largest recommended N for this T/Ne is ' & ,i4/)")ntest endif C C----------------------------------------------------------------------- C C *** START DENSITY LOOP C C----------------------------------------------------------------------- C call cpu_time(timei) C DEN0=LOG10(DENE) C DO NE=1,NDEN C TE=DEN0+(NE-1)*DENINC DENE=DTEN**TE c c output headers c if(bout)then debye=tee/(tkay*8*pi*dene*bohra**3) !debye radius in a.u. debye=sqrt(debye) if(iout.gt.0.or.ne.gt.1)then write(iout,20010)tee,dene write(iout,20020)debye,debye*bohra endif if(iout7.gt.0)then write(iout7,20030)tee,dene write(iout7,20040)debye,debye*bohra endif if(ld.eq.1)then if(iout.eq.0)then write(iout,20050) write(iout,20060) else write(iout,20070) write(iout,20080) endif if(iout7.gt.0)then !summed write(iout7,20090) write(iout7,20100) endif else !quadrupole+, resolved onky - maybe summed endif nold=-1 endif C C----------------------------------------------------------------------- C C *** START N-LOOP C C----------------------------------------------------------------------- C DO N=NI,NF,NINC C N1=N-1 C C----------------------------------------------------------------------- C C *** START L-LOOP C C----------------------------------------------------------------------- C IF(LINC.GT.0)THEN LF=MIN(LF0,N1-LDX,N1) ELSE LI=MIN(LI0,N1-LDX,N1) ENDIF C DO L1=LI,LF,LINC C TLIFE=D1P10 !USER CAN SUPPLY RYD LIFETIME IN SEC C DELTAE=D1M20 !USER CAN SUPPLY RYD DELTA-E IN RYDS C IERR=0 !INIT AS .LT. 0 JUST RE-SETS C L10=L1 crev if(ld0.eq.0.and.ldm.lt.0)l10=-l10-1!flag backwards recursion C C CALL PSMX(PMASS,TMASS,ZP,ZT,N,L10,LD0,TEE,DENE,TLIFE,DELTAE & ,QPSMX,QPS64X,IERR) C C IF(IERR.GT.0)THEN WRITE(0,*)'IERR=',IERR STOP '*** FAILURE IN SR.PSMX' ENDIF C C----------------------------------------------------------------------- c c tbd: write info on R_delta-e/life(v_rms) if exists e.g. .lt. 2*R_Debye c c write-out QPS64X,QPSMX rate coefficients (cm^3/s) in illustrative ways c *** the user will want to manipulate and store them according to their c own needs rather than slavishly following what is in sr.output.*** c if(bout)then c call output(iout,iout7,bsum,brslv,nold,pmass,tmass,n,l1 & ,ld0,ld,ldm,linc,tee,debye,tlife,deltae,qpsmx,qps64x) c nold=n endif C C----------------------------------------------------------------------- C ENDDO !END L-LOOP C C----------------------------------------------------------------------- C ENDDO !END N-LOOP C C----------------------------------------------------------------------- C ENDDO !END DENSITY LOOP C C----------------------------------------------------------------------- C C OPTIONAL RETURN TO INITIAL STATE: C IERR=-1 C CALL PSMX(PMASS,TMASS,ZP,ZT,N,L1,LD0,TEE,DENE,TLIFE,DELTAE & ,QPSMX,QPS64X,IERR) C C----------------------------------------------------------------------- c if(iout.gt.0)call flush(iout) !case user exits ungracefully if(iout7.gt.0)call flush(iout7) !case user exits ungracefully c call cpu_time(timef) times=timef-timei write(*,"(/1pe9.2,' secs')")times C GO TO 10 !LETS DO IT ALL AGAIN C C----------------------------------------------------------------------- C 20010 format(/'# T=',1p,e9.2,' K ',' N=',e9.2,' /cm**3') 20020 format(/'# R_Debye= ',1p,e9.2,' au',3x,1p,e8.2,' cm'/) 20030 format(/'# T=',1p,e9.2,' K ',' N=',e9.2,' /cm**3') 20040 format(/'# R_Debye= ',1p,e9.2,' au',3x,1p,e8.2,' cm'/) 20050 format(/20x,'RATE COEFFS (CM**3/S)') 20060 format(3x,'N',2x,'L1->L2',3x,'L',4x & ,'PSM',6x,'PS64',3X,'(PSM-PS64)/PSM') 20070 format(/'#',19x,'RATE COEFFS (CM**3/S)') 20080 format('#',2x,'N',2x,'L1->L2',3x,'L',4x & ,'PSM',6x,'PS64',3X,'(PSM-PS64)/PSM') 20090 format(/'#',15x,'RATE COEFFS (CM**3/S)') 20100 format('#',2x,'N',3x,'L',8x & ,'PSM',6x,'PS64',3X,'(PSM-PS64)/PSM') C C----------------------------------------------------------------------- C END PROGRAM TPSMX C C*********************************************************************** C subroutine output(iout,iout7,bsum,brslv,nold,pmass,tmass,n,l1 & ,ld0,ld,ldm,linc,tee,debye,tlife,deltae,qpsmx,qps64x) c c----------------------------------------------------------------------- c c sr.output illustrates example output based-on analytic qpsmx & qps64x c rate coefficients. c the remaining input is described/defined by the calling routine. c c e.g. c c 1/ the orginal qpsmx,qps64x is written, it may be final-state resolved c or unresolved by final-state. c c 2/ if resolved then l->l' can be summed-over final l'=l+/-1. c c 3/ if unresolved then resolved data can be determined by recursion c from l=0 or l=n-1, if the input l-data *starts* there for a given n. c c note, starting at l=0 and recurring-up is more accurate c overall than starting at l=n-1 and recurring down. c c In addition, fractional differences are constructed between analystic c qpsmx and qps64x. c c----------------------------------------------------------------------- c implicit real*8(a-h,o-z) implicit integer*4 (i-n) c logical bsum,brslv c parameter (dzero=0.0d0) parameter (d1p10=1.d10) parameter (d1m20=1.d-20) c parameter (tkay=1.5789d5) !kelvin to rydbergs parameter (bohra=5.29177d-9) !cm to a.u. parameter (bohrt=2.41889d-17) !secs parameter (emass=5.485798d-4) !electron mass a.m.u c parameter (tolp =2.d-1) !flag when (psmx-ps64x)/psmx.gt. c parameter (xl=0.72) !p&s (30) parameter (xd=1.12) !p&s (29) parameter (idtest=2) c c----------------------------------------------------------------------- c rmass=pmass*tmass/(pmass+tmass) !reduced mass rmass=rmass/emass vrms=sqrt(tee/(rmass*tkay)) !v_rms c dlife=min(tlife,d1p10) dlife=xl*dlife*vrms/bohrt !tlife in sec c ddele=max(deltae,d1m20) ddele=2*xd*vrms/ddele !deltae in ryd c dtest=min(dlife,ddele) c if(dtest.lt.idtest*debye)then !write for info, at least if(dtest.lt.dlife)then write(iout,"(/1x,'R_delta-E(RMS)=' & ,1p,e9.2,' au',3x,1p,e8.2,' cm'/)")ddele,ddele*bohra else write(iout,"(/1x,'R_life(RMS) =' & ,1p,e9.2,' au',3x,1p,e8.2,' cm'/)")dlife,dlife*bohra endif endif c c initialze for summing-over, if resolved c qqm=1 !flag if(bsum.and.n.ne.nold)then qqm0=0 if(brslv.or.ld0.ne.0)then qpsmx0=0 qps64x0=0 endif endif c c basic resolved c if(ld0.eq.0)then c ps is unresolved but qm is still to be summed, so skip until later elseif(ld.eq.1)then !all resolved diffp=dzero if(qpsmx.ne.dzero)diffp=(qpsmx-qps64x)/qpsmx write(iout,"(3i4,5x,1p,2e9.2,1(2x,e9.2))")n,l1,l1+ld0 & ,qpsmx,qps64x,diffp c flag numerics if(abs(diffp).gt.tolp)then write(iout, & "('*** PS64 INACCURACY: (PSM-PS64)/PSM =' & ,1P,E9.2)")diffp endif else endif c c sum resolved (if any) and write-out totals c if(bsum)then !sum resolved c l11=2*l1+1 !symetrize qpsmx=qpsmx*l11 qps64x=qps64x*l11 c c ps may be explicitly unresolved. if so, then c if we start at l1=0 or n-1 we can use recurrence to obtain resolved ps c which is consistent with the unresolved. c lp=l1 if(ldm.lt.0)lp=lp-ld !min(l1,l2) lw=lp if(linc.lt.0)lw=lw+ld lw1=2*lw+1 c if(qqm0.ne.dzero)then if(ld0.ne.0)then !resolved sum qpsmx_sum=(qpsmx+qpsmx0)/lw1 qps64x_sum=(qps64x+qps64x0)/lw1 else !unresolved if(l11.ne.lw1)stop 'l11.ne.lw1' qpsmx_sum=qpsmx/lw1 qps64x_sum=qps64x/lw1 if(brslv)then !can apply recurrence qpsmx=qpsmx-qpsmx0 !to get consistent resolved ps qps64x=qps64x-qps64x0 endif endif else go to 250 endif c 150 continue c iout0=(iout+abs(ld0))*iout7 if(ld.eq.1)then !dipole if(brslv)then diffp=dzero if(qpsmx.ne.dzero)diffp=(qpsmx-qps64x)/qpsmx write(iout,"(3i4,5x,1p,2e9.2,1(2x,e9.2))") & n,l1,l1+ldm,qpsmx/lw1,qps64x/lw1,diffp endif diffp=dzero if(qpsmx_sum.ne.dzero)diffp=(qpsmx_sum-qps64x_sum)/qpsmx_sum if(iout0.eq.0)then write(iout0,"(i4,8x,i4,1x,1p,2e9.2,1(2x,e9.2))")n,lw & ,qpsmx_sum,qps64x_sum,diffp endif if(iout7.gt.0)then write(iout7,"(i4,i4,4x,1x,1p,2e9.2,1(2x,e9.2))")n,lw & ,qpsmx_sum,qps64x_sum,diffp endif c flag numerics if(ld0.eq.0)then if(abs(diffp).gt.tolp.and.ld0.eq.0)then write(iout, & "('*** PS64 INACCURACY: (PSM-PS64)/PSM =' & ,1P,E9.2)")diffp endif endif c else !quadrupole+ endif c 250 continue c if(lp+ld.eq.n-1)then !pick-up last/first single lp=lp+ld lw=lp lw0=lw1 lw1=2*lw+1 if(ld0.ne.0.or.ldm.lt.0.or.brslv)then qpsmx_sum=qpsmx/lw1 qps64x_sum=qps64x/lw1 if(brslv)lw1=lw0 go to 150 else !ld0=0 has one more l-loop for ps qqm0=qqm !for flag endif elseif(lp.eq.0)then !pick-up first/last single lp=-1 lw=0 lw1=1 if(ld0.ne.0.or.ldm.gt.0.or.brslv)then qpsmx_sum=qpsmx qps64x_sum=qps64x if(brslv)lw1=3 go to 150 else !ld0=0 has one more l-loop for ps qqm0=qqm !for flag endif else qqm0=qqm !for flag if(brslv.or.ld0.ne.0)then qpsmx0=qpsmx qps64x0=qps64x endif endif c endif c return c end subroutine output C C*********************************************************************** C SUBROUTINE PSMX(PMASS,TMASS,ZP,ZT,N,L10,LD0,TEE0,DENE0 & ,TLIFE0,DELTAE0,QPSMX,QPS64X,IERR) C C----------------------------------------------------------------------- C C SR.PSMX CALCULATES MAXWELLIAN RATE COEFFICIENTS FOR L-CHANGING C COLLISIONS, I.E. THOSE WHERE THE TRANSITION ENERGY DELTA-E IS C APPROX ZERO. E.G. N L1->N L2, WHERE L2=L1+LD0 AND LD=ABS(LD0)=1 C FOR DIPOLE. LD=0 FLAGS UNRESOLVED (DIPOLE) PS, C RESULTS FOR L1->L1+LD ARE OBTAINED FROM L1+LD->L1 AND THE C APPLICATION OF THE RECIPROCITY RELATION, IN THIS PRODUCTION CODE. C C DIPOLE (ONLY) RESULTS ARE PROVIDED FOR THE C ORIGINAL PENGELLY SEATON (PS64) APPROXIMATION C PENGELLY & SEATON MNRAS 127, 165 (1964) C AND FOR THE C MODIFIED PENGELLY SEATON (PSM20) APPROXIMATION DETAILED BY C BADNELL ET AL MNRAS XXX , XXX (2021), WHICH IMPROVES-UPON THE C ORIGINAL MODIFIED PENGELLY SEATON (PSM17) INTRODUCED BY C GUZMAN ET AL MNRAS 464, 312 (2017) C ALL RESULTS GIVEN ARE EVALUATED FROM THE ANALYTIC FORMULAE OF REFS. C C INPUT (COMPULSORY): C C PMASS: IS THE PROJECTILE MASS IN AMU. E.G. =1.0072764 FOR A PROTON. C C TMASS: IS THE TARGET MASS IN AMU, CONTAINING A RYDBERG ELECTRON C .E.G. =1.007825 FOR THE HYDROGEN ATOM. C C ZP: IS THE CHARGE OF THE PROJECTILE E.G. PROTON=1, ALPHA PARTICLE=2 C C ZT: IS THE CHARGE OF THE TAGRET **AS SEEN BY THE RYDBERG ELECTRON** C SO ZT=1 FOR NEUTRAL HYDROGEN, AND ZT=2 FOR HE+. C C N: IS THE PRINCIPAL QUANTUM NUMBER OF THE RYDBERG ELECTRON, C RESTRICTED TO BE .LE. 1000 FOR THE PRODUCTION CODE. C C L10:IS THE INITIAL ORBITAL ANGULAR MOMENTUM OF THE RYDBERG ELECTRON C C LD0:DENOTES THE 2**LD-POLE ELECTRIC MULTIPOLE TRANSITION OPERATOR. C WHERE LD=ABS(LD0)=1 FOR DIPOLE AND LD=0 FLAGS UNRESOLVED PS. C C TEE: IS THE ELECTRON TEMPRATURE IN KELVIN USED TO CALCULATE THE C MAXWELLIAN RATE COEFFICIENTS. C C DENE: IS THE *ELECTRON* NUMBER DENSITY (*NOT* HYDROGEN DENSITY) C IN /CM^3 USED TO DETERMINE THE DEBYE CUT-OFF FOR THE C LARGE IMPACT PARAMETER. C C IERR: .LT.0 FLAGS DE-ALLOCATION OF PREVIOUSLY ALLOCATED ARRAYS AND C RETURNS THE ROUTINE TO ITS ORIGINAL UNINITIALIZED STATE. C .GE.0 IS OTHERWISE IGNORED. C C INPUT (OPTIONAL): C C TLIFE: IS THE LIFETIME, IN SECONDS, OF THE RYDBERG NL-STATE. C C DELTAE: IS THE DE/EXCITATION ENERGY, IN RYD, OF THE L-CHANGING C TRANSITION, .LE. ZERO IS IGNORED. C C TLIFE AND/OR DELTAE PROVIDE AN ALTERNATIVE EFFECTIVE CUT-OFF, C FOLLOWING PENGELLY & SEATON MNRAS 127, 165 (1964). C THE MINIMUM CUT-OFF VALUE IS USED. C C N.B. THE INPUT (EXCLUDING IERR) IS NOT REDEFINED. C C OUTPUT C C QPS64X: IS THE ORIGINAL PENGELLY-SEATON (PS64) RATE COEFFICIENT C FOR THE TRANSITION N L1->L1-LD IN CM^3/S (ANALYTIC). C C QPSMX: IS THE MODIFIED PENGELLY-SEATON (PSM20) RATE COEFFICIENT C FOR THE TRANSITION N L1->L1-LD IN CM^3/S (ANALYTIC). C C IERR: .GT. ZERO FLAGS FAILURE, USUALLY DUE TO USER INPUT. C C----------------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) C LOGICAL BMOD,BPLTOT,BREV,BMIN,BTWO,BTWOH,BCUT,BFIRST & ,bexp,bout,bex C PARAMETER (MXNVAL=10000) !MAX PRINCIPAL QUANTUM NO. ALLOW C C SET THE USUAL CONSTANTS C PARAMETER (IZERO=0) PARAMETER (IONE=1) C PARAMETER (ITWO=2) !CONFLICTS WITH BTWO VARIABLE!! C PARAMETER (DZERO=0.0D0) PARAMETER (DHALF=0.5D0) PARAMETER (D3QRT=0.75D0) PARAMETER (D3HALF=1.5D0) PARAMETER (DONE=1.0D0) PARAMETER (DTWO=2.0D0) PARAMETER (DTHREE=3.0D0) C PARAMETER (D1P10=1.D10) PARAMETER (D1M20=1.D-20) C PARAMETER (EULER=0.57721566490153286D0) C C SET SOME I/O CONVERSION FACTORS C PARAMETER (EMASS=5.485798D-4) !ELECTRON MASS A.M.U PARAMETER (TKAY=1.5789D5) !KELVIN TO RYDBERGS PARAMETER (BOHRA=5.29177D-9) !CM TO A.U. PARAMETER (BOHRT=2.41889D-17) !SECS C PARAMETER (CQ=BOHRA**3/BOHRT) !a_0^3/tau_0 C PARAMETER (T2NEBIG=1.D+25) !MAX ALLOWED VALUE OF T**2/Ne PARAMETER (TOLT2NE=1.D-15) !MIN ALLOWED VALUE OF T**2/Ne C SAVE C DATA BFIRST/.TRUE./,ierrout/0/,iun6/-1/,debyo/dzero/ C C----------------------------------------------------------------------- C IF(IERR.LT.0)THEN !ALLOW USER TO RE-SET BFIRST=.TRUE. IERR=0 RETURN ENDIF C C----------------------------------------------------------------------- C IF(BFIRST)THEN C BFIRST=.FALSE. C C----------------------------------------------------------------------- C C SET SOME CONSTANTS ETC C PI=ACOS(-DONE) RTPI=SQRT(PI) C C*********************************************************************** C************ USER JOE SHOULD NOT CHANGE ANYTHING BELOW HERE ********** C*********************************************************************** C C----------------------------------------------------------------------- C C I/O SET-UP C (SET UNIT NUMBER NEGATIVE TO SUPPRESS WRITES THERE) C iout=-4 !=0 rate coeffs to screen bout=iout.ge.0 !more details/debug than tpsmx.out in main if(iout.eq.6)inquire(file='tpsmx.out',number=iun6) !=6 to outfile if(iout.gt.0.and.iun6.ne.6)open(iout,file='psmx.out') C C----------------------------------------------------------------------- C C ALLOW FOR VARIATION FROM USUAL PENGELLY & SEATON (1964): C C BMOD=.T. THEN PHALF=1/2 PMOD=2/3 GIVES MODIFIED; FIXED B-RATIO=1/2. C CAN ALSO SWITCH SMALL R BEHAVIOR FROM ORIGINAL R/R_C TO R/R_1 C C BMOD=.F. THEN PHALF=1/2 PMOD=1 GIVES USUAL P&S; L-DEPENDENT B-RATIO. c can be used to test effect of btwo and bmin on standard P&S c since "psm" is then original ps64, save exact exp. int. used. C BMOD=.TRUE. ! ***DO NOT CHANGE*** C IF(BMOD)THEN !***NOTE CHANGE IN DEF OF PHALF! C PHALF=DHALF !R1 DEFINED BY P=PONE=PHALF*B.R. IPOW=1 !IF NON-INTEGER DEFINE AS REAL*8 PMOD=DTWO/(2+IPOW) !THEN P=PONE*(R/R1)**IPOW :R1 NOT CODED ***' C C----------------------------------------------------------------------- C C BMIN=.T./.F. FLAGS +/- CONTRIB FROM E=0-E_MIN FOR (MODIFIED) P&S C bmin=.False. !STANDARD P&S NEGLECTS 0-E_EMIN c BMIN=bmin.or.BMOD !PSM INCLUDES IT C C----------------------------------------------------------------------- c c bexp=.f. PS64 uses original P&S approx of exponential integral c bexp=.t. PS64 uses exact exponential integral E_min-infty. c not referenced by PSM (even if BMOD=.False.) c bexp=.False. C if(bout)then if(bexp)then write(iout,*)'*** SR.PSMX: ' & ,'Exact Exponential integral in use for PS64' else write(iout,*)'*** SR.PSMX: ' & ,'Approx. Exponential integral in use for PS64 (only)' endif endif C C----------------------------------------------------------------------- C C BPLTOT=.F. GENERATES FINAL STATE RESOLVED DATA WITH APPROPRIATE C BRANCHING RATIOS AND MATCHING POINTS. C C BPLTOT=.T. GENERATES TOTAL UNRESOLVED BY FINAL STATE DATA SO UNIT C BRANCHING AND MATCHING POINT BASED ON TOTAL PROBABILTY. C NOTE, THIS DOES *NOT* EQUAL L-SUMMED MODIFIED P&S, WHICH C USES DIFFERENT MATCHING POINTS FOR L-1,L+1; NOR ORGINAL C P&S DETERMINED BY RECIPROCITY. C BPLTOT=LD0.EQ.0 C if(bout.and.bpltot)write(iout,*) x'*** NOTE: PENGELLY & SEATON USES UNRESOLVED DIPOLE TERM ***' C C----------------------------------------------------------------------- C C BREV=.T. CALCULATE L->L+1 FROM L+1->L AND RECIPROCITY C =.F. CALCULATE L->L+1 EXPLICITLY *AND* WARN - NOT FOR PRODUCTION C BREV=.TRUE. ! ***DO NOT CHANGE*** C IF(.NOT.BREV.AND.LD0.GT.0)THEN if(bout)write(iout,*) & '*** SR.PSMX: PSM WARNING, FOR L -> L+1 YOU SHOULD CALCULATE' & ,' L+1 -> L AND APPLY DETAILED BALANCE ***' ENDIF C C----------------------------------------------------------------------- C C ENERGY DEPENDENT CUT-OFF WILL EXCEED DEBYE FOR LARGE ENOUGH ENERGY. C C BTWO=.TRUE.COMBINES USE OF BOTH, SWITCHING OVER AT THE MATCHING ENERGY C SEE GUZMAN ET AL (2017) OR BADNELL ET AL (2021) FOR DETAILS C THIS IS APPLIED TO ALL NUMERICAL AND ANALYTIC PSM METHODS. C (IF R(E=V_RMS).GT.IDTEST*R(DEBYE) IT IS RE-SET .FALSE.) C C BTWO=.FALSE. USES DEBYE OR ENERGY DEPENDENT AT ALL ENERGIES, C DEPENDING ON WHICH IS THE SMALLER AT V_RMS. C btwo=.False. !STANDARD P&S USES ONE OR OTHER C BTWO=BTWO.OR.BMOD !PSM USES COMBINED C if(.not.btwo.and.bout)write(iout,*)'*** ATTENTION: CUT-OFF IS ' & ,'DEBYE *OR* LIFETIME/DELTA-E **ONLY**, *NOT* COMBINED ***' C IF(BTWO)THEN IDTEST=2 !USE PURE DEBYE WHEN R(VRMS).GT.IDTEST*R(DEBYE) ELSE IDTEST=1 !AS NO COMBINED CASE POSSIBLE ENDIF C BTWOH=BTWO !HOLD TO ALLOW BTWO RE-SET .F. C C SET XL & XD FOR POSSIBLE LIFETIME & DELTA-E USAGE C x0=d3qrt !use 3*V_RMS/4 X0=DONE !P&S C X=0.72 !P&S (30) XL=X*X0 C X=1.12 !P&S (29) XD=X*X0 C C----------------------------------------------------------------------- C C BCUT=.T. APPLIES LARGE I.P. CUT-OFF TO ALL DELTA-L TRANSITIONS. C C BCUT=.F. APPLIES TO DIPOLE ONLY; (EFFECTIVE) COLLISION STRENGTHS C FOR DELTA-L.GT.1 ARE THEN THEIR CONSTANT LIMITING VALUE. C "ONLY" REQUIRED FOR RCUT.LT.~0.1 CM, I.E. C FOR H, VERY LOW-T/HIGH DENSITY; C FOR HE, AND/OR HIGHLY NON-DEGENERATE. C C RCUT0 .GT. 0 IS ANY USER-SUPLLIED TEST CUT-OFF (IN A.U.) C .LT. 0 (ANY) ALLOWS REDEFINITION BY DEBYE/LIFETIME/DELTA-E C BCUT=.TRUE. C RCUT0=-DONE !ALLOWS CUT-OFF TO BE REDEFINED C BCUT=RCUT0.NE.DZERO.AND.(BCUT.OR.ABS(LD0).LE.1) !SET I.P. CUT-OFF C IF(.NOT.BCUT)THEN IF(ABS(LD0).LE.1)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: MUST FLAG AN I.P. CUT-OFF FOR DIPOLE TRANSITIONS' IERR=12 RETURN ENDIF if(bout)WRITE(iout,*) & '*** SR.PSMX: ANY I.P. CUT-OFF IGNORED FOR DELTA-L=',ABS(LD0) ENDIF C C----------------------------------------------------------------------- C C*********************************************************************** C************ USER JOE SHOULD NOT CHANGE ANYTHING ABOVE HERE ********** C************ AND CERTAINLY SHOULD NOT CHANGE ANYTHING BELOW ********** C*********************************************************************** C C ENDIF !END INITIALIZATION C C----------------------------------------------------------------------- C C INITIALIZE OUTPUT INCASE ABORT C QPS64X=DZERO QPSMX=DZERO C C----------------------------------------------------------------------- C C CHECK USER INPUT C IERR=0 C IF(PMASS.LE.DZERO)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: PROJECTILE MASS MUST BE POSITIVE' IERR=101 RETURN ENDIF C IF(NINT(PMASS).LT.IONE)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: PROJECTILE MASS MUST BE IN A.M.U.' IERR=102 RETURN ENDIF C IF(TMASS.LE.DZERO)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: TARGET MASS MUST BE POSITIVE' IERR=103 RETURN ENDIF C IF(NINT(TMASS).LT.IONE)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: TARGET MASS MUST BE IN A.M.U.' IERR=104 RETURN ENDIF C IF(ZP.LE.DZERO)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: PROJECTILE CHARGE MUST BE POSITIVE' IERR=201 RETURN ENDIF C IF(ZT.LE.DZERO)THEN !ZT IS CHARGE AS SEEN BY RYDBERG if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: TARGET CHARGE MUST BE NON-NEGATIVE' IERR=202 RETURN ENDIF C IF(N.LT.0)THEN if(ierrout.ge.0)WRITE(ierrout,*) WRITE(0,*)'*** SR.PSMX: N MUST BE NON-NEGATIVE' IERR=301 RETURN ENDIF C IF(N.GT.MXNVAL)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: N (>MXNVAL) TOO LARGE FOR PRODUCTION CODE' IERR=302 RETURN ENDIF C L1=L10 crev if(l10.lt.0)l1=-l1-1 !re-set l1 (backwards recursion C IF(L1.LT.0)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: L MUST BE NON-NEGATIVE' IERR=401 RETURN ENDIF C IF(L1.GE.N)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: L MUST BE LESS THAN N' IERR=402 RETURN ENDIF C IF(ABS(LD0).GT.3)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: MULTIPOLE (>3) TOO LARGE FOR PRODUCTION CODE' IERR=403 RETURN ENDIF C IF(TEE0.LE.DZERO)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: TEMPERATURE (KELVIN) MUST BE POSITIVE' IERR=501 RETURN ENDIF C IF(DENE0.LE.DZERO)THEN if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: DENSITY (/CM3) MUST BE POSITIVE' IERR=502 RETURN ENDIF C C TESTS ON USER TLIFE0 & DELTAE0? C C... C C----------------------------------------------------------------------- C C INPUT O.K. SO NOW CONVERT UNITS ETC C ZP2=ZP*ZP !PROJECTILE CHARGE SQUARED ZT2=ZT*ZT !TARGET CHARGE SQUARED C RMASS=PMASS*TMASS/(PMASS+TMASS) !REDUCED MASS RMASS=RMASS/EMASS C TEE=TEE0 DENE=DENE0 c t2ne=tee*tee/dene c if(t2ne.lt.tolt2ne)then write(0,"(/' *** SR.PSMX: Temperature/Density ratio too small' & ,' - set T**2/Ne =',1p,e9.2,' .gt.',e9.2)")t2ne,tolt2ne ierr=1 return endif c if(t2ne.gt.t2nebig)then write(0,"(/' *** SR.PSMX: Temperature/Density ratio too large' & ,' - set T**2/Ne =',1p,e9.2,' .lt.',e9.2)")t2ne,t2nebig ierr=1 return endif C TEE=TEE/TKAY !KELVIN TO RYD C DENE=DENE*BOHRA**3 !CM TO A.U. C VRMS2=TEE/RMASS VRMS=SQRT(VRMS2) !V_RMS C CQQ=CQ*2*RTPI/VRMS !FOR PROBABILITY TO RATE COEFF. C TCM=2*RMASS*RMASS/CQQ !AS SEATON INCLUDES RMASS^(-3/2) C TCM=TCM/(AN*AN) !FOR UPSILON/N^4 C C----------------------------------------------------------------------- C C DEBYE CUT-OFF FOR USER T/NE (PRINCIPALLY FOR DIPOLE TRANSITIONS) C C----------------------------------------------------------------------- C DEBYE2=TEE/(8*PI*DENE) DEBYE=SQRT(DEBYE2) !DEBYE RADIUS IN A.U. c if(debye.ne.debyo)then if(bout)write(iout,"(/' R_Debye= ',6x,1p,e9.2,' au',3x,1p & ,e8.2,' cm')")debye,debye*bohra endif c if(bout.and.debye.ne.debyo)then if(ld0.eq.0)then write(iout,"(/16x,'RATE COEFFS (CM**3/S)')") write(iout,"(3x,'N',3x,'L',8x & ,'PSM',6x,'PS64',3X,'(PSM-PS64)/PSM')") elseif(abs(ld0).eq.1)then write(iout,"(/20x,'RATE COEFFS (CM**3/S)')") write(iout,"(3x,'N',2x,'L1->L2',3x,'L',4x & ,'PSM',6x,'PS64',3X,'(PSM-PS64)/PSM')") else endif endif c debyo=debye C C----------------------------------------------------------------------- C C SELECT CUT-OFF: C BTWO=BTWOH !RE-INSTATE C ITWO=1 !INITIALIZE FOR ENERGY IND DEBYE C C RCUT0.GT.0 IS ANY USER GLOBAL (IND. NL) ENERGY IND. I.P. CUT-OFF C FLAG RCUT0 NEGATIVE TO USE MIN OF DEBYE/DELTAE/LIFETIME C IF(RCUT0.LT.DZERO)THEN C RCUT0=-DEBYE !INITIALIZE FOR DEBYE C C----------------------------------------------------------------------- C C LIFETIME CUT-OFF BASED ON USER SUPPLIED TLIFE0 FOR THIS NL C C----------------------------------------------------------------------- C TLIFE=MIN(TLIFE0,D1P10) TLIFE=TLIFE/BOHRT !SEC TO A.U. DLIFE=XL*TLIFE*VRMS C C----------------------------------------------------------------------- C C ENERGY CUT-OFF BASED ON USER SUPPLIED DELTAE0 IN RYDS FOR NL->L+LD0 C C----------------------------------------------------------------------- C DELTAE=MAX(DELTAE0,D1M20) beta=5*deltae*n*n beta=5*beta*beta if(tee.lt.beta)then if(ierrout.ge.0)WRITE(ierrout,*) & '*** SR.PSMX: Delta-e too large for kT, increase kT to' & ,nint(beta*tkay),' K' IERR=13 GO TO 999 endif C DDELE=2*XD*VRMS/DELTAE !DELTAE IN RYD C C----------------------------------------------------------------------- C C USE MIN OF ENERGY IND. DEBYE AND ENERGY DEP. LIFETIME/DELTA-E. C C----------------------------------------------------------------------- C DTEST=MIN(DLIFE,DDELE) C IF(DTEST*DTEST.LT.-IDTEST*RCUT0*DTEST)THEN if(bout)then if(dtest.lt.dlife)then write(iout,"(/' R_delta-E(RMS)=',1p,e9.2,' au',3x,1p,e8.2 & ,' cm')")ddele,ddele*bohra else write(iout,"(/' R_life(RMS) =',1p,e9.2,' au',3x,1p,e8.2 & ,' cm')")dlife,dlife*bohra endif endif RCUT0=-DTEST !HOLD FLAG ITWO=2 !USE E-DEP RCUT (BTWO=.T.OR.F.) ELSE BTWO=.FALSE. !AS PURE DEBYE ENDIF C ELSE C BTWO=.FALSE. C ENDIF C C----------------------------------------------------------------------- c c define an effective deltae based-on rcut (for poss. future use.) c if(rcut0.ne.dzero)deltae=2*xd*vrms/abs(rcut0) C C----------------------------------------------------------------------- C C WE HAVE FINAL RCUT CUT-OFF C RCUT=ABS(RCUT0) RCBAR=RCUT*RCUT !SQUARED c rmean=(3*n*n-l1*(l1+1))/(2*zt) !mean radius if(rcut*rcut.lt.rcut*rmean)then if(rcut.lt.debye*0.999D0)then !warn if(ierrout.ge.0) & WRITE(ierrout,"('*** SR.PSMX: R_cut(RMS)=',1p,e9.2,' au',3x, & '.LT. Mean radius=',1p,e9.2,' au'/)")rcut,rmean else !abort if(ierrout.ge.0) & WRITE(ierrout,"('*** SR.PSMX: R_Debye=',1p,e9.2,' au',3x, & '.LT. Mean radius=',1p,e9.2,' au'/)")rcut,rmean IERR=14 GO TO 999 endif endif C C DEFINE MATCHING ENERGY FOR COMBINED LIFE/DELTA & DEBYE: C RCUT(DEBYE)=RCUT(LIFETME) C IF(BTWO)THEN EEC=DEBYE2/RCBAR !IS /KT c write(0,*)'eec=',eec ELSE !JUST INITIALIZE IF(ITWO.EQ.1)THEN !DEBYE EEC=DZERO ELSE !LIFE/DELTA-E EEC=D1P10 ENDIF ENDIF C C----------------------------------------------------------------------- C C NL SET-UP C C----------------------------------------------------------------------- C AN=N*N C LD=ABS(LD0) LD=MAX(LD,IONE) C IF(BREV.AND.LD0.GT.0)THEN !IF L->L+1 THEN LDM=-LD !CALCULATE L+1->L L1=L1+LD !AND APPLY RECIPROCITY ELSE LDM=SIGN(LD,LD0) !SO L->L+1 CASE LD0=0 if(l10.lt.0)ldm=-ld !test l->l-1 ENDIF C L11=2*L1+1 W1=L11 !INITIAL L STAT WEIGHT L2=L1+LDM L22=2*L2+1 W2=L22 !FINAL L STAT WEIGHT W12=W1/W2 LGT=MAX(L1,L2) C C----------------------------------------------------------------------- C C INITIALIZE FOR PS: C C----------------------------------------------------------------------- C IF(LD.EQ.1)THEN !SET-UP FOR PSM & PS64 SLINIJ=9*AN*LGT*(AN-LGT*LGT)/(2*ZT2) !RESOLVED SLINI=9*AN*(AN-L1*(L1+1)-1)/(2*ZT2) !UNRESOLVED IF(BPLTOT)THEN !SET HISTORIC SLIN FOR P&S, COULD JUST SET PPS0 W=DONE SLIN=SLINI BR=DONE ELSE W=L11 SLIN=SLINIJ IF(BMOD)THEN BR=DHALF !GUZMAN ET AL PSM c br=slinij/slini !summers, storey etc c br=br/w ELSE BR=SLINIJ/SLINI BR=BR/W ENDIF ENDIF PPS0=2*ZP2*SLIN/(3*W) !includes AS spin factor 2 RONE=PPS0/BR RONE=RONE/PHALF !=(R1*V)**2 C PONE=PPS0/RONE !=PHALF*BRANCHING RATIO R1BAR=RONE/VRMS2 !SQUARED EMIN=R1BAR/RCBAR !R1=RCUT (IS /KT) if(itwo.eq.2)emin=sqrt(emin) !energy-dependent emin EMIND=R1BAR/DEBYE2 !DEBYE EMIN/KT ENDIF C C----------------------------------------------------------------------- C C NOW EVALUATE ANALYTIC PSM & PS64 RATE COEFFICIENTS C C----------------------------------------------------------------------- C IF(LD.EQ.1)THEN !DIPOLE C C CONTRIB E_MIN->\infty C E1=ENINT(IONE,.FALSE.,EMIN) !USE EXACT FOR E_MIN NOT SMALL C TX=EXP(-EMIN) TT=PMOD*TX+ITWO*E1 !DEBYE OR LIFE/DELTA-E, ALL PS C IF(BTWO)THEN !COMBINED IF(EEC.GT.EMIN)THEN !ALL LIFE/DELTA-E TT=TT-ENINT(IONE,.FALSE.,EEC) ELSE !TRUNCATE BY DEBYE (UNLIKELY) TX=EXP(-EMIND) TT=PMOD*TX+ENINT(IONE,.FALSE.,EMIND) ENDIF ENDIF C C CONTRIB 0->E_MIN C IF(BMIN)THEN C IF(IMIN.EQ.1)THEN C C PSM20 IF(BTWO)THEN !COMBINED C IF(EEC.LT.EMIN)THEN !TRUNCATE BY DEBYE (UNLIKELY) c c write(0,*)"Note, matching energy E_c .lt. E_min" c & ,eec,emin TC=EXP(-EEC) T=6*(DONE-TC*(DONE+(DONE+EEC/4)*EEC))/EMIN**3 RTED=SQRT(EMIND) RTEC=SQRT(EEC) T0=D3QRT*RTPI*(ERF(RTED)-ERF(RTEC)) T0=T0/RTED-TX*(D3HALF+EMIND) T0=T0/EMIND+T C ELSE !ALL LIFE/DELTA-E C T0=(DONE-TX*(DONE+EMIN))/EMIN T0=DTWO*T0/EMIN-TX T0=DTHREE*T0/EMIN-TX c c write(0,*)"Note, matching energy E_c .gt. E_min" c & ,eec,emin ENDIF C ELSE !ONE OR OTHER C IF(ITWO.EQ.1)THEN !DEBYE C RTEM=SQRT(EMIN) T0=D3QRT*RTPI*ERF(RTEM)/RTEM-TX*(D3HALF+EMIN) T0=T0/EMIN C ELSE !LIFE/DELTA-E C T0=(DONE-TX*(DONE+EMIN))/EMIN T0=DTWO*T0/EMIN-TX T0=DTHREE*T0/EMIN-TX C ENDIF C ENDIF C ELSE C C PSM17 C IF(BTWO)THEN !COMBINED C IF(EEC.LT.EMIN)THEN !TRUNCATE BY DEBYE (UNLIKELY) c c write(0,*)"Note, matching energy E_c .lt. E_min" c & ,eec,emin TC=EXP(-EEC) T=(DTWO-TC*(DTWO+EEC))/EMIN T0=T/EMIN-TX*(DONE+DONE/EMIND) C ELSE !ALL LIFE/DELTA-E C T0=(DONE-TX*(DONE+EMIN))/EMIN T0=DTWO*T0/EMIN-TX c c write(0,*)"Note, matching energy E_c .gt. E_min" c & ,eec,emin ENDIF C ELSE !ONE OR OTHER C T0=(DONE-TX*(DONE+EMIN))/EMIN !DEBYE IF(ITWO.EQ.2)T0=DTWO*T0/EMIN-TX !LIFE/DELTA-E C ENDIF C ENDIF C T0=PMOD*T0 TT=TT+T0 C ENDIF C C STANDARD PS64 OMITS 0-E_MIN AND DOES NOT COMBINE CUT-OFFS c also, it approximates the exponential integral, although c its validity can be extended somewhat by not doing so, hence c if(bexp)then tt64=tx+itwo*e1 else tt64=done-itwo*(euler+log(emin)) endif C QPSMX=CQQ*PPS0*TT UPSMX=QPSMX*TCM*L11 QPS64X=CQQ*PPS0*TT64 UPS64X=QPS64X*TCM*L11 C IF(LDM*LD0.LT.IZERO)THEN !APPLY RECIPROCITY FOR PRODUCTION QPSMX=QPSMX*W12 QPS64X=QPS64X*W12 ENDIF c if(bout)then diffp=dzero if(qpsmx.ne.dzero)diffp=(qpsmx-qps64x)/qpsmx if(ld0.eq.0)then write(iout,"(2i4,5x,1p,2e9.2,1(2x,e9.2))")n,l1 & ,qpsmx,qps64x,diffp c & ,upsmx,ups64x,diffp else write(iout,"(3i4,5x,1p,2e9.2,1(2x,e9.2))")n,l1,l2 & ,qpsmx,qps64x,diffp c & ,upsmx,ups64x,diffp endif endif C ELSEIF(LD.EQ.2)THEN !QUADRUPOLE C ENDIF c if(bout)call flush(iout) !case user exits ungracefully C C----------------------------------------------------------------------- C 999 CONTINUE C RETURN C END SUBROUTINE PSMX C C*********************************************************************** C REAL*8 FUNCTION ENINT(N,BHALF,X) C C----------------------------------------------------------------------- C C EVALUATES EXPONENTIAL INTEGRAL C N R BADNELL: ADAPTED FROM NUMERICAL RECIPES AND ABRAMOWITZ & STEGUN C C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) C LOGICAL BHALF C PARAMETER (DZERO=0.0D0) PARAMETER (DONE=1.0D0) PARAMETER (DTWO=2.0D0) PARAMETER (BIG=1.D30) PARAMETER (TOLE=1.D-5) PARAMETER (MAXIT=100) PARAMETER (EULER=0.5772156649D0) C C GENERATE INTEGER (BHALF.EQ.FALSE.) EXPONENTIAL INTEGRAL ORDER N C GENERATE HALF-INTEGER (BHALF.EQ.TRUE.) EXPONENTIAL INTEGRAL ORDER N+.5 C IF(N.LT.0.OR.X.LT.DZERO.OR.X.EQ.DZERO.AND.(N.EQ.0.OR.N.EQ.1))THEN WRITE(0,*)'ILLEGAL ARGUMENT IN FN.ENINT: N,X=',N,X ENINT=DZERO RETURN ENDIF C NM=N-1 C H=DZERO IF(BHALF)H=DONE/DTWO C IF(N.EQ.0.AND..NOT.BHALF)THEN EXPINT=EXP(-X)/X ELSEIF(X.EQ.DZERO)THEN EXPINT=DONE/(NM+H) ELSEIF(X.GT.DONE)THEN !CONTINUED FRACTION BX=X+N+H CX=BIG DX=DONE/BX HX=DX DO I=1,MAXIT AX=-I*(NM+H+I) BX=BX+DTWO DX=DONE/(AX*DX+BX) CX=BX+AX/CX DEL=CX*DX HX=HX*DEL IF(ABS(DEL-DONE).LT.TOLE)THEN EXPINT=HX*EXP(-X) GO TO 100 ENDIF ENDDO c write(0,*)'FN.ENINT continued fraction not converged' EXPINT=-DONE ELSE !SERIES FACT=DONE IF(BHALF)THEN PI=ACOS(-DONE) FACH=-SQRT(PI*X)/(N-H) ENDIF IF(NM.LE.0)THEN IF(BHALF)THEN IF(NM.EQ.0)THEN EXPINT=FACH+DTWO ELSE EXPINT=SQRT(PI/X)-DTWO ENDIF ELSE EXPINT=-LOG(X)-EULER ENDIF ELSE EXPINT=DONE/(NM+H) ENDIF DO I=1,MAXIT FACT=-FACT*X/I IF(BHALF)FACH=-FACH*X/(I-H) IF(I.NE.NM)THEN DEL=-FACT/(I-NM-H) ELSE IF(BHALF)THEN DEL=FACH+DTWO*FACT ELSE !DIGAMMA PSI=-EULER DO II=1,NM PSI=PSI+DONE/II ENDDO DEL=FACT*(-LOG(X)+PSI) ENDIF ENDIF EXPINT=EXPINT+DEL IF(ABS(DEL).LT.ABS(EXPINT)*TOLE)GO TO 100 ENDDO c write(0,*)'ENINT series not converged' EXPINT=-DONE ENDIF C 100 ENINT=EXPINT C RETURN C END FUNCTION ENINT C C***********************************************************************