function bamabs,w,abund=abund,ikeV=ikeV,Fano=Fano,noHeH=noHeH,$ verbose=verbose, _extra=e ;+ ;function bamabs ; return photoelectric absorption cross-sections in energy range ; 30 eV - 10 keV using the polynomial fit coefficients determined ; by Monika Balucinska-Church and Dan McCammon (Balucinska-Church, ; M.\ and McCammon, D.\ 1992, ApJ 400, 699). This is an update ; to Morrison & McCammon 1983. ; ;syntax ; sigabs=bamabs(w,abund=abund,/ikeV,/Fano,verbose=v) ; ;parameters ; w [INPUT; required] photon energy values at which to compute ; the absorption ; * default units are [Ang], unless IKEV is set, when they are ; assumed to be [keV] ; ;keywords ; abund [INPUT] abundances relative to H=1 ; * default is to use Anders & Grevesse 1982 ; ikeV [INPUT] if set, assumes that W are in units of keV ; Fano [INPUT] if set, the 4 strongest auto-ionizing resonances ; of He I are included; numbers come from Oza (1986, Phys ; Rev A, 33, 824), and Fernley et al. (1987, J.Phys B 20, 6457). ; noHeH [INPUT] if set to any number other than 1 or 2, excludes ; H and He from the cross-sections ; * if set to 1, only excludes H ; * if set to 2, only excludes He ; verbose [INPUT] controls chatter ; _extra [JUNK] here only to avoid crashing the program ; ;subroutines ; GETABUND ; INICON ; ;restrictions ; works only in energy range 30 eV - 10 keV ; be warned that the edges are idealized and thus unrealistic ; for high resolution spectra because they do not include ; any resonance structure ; ;history ; vinay kashyap (OctMM; based on TOTLXS.FOR and XSCTNS.FOR of ; Balucinska-Church & McCammon, obtained from ADC catalog VI/62A, ; ftp://adc.gsfc.nasa.gov/pub/adc/archives/catalogs/6/6062A/ ) ; bug corrections (Brian Kern, 5Dec2001): Ne cross-section was ; missing "^3"; Fano calc had EPSj[j] (VK; Dec2001) ; force calculations only in valid energy range; added keyword ; NOHEH (VK; Jun03) ; bug correction: "exclude=0" means "include" (VK; Jul03) ; bug correction: lambda selection for He (Andy Ptak; VK Jun05) ; edge lambda correction: O, N, C changed acc. to numbers in ; Atomic Data Booklet (2nd Edition), Thompson, A., et al., 2001, ; LBL/PUB-490 Rev. 2 (Jeremy Drake; VK Sep2006) ;- ; usage ok='ok' & np=n_params() & nw=n_elements(w) if np eq 0 then ok='Insufficient parameters' else $ if nw eq 0 then ok='W is undefined' if ok ne 'ok' then begin print,'Usage: sigabs=bamabs(w,abund=abund,/ikeV,/Fano,noHeH=noHeH,verbose=v)' print,' return photoelectric absorption cross-sections following' print,' Balucinska-Church & McCammon 1992 update to Morrison & McCammon 1983' if np ne 0 then message,ok,/info return,-1L endif ; inputs nrg=[float(abs(w))] nab=n_elements(abund) if nab lt 30 then abund=getabund('anders & grevesse') if n_elements(ikeV) eq 0 then ikeV=0 if not keyword_set(ikeV) then begin ;convert from [Ang] to [keV] oo=where(nrg gt 0,moo) if moo gt 0 then nrg[oo]=12.3985/nrg[oo] else return,0*w endif v=0 & if keyword_set(verbose) then v=fix(abs(verbose[0])) > 0 exclH=0 & exclHe=0 if keyword_set(noHeH) then begin exclH=1 & exclHe=1 if noHeH[0] eq 1 then exclHe=0 ;only exclude H if noHeH[0] eq 2 then exclH=0 ;only exclude He endif ; initialize inicon,amu=amu,atom=atom,fundae=f ;H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Cl, Ar, Ca, Cr, Fe, Ni useZ=[1, 2, 6, 7, 8, 10, 11, 12, 13, 14,16, 17, 18, 20, 24, 26, 28] nZ=n_elements(useZ) EE = nrg*1e3 ;keV -> eV ;ok=where(EE gt 0 and EE le 1e4,mok) & elog=fltarr(nw)-90. ok=where(EE ge 30 and EE le 1e4,mok) & elog=fltarr(nw)-90. if mok gt 0 then elog[ok]=alog(EE[ok]) ; output sig=0.*nrg ;photoelectric absorption cross-section in [cm^2] for i=0L,nZ-1L do begin ;{step through encoded elements iZ=useZ[i] & aw=(amu.(iZ-1L))[0] & elem=atom[iZ-1L] case elem of 'H': begin ;HYDRO (ANGIE BETKER, JEFF BLOCH, MBC 1991) hydro=fltarr(nw) if mok gt 0 then begin x=21.46941 + (3. - 2.060152)*Elog - 0.1492932*Elog*Elog +$ 5.4634294e-3*(Elog^3) hydro[ok]=Exp(X[ok])/(EE[ok]^3) endif ; sig=sig + aw*hydro/f.NA*abund[iZ-1L]*(1-exclH) end 'He': begin ;HELIUM (MBC 1997) helium=fltarr(nw) C1=[-4.7416, 14.8200, -30.8678, 37.3584, -23.4585, 5.9133] ;polynomial coefficients for Yan et al data IP=n_elements(C1) Q=[2.81, 2.51, 2.45, 2.44] ;parameters Q for resonances (Fernley et al. 1987) EION=24.58 ;ionization edge in eV NU=[1.610, 2.795, 3.817, 4.824] ;parameters NU for resonances (Oza 1986) GAMMA=[2.64061E-3, 6.20116E-4, 2.56061E-4, 1.320159E-4] ;parameters GAMMA for resonances (Oza 1986) ; lambda=fltarr(nw)+12398.54 if mok gt 0 then begin ;(non-zero energies lambda[ok]=12398.54/EE[ok] X=EE/EION oo=where(lambda ge 1.239854 and lambda le 503.97,moo) ; was "or", which, as Andy Ptak points out, is always true y=fltarr(nw)+1. & sig_yan=fltarr(nw) if moo gt 0 then begin for j=1L,ip do y[oo]=y[oo]+c1[j-1L]/(x[oo]^(j/2.)) sig_yan[oo]=733.E-24/(EE[oo]/1e3)^(3.5) * Y[oo] if keyword_set(Fano) then begin ;(Fano profile correction ;FANO (MBC 1993) nf=n_elements(Q) for j=0L,nf-1L do begin eps = 911.2671/LAMBDA epsj = 3.0 - 1./(NU[j]^2) + 1.807317 XX = 2.0*(EPS-EPSj)/GAMMA[j] fano_factor=(XX-Q[j])^2/(1.+XX^2) sig_yan = sig_yan * fano_factor endfor endif ;Fano) endif helium = sig_yan * f.NA/aw endif ;MOK>0) ; sig = sig + aw*helium/f.NA*abund[iz-1L]*(1-exclHe) end 'C': begin ;CARBON carbon=fltarr(nw) if mok gt 0 then begin ;(non-zero energies ;o1=where(EE gt 0 and EE lt 284.,mo1) ;o2=where(EE ge 284. and EE le 1e4,mo2) o1=where(EE gt 0 and EE lt 284.2,mo1) ;JJD: cf. ADB rev2 o2=where(EE ge 284.2 and EE le 1e4,mo2) ;JJD: cf. ADB rev2 X=fltarr(nw)-200. if mo1 gt 0 then X[o1] =$ 8.74161 + (7.13348*Elog[o1]) + (-1.14604*Elog[o1]^2) +$ (0.0677044*Elog[o1]^3) if mo2 gt 0 then X[o2] =$ 3.81334 + (8.93626*Elog[o2]) + (-1.06905*Elog[o2]^2) +$ (0.0422195*Elog[o2]^3) carbon[ok]=exp(x[ok])/(ee[ok]^3) endif ;MOK>0) ; sig = sig + aw*carbon/f.NA*abund[iz-1L] end 'N': begin ;NITRO nitro=fltarr(nw) if mok gt 0 then begin ;(non-zero energies ;o1=where(EE gt 0 and EE lt 401.,mo1) ;o2=where(EE ge 401. and EE le 1e4,mo2) o1=where(EE gt 0 and EE lt 409.9,mo1) ;JJD: cf. ADB rev2 o2=where(EE ge 409.9 and EE le 1e4,mo2) ;JJD: cf. ADB rev2 X=fltarr(nw)-200. if mo1 gt 0 then X[o1] =$ 9.24058 + (7.02985*Elog[o1]) + (-1.08849*Elog[o1]^2) +$ (0.0611007*Elog[o1]^3) if mo2 gt 0 then X[o2] =$ -13.0353 + (15.4851*Elog[o2]) + (-1.89502*Elog[o2]^2) +$ (0.0769412*Elog[o2]^3) ; nitro[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*nitro/f.NA*abund[iz-1L] end 'O': begin ;OXYGEN oxygen=fltarr(nw) if mok gt 0 then begin ;(non-zero energies ;o1=where(EE gt 0 and EE lt 531.7,mo1) ;o2=where(EE ge 531.7 and EE le 1e4,mo2) o1=where(EE gt 0 and EE lt 543.1,mo1) ;JJD: cf. ADB rev2 o2=where(EE ge 543.1 and EE le 1e4,mo2) ;JJD: cf. ADB rev2 X=fltarr(nw)-200. if mo1 gt 0 then X[o1] =$ 2.57264 + (10.9321*Elog[o1]) + (-1.79383*Elog[o1]^2) +$ (0.102619*Elog[o1]^3) if mo2 gt 0 then X[o2] =$ 16.53869 + (0.6428144 + 3.)*Elog[o2] - 0.3177744*Elog[o2]^2 +$ 7.9471897e-3*(Elog[o2]^3) ; oxygen[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*oxygen/f.NA*abund[iz-1L] end 'Ne': begin ;NEON neon=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 867.,mo1) o2=where(EE ge 867. and EE le 1e4,mo2) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] =$ -3.04041 + (13.0071*Elog[o1]) + (-1.93205*Elog[o1]^2) +$ (0.0977639*Elog[o1]^3) if mo2 gt 0 then X[o2] =$ 17.6007 + (3.29278*Elog[o2]) + (-0.263065*Elog[o2]^2) +$ (5.68290E-3*Elog[o2]^3) ; neon[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*neon/f.NA*abund[iz-1L] end 'Na': begin ;SODIUM sodium=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 1071.7,mo1) o2=where(EE ge 1071.7 and EE le 1e4,mo2) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] =$ -2737.598 + (2798.704 + 3.)*Elog[o1] - 1009.892*Elog[o1]^2 +$ 87.16455*(Elog[o1]^3) + 43.20644*(Elog[o1]^4) -$ 15.27259*(Elog[o1]^5) + 2.180531*(Elog[o1]^6) -$ 0.1526546*(Elog[o1]^7) + 4.3137977E-3*(Elog[o1]^8) if mo2 gt 0 then X[o2] = $ 1.534019 + (6.261744 + 3.)*Elog[o2] - 0.9914126*Elog[o2]^2 +$ 3.5278253E-2*(Elog[o2]^3) ; sodium[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*sodium/f.NA*abund[iz-1L] end 'Mg': begin ;MAGNES magnes=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 49.45,mo1) o2=where(EE ge 49.45 and EE lt 1303.4,mo2) o3=where(EE ge 1303.4 and EE le 1e4,mo3) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] = 7.107172 + (0.7359418 + 3.)*Elog[o1] if mo2 gt 0 then X[o2] = $ -81.32915 + (62.2775 + 3.)*Elog[o2] - 15.00826*Elog[o2]^2 +$ 1.558686*(Elog[o2]^3) - 6.1339621E-2*(Elog[o2]^4) if mo3 gt 0 then X[o3] = $ -9.161526 + (10.07448 + 3.)*Elog[o3] - 1.435878*Elog[o3]^2 +$ 5.2728362E-2*(Elog[o3]^3) ; magnes[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*magnes/f.NA*abund[iz-1L] end 'Al': begin ;ALUM alum=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 72.78,mo1) o2=where(EE ge 72.78 and EE lt 1559.9,mo2) o3=where(EE ge 1559.9 and EE le 1e4,mo3) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] =$ 26.90487 + (3. - 9.135221)*Elog[o1] + 1.175546*Elog[o1]^2 if mo2 gt 0 then X[o2] = -38.1232 + 29.5161*Elog[o2] -$ 4.45416*Elog[o2]^2 + 0.226204*Elog[o2]^3 if mo3 gt 0 then X[o3] =$ 14.6897 + 4.22743*Elog[o3] - 0.344185*Elog[o3]^2 +$ 8.18542E-3*Elog[o3]^3 ; alum[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*alum/f.NA*abund[iz-1L] end 'Si': begin ;SILICN silicn=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 100.6,mo1) o2=where(EE ge 100.6 and EE lt 1840.0,mo2) o3=where(EE ge 1840.0 and EE le 1e4,mo3) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] = $ -3.066295 + (7.006248 + 3.)*Elog[o1] - 0.9627411*Elog[o1]^2 if mo2 gt 0 then X[o2] = $ -182.7217 + (125.061 + 3.)*Elog[o2] - 29.47269*Elog[o2]^2 +$ 3.03284*(Elog[o2]^3) - 0.1173096*(Elog[o2]^4) if mo3 gt 0 then X[o3] = $ -33.39074 + (18.42992 + 3.)*Elog[o3] -$ 2.385117*Elog[o3]^2 + 8.887583e-2*(Elog[o3]^3) ; silicn[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*silicn/f.NA*abund[iz-1L] end 'S': begin ;SULFUR sulfur=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 165.0,mo1) o2=where(EE ge 165.0 and EE lt 2470.5,mo2) o3=where(EE ge 2470.5 and EE le 1e4,mo3) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] = $ 598.2911 + (3. - 678.2265)*Elog[o1] + 308.1133*Elog[o1]^2 -$ 68.99324*(Elog[o1]^3) + 7.62458*(Elog[o1]^4) - 0.3335031*(Elog[o1]^5) if mo2 gt 0 then X[o2] = $ 3994.831 + (3. - 3693.886)*Elog[o2] + 1417.287*Elog[o2]^2 -$ 287.9909*(Elog[o2]^3) + 32.70061*(Elog[o2]^4) -$ 1.968987*(Elog[o2]^5) + 4.9149349E-2*(Elog[o2]^6) if mo3 gt 0 then X[o3] = $ -22.49628 + (14.24599+ 3.)*Elog[o3] - 1.848444*Elog[o3]^2 +$ 6.6506132E-2*(Elog[o3]^3) ; sulfur[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*sulfur/f.NA*abund[iz-1L] end 'Cl': begin ;CHLORN chlorn=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 202.0,mo1) o2=where(EE ge 202.0 and EE lt 2819.6,mo2) o3=where(EE ge 2819.6 and EE le 1e4,mo3) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] = $ 6253.247 +(3. - 8225.248)*Elog[o1] + 4491.675*Elog[o1]^2 -$ 1302.145*(Elog[o1]^3) + 211.4881*(Elog[o1]^4) -$ 18.25547*(Elog[o1]^5) + 0.6545154*(Elog[o1]^6) if mo2 gt 0 then X[o2] = $ -233.0502 + (143.9776 + 3.)*Elog[o2] - 31.12463*Elog[o2]^2 +$ 2.938618*(Elog[o2]^3) - 0.104096*(Elog[o2]^4) if mo3 gt 0 then X[o3] = $ -23.74675 + (14.50997 + 3.)*Elog[o3] - 1.857953*Elog[o3]^2 +$ 6.6208832E-2*(Elog[o3]^3) ; chlorn[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*chlorn/f.NA*abund[iz-1L] end 'Ar': begin ;ARGON argon=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 245.0,mo1) o2=where(EE ge 245.0 and EE lt 3202.9,mo2) o3=where(EE ge 3202.9 and EE le 1e4,mo3) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] = $ -330.3509 + (267.7433 + 3.)*Elog[o1] - 78.90498*Elog[o1]^2 +$ 10.35983*(Elog[o1]^3) - 0.5140201*(Elog[o1]^4) if mo2 gt 0 then X[o2] = $ -5.71870 + (8.85812*Elog[o2]) + (-0.307357*Elog[o2]^2) +$ (0.00169351*(Elog[o2]^3)) + (-0.0138134*(Elog[o2]^4)) +$ (0.00120451*(Elog[o2]^5)) if mo3 gt 0 then X[o3] = $ 19.1905 + (2.74276*Elog[o3]) + (-0.164603*Elog[o3]^2) +$ (0.00165895*Elog[o3]^3) ; argon[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*argon/f.NA*abund[iz-1L] end 'Ca': begin ;CALC calci=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 349.31,mo1) o2=where(EE ge 349.31 and EE lt 4038.1,mo2) o3=where(EE ge 4038.1 and EE le 1e4,mo3) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] = $ -873.972 + (865.5231 + 3.)*Elog[o1] - 339.678*Elog[o1]^2 +$ 66.83369*(Elog[o1]^3) - 6.590398*(Elog[o1]^4) + 0.2601044*(Elog[o1]^5) if mo2 gt 0 then X[o2] = $ -3449.707 + (2433.409 + 3.)*Elog[o2] - 682.0668*Elog[o2]^2 +$ 95.3563*(Elog[o2]^3) - 6.655018*(Elog[o2]^4) + 0.1854492*(Elog[o2]^5) if mo3 gt 0 then X[o3] = $ 18.89376 + (3. - 0.2903538)*Elog[o3] - 0.1377201*Elog[o3]^2 ; calci[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*calci/f.NA*abund[iz-1L] end 'Cr': begin ;CHROM chrom=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 598.0,mo1) o2=where(EE ge 598.0 and EE lt 691.0,mo2) o3=where(EE ge 691.0 and EE lt 5988.8,mo3) o4=where(EE ge 5988.8 and EE le 1e4,mo4) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] = $ -0.4919405 + (12.66939 + 3.)*Elog[o1] - 5.199775*Elog[o1]^2 +$ 1.086566*(Elog[o1]^3) - 0.1196001*(Elog[o1]^4) +$ 5.2152011E-3*(Elog[o1]^5) if mo2 gt 0 then X[o2] = 27.29282 +(3. - 2.703336)*Elog[o2] if mo3 gt 0 then X[o3] = -15.2525 + (13.23729 + 3.)*Elog[o3] -$ 1.966778*Elog[o3]^2 + 8.062207E-2*(Elog[o3]^3) if mo4 gt 0 then X[o4] = 8.307041 + (2.008987 + 3.)*Elog[o4] -$ 0.2580816*Elog[o4]^2 ; chrom[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*chrom/f.NA*abund[iz-1L] end 'Fe': begin ;IRON iron=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 707.4,mo1) o2=where(EE ge 707.4 and EE lt 7111.2,mo2) o3=where(EE ge 7111.2 and EE le 1e4,mo3) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] = $ -15.07332 + (18.94335 + 3.)*Elog[o1] - 4.862457*Elog[o1]*Elog[o1] +$ 0.5573765*(Elog[o1]^3) - 3.0065542E-2*(Elog[o1]^4) +$ 4.9834867E-4*(Elog[o1]^5) if mo2 gt 0 then X[o2] = $ -253.0979 + (135.4238 + 3.)*Elog[o2] - 25.47119*Elog[o2]*Elog[o2] +$ 2.08867*(Elog[o2]^3) - 6.4264648E-2*(Elog[o2]^4) if mo3 gt 0 then X[o3] = $ -1.037655 + (4.022304 + 3.)*Elog[o3] - 0.3638919*Elog[o3]*Elog[o3] ; iron[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*iron/f.NA*abund[iz-1L] end 'Ni': begin ;NICKEL nickel=fltarr(nw) if mok gt 0 then begin ;(non-zero energies o1=where(EE gt 0 and EE lt 853.6,mo1) o2=where(EE ge 853.6 and EE lt 8331.6,mo2) o3=where(EE ge 8331.6 and EE le 1e4,mo3) X=fltarr(nw)-200. if mo1 gt 0 then X[o1] = $ -7.919931 + (11.06475 + 3.)*Elog[o1] - 1.935318*Elog[o1]^2 +$ 9.3929626e-2*(Elog[o1]^3) if mo2 gt 0 then X[o2] = $ 3.71129 + (8.45098*Elog[o2]) + (-0.896656*Elog[o2]^2) +$ (0.0324889*Elog[o2]^3) if mo3 gt 0 then X[o3] = 28.4989 + (0.485797*Elog[o3]) ; nickel[ok] = EXP(X[ok])/(EE[ok]^3) endif ;MOK>0) ; sig = sig + aw*nickel/f.NA*abund[iz-1L] end else: message,elem+': cross-sections not available',/info endcase if v gt 0 then begin if nw eq 1 then print,elem+': ',sig[0],' @'+strtrim(nrg[0],2) if v gt 10 and mok gt 1 then plot,w[ok],sig[ok],/yl,title=elem,psym=3 endif endfor ;i=0,nZ-1} return,sig end