c func-photo.f c functions and routines to compute photoionation rates c routines: c FUNCTION PH_Rate(k,j) - compute photoionization rates c FUNCTION thresholdE(k,nek,s) - obtain threshold energy c SUBROUTINE rate_integral(alpha) - performs the rate integrate c FUNCTION funk(E) - the integrand of the rate integral c FUNCTION FlogJEsed(E) - obtain 2nd derivitives of SED c FUNCTION Flogsig(E) - obtain 2nd derivitives of cross sections c c......................................................................... c DOUBLE PRECISION FUNCTION PH_Rate(k,j) c c author: Chris Churchill c this routine has two purposes: c (1) compute the partial photoionization rate for species k in c ionization stage j (by partial we mean for each individual shell); c store these in the array R_phs(s), where "s" is the shell c number (see below); the partial rates will be used to compute the c Auger rates c (2) compute the total photoionization rate (summed over all shells) c for single electron ejection (no Auger electrons); return this c rate as PH_Rate c the partial photoionization rate coefficient is determined by c integrating the ionizing spectrum modulated by the partial c photoionization cross section (units cm^2) over energy from the c threshold energy for the shell to the maximum energy of the c ionizing spectrum c the shell numbers are indexed (by s) as follows c s=1 1s (l=0) c s=2 2s (l=0) c s=3 2p (l=1) c s=4 3s (l=0) c s=5 3p (l=1) c s=6 3d (l=2) c s=7 4s (l=0) c the partial cross sections are obtained via a call to a modified c version of the code, phfit2.f, written by D. Verner; it has now c been thoroughly commented and modified; mods include changing it c from a subroutine call to a function call (thus the argument c returned, the cross section, has been removed from the call c parameter list); the module is now called phfit2-mod.f, see c comments therein for details (especially about the bug I found!) c after obtaining the cross section as a function of energy, we then c perform the integration c ratekj(s) = 4*pi * integral{sigkj(s,E)*J(E)*dE/E} E=E0->Emax c to obtain the rate coefficient (units cm^3/s)m where J(E) is the c mean intensity of the ionizing spectrum (per unit energy); the c integration is performed in routine rate_integral; in practice we c perform the integration in parts c we assume k such that the species numbers k go as 1=H, 2=He, etc, c and the ionization stages go as 1=neutral, 2=single, 3=double c ionized etc.; the ion is assumed to be in its ground state; our c formalism is to index each ion by (k,j), whereas Verner's notation c for his phfit2 routine indexes each element as (nz,ne); the c translation is c OUR'S VERNER'S c k nz atomic number c nek ne number of electrons bound to ion c s is shell number c such that Verner's indices in terms of our indices are c nz = k = atomic number c ne = nek = k - j + 1 = number of bound electrons c is = s c c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c implicit none include 'rates.h' integer i,j,k,nek,s,nout,Nsplit double precision logElo,logEhi,dlogE,Elo,E,pow double precision Eth(NShmax) double precision alphakj,aphsum double precision phfit2,thresholdE include 'rates.com' include 'com-photo.com' include 'com-auger.com' c we need the threshold energies for each shell, so we include the c ph1 common block from Verner's block data; we also need the total c number of populated shells, so include the ntot common block double precision ph1 COMMON/ph1/ ph1(6,30,30,7) integer ninn,ntot common/ninn/ ninn(30) ! inner shell COMMON/ntot/ ntot(30) ! outer shell c end declarations c ........... c SET UP SHOP c ........... c determine the number of bound electrons, nek, and the index of the c most outer shell populated by one or more electron(s) for this c number of electrons, nout; nout is determined from Verner's ntot c array; it is assumed the ion is in the ground state nek = k - j + 1 nout = ntot(nek) c obtain the threshold energies for the integration for each shell c of k,j DO 11 s=1,nout Eth(s) = thresholdE(k,nek,s) 11 CONTINUE c ..................... c BEGIN THE INTEGRATION c ..................... c integrate over cross section and ionizing spectral distribution to c obtain the photoionization rate coefficient one shell at a atime DO 01 s=1,nout aphsum = 0.0d0 c we must trap the special cases for the KI, CaI and CaII electron c configurations (4s1, 4s2) since the shell s=6 (3d) is empty, so c null and advance to the next shell IF (s.eq.6) then IF ((k.eq.19.AND.nek.eq.19).OR. ! KI (4s1) & (k.eq.20.AND.nek.eq.19).OR. ! CaII (4s1) & (k.eq.20.AND.nek.eq.20)) then ! CaI (4s2) R_phs(s) = 0.0d0 GOTO 01 ! advance to next shell ELSE GOTO 02 ! carry on END IF END IF c set up the energy grid for the partial cross section calculation; c Emin is the current shell energy and is the lower energy limit c over which the integration for the rate is performed; Emax is c global and has been returned from routine getJEofz (the maximum c energy of the ionizing spectrum) 02 logEmin = log10(Eth(s)) c if the shell energy is greater than the maximum energy of the SED c array, then then the contribution from this shell is null, skip to c next shell (the shell energies decline as shell number incerases) IF (logEmin.ge.logEsed(NEsed)) GOTO 01 c if the shell energy is less than SED min energy then we have a c situation where ionization from this shell can still occur ans so c it will contribute to the photoionization rate, but we are c missing the lowest energy photons from the SED array, we fix this c by resetting the minimum energy of the integration to the lowest c energy in the SED array; note however, that this represent some c missed contribution to the photoionization rate in principle, but c not from a computational standpoint (i.e., in the real world the c SED energy does not cut off to absolute zero.. however, as can be c seen with stellar only SEDS, the energy drops so rapidly that c this treatment is still an excellent approximation) IF (logEmin.lt.logEsed(1)) logEmin = logEsed(1) c good to go, set up the energy grid dlogE = 0.01 logElo = logEmin logEhi = logEmax + Ebuff NExsec = (logEhi-logElo)/dlogE + 1 c be very careful to get the exact edge, for it may be a sharp c edge and we want to be sure we catch all of it for the c integration; the -18 conversts from Mb to cm^2 Elo = Eth(s) logEsig(1) = log10(Elo) logsigkj(1) = log10(phfit2(k,nek,s,Elo)) - 18.0d0 c populated the partial cross section over energy DO 03 i=2,NExsec pow = float(i-1)*dlogE E = Elo*(10.0d0**pow) logEsig(i) = log10(E) logsigkj(i) = log10(phfit2(k,nek,s,E)) - 18.0d0 03 CONTINUE c obtain the global 2nd derivatives of logsigkj(E) for future c interpolating along the energy direction during the integration c to obtain the rate coefficient CALL Flogsig_derivatives c integrate beyond the current shell edge one energy decade at a c time; this splitting of the integral speeds up convergence because c it is difficult to numerically integrate a function over several c (4-6!) orders of magnitude in energy; variable Nsplit determines c the number of full energy decades to loop over logEmax = logEsed(NEsed) dlogE = logEmax - logEmin Nsplit = int(dlogE) IF (Nsplit.eq.0) Nsplit = 1 DO 04 i=1,Nsplit logEmax = logEmin + 1.0d0 CALL rate_integral(alphakj) aphsum = aphsum + alphakj logEmin = logEmax 04 CONTINUE c now integrate over the final energy decade (this energy interval c will not be a full decade- and the contribution is usually c negligible, though that depends upon the ioninzing SED) logEmin = logEmax logEmax = logEsed(NEsed) CALL rate_integral(alphakj) aphsum = aphsum + alphakj c store the rate for this shell R_phs(s) = aphsum 01 CONTINUE ! increment to next shell c .......................................................... c COMPUTE PHOTOIONIZATION RATE FOR SINGLE ELECTRON EJECTIONS c .......................................................... c the photoionization rate for single electron ejection is the sum c over all the shells with each shell weighted by the probability of c yielding a single electron (no Auger electrons) PH_Rate = 0.0d0 DO 05 s=1,nout PH_Rate = PH_Rate + W_auger(k,j,s,1)*R_phs(s) 05 CONTINUE c return RETURN END c c......................................................................... c DOUBLE PRECISION FUNCTION thresholdE(k,nek,s) c find the threshold E; this is due to some freaking esoteric logic c in Verner's phfti2 code (which came uncommented!). the need for c it is explained in the introduction to his 1996 paper; however, it c is not clear how it is implemented in his code c c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c implicit none include 'rates.h' integer k,nek,s,is double precision phfit2 double precision E,x include 'com-photo.com' c we need the threshold energies for each shell, so we include the c ph1 common block from Verner's PH1 block data; we also need the c inner shells so we include Verner's NINN block data double precision ph1 COMMON/ph1/ ph1(6,30,30,7) integer ninn common/ninn/ ninn(30) ! inner shell c obtain the threshold for this shell as employed by Verner is = s E = ph1(1,k,nek,is) x = phfit2(k,nek,is,E) IF (x.gt.0.0) then thresholdE = E ELSE is = ninn(nek) thresholdE = ph1(1,k,nek,is) END IF RETURN END c c......................................................................... c SUBROUTINE rate_integral(alpha) c integrate! we integrate the function FUNK from E1 to E2 c the integrand FUNK is the product of the photoionization cross c section and the mean intensity of the ionizing spectrum (SED) c we use the Numerical Recipe integrator called QROMB c c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c implicit none include 'rates.h' external funk double precision funk double precision E1,E2,alpha include 'com-photo.com' E1 = 10.0d0**logEmin E2 = 10.0d0**logEmax CALL qromb(funk,E1,E2,alpha) alpha = 4.0d0*pi*alpha RETURN END c c......................................................................... c DOUBLE PRECISION FUNCTION funk(E) c the integrand, which is the product of the photoionization cross c section and the mean intensity of the ionizing spectrum (SED) c since we have stored the cross sections and the SED in log10 units c we exponentiate; recall that the cross section and the SED are c stored in arrays with dependent variable log10(E) at intervals of c E; thus, we require interpolation to obtain their values at c arbitrary E c the function Flogsig(E) performs cubic spline interpolation in c order to return log10 of the photoionization cross section at c arbitrary energy E c the function FlogJEsed(E) performs cubic spline interpolation in c order to return log10 of the ionizing SED at arbitrary energy E c c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c implicit none include 'rates.h' double precision E,sigE,JE double precision Flogsig,FlogJEsed include 'com-photo.com' sigE = 10.0d0**Flogsig(E) JE = 10.0d0**FlogJEsed(E) funk = sigE*JE/E RETURN END c c......................................................................... c DOUBLE PRECISION FUNCTION FlogJEsed(E) c this function returns the value of log(JEsed) at arbitrary energy c E (eV) for the integrating routine by spline interpolation c this is a converted Numerical Recipe (splint.f) c c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c IMPLICIT DOUBLE PRECISION (A-H,O) include 'rates.h' include 'com-photo.com' X = log10(E) KLO=1 KHI=NEsed 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF (logEsed(K).GT.X) THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF DH = logEsed(KHI) - logEsed(KLO) IF (DH.EQ.0.) THEN WRITE(6,*) 'ERROR(FlogJE): E cannot be bracketed due to' WRITE(6,*) 'bad logEsed input at indices KLO,KHI= ',KLO,KHI WRITE(6,*) 'E(KLO),JE(KLO) = ',logEsed(KLO),logJEsed(KLO) WRITE(6,*) 'E(KHI),JE(KHI) = ',logEsed(KHI),logJEsed(KHI) STOP END IF A = (logEsed(KHI)-X)/DH B = (X-logEsed(KLO))/DH Y = A*logJEsed(KLO)+B*logJEsed(KHI) + * ((A**3-A)*d2JEsed(KLO)+(B**3-B)*d2JEsed(KHI))*(DH**2)/6.0d0 FlogJEsed = Y RETURN END c c......................................................................... c DOUBLE PRECISION FUNCTION Flogsig(E) c this function returns the value of log(sigkj) at arbitrary energy c E for the integrating routine by spline interpolation c this is a converted Numerical Recipe (splint.f) c c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c IMPLICIT DOUBLE PRECISION (A-H,O-Z) include 'rates.h' include 'com-photo.com' X = log10(E) KLO=1 KHI=NExsec 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF (logEsig(K).GT.X) THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF DH = logEsig(KHI) - logEsig(KLO) IF (DH.EQ.0.) THEN WRITE(6,*) 'ERROR(Flogsig): E cannot be bracketed due to' WRITE(6,*) 'bad logEsig input at indices KLO,KHI= ',KLO,KHI STOP END IF A = (logEsig(KHI)-X)/DH B = (X-logEsig(KLO))/DH Y = A*logsigkj(KLO)+B*logsigkj(KHI) + * ((A**3-A)*d2sigdE2(KLO)+(B**3-B)*d2sigdE2(KHI))*(DH**2)/6.0d0 Flogsig = Y RETURN END