function ionabs,w,abund=abund,ionfrac=ionfrac,vfkydir=vfkydir,ikev=ikev,$
nconsec=nconsec,DEM=DEM,logt=logt,verbose=verbose,noHeH=noHeH,$
icrstr=icrstr, _extra=e
;+
;function ionabs
; returns the photoelectric absorption cross-sections for a given
; photon energy for specified chemical composition, and ion fractions.
; Uses ground state photoionisation cross-sections computed using the
; fortran code of Verner, D.A., Ferland, G.J., Korista, K.T., & Yakovlev,
; D.G., 1996, ApJ, 465, 487 and is supplemented with high resolution
; cross-sections in the vicinity of 1s-2p resonances for
; O (from Garcia et al., 2005, ApJS, 158, 68),
; Ne (from Juett et al., 2006, ApJ, 648, 1066), and
; C (from Hasogle et al, 2010, ApJ, 724, 1296).
;
;syntax
; sigabs=ionabs(w,abund=abund,ionfrac=ionfrac,vfkydir=vfkydir,/ikeV,$
; nconsec=nconsec,DEM=DEM,logt=logt,verbose=verbose,noHeH=noHeH,$
; icrstr=icrstr,chidir=chidir,eqfile=eqfile)
;
;parameters
; w [INPUT; required] photon energy values at which to compute
; the absorption cross-section
; * 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 Grevesse & Sauval
; ionfrac [I/O] the ion fraction for all elements
; * must be a 2D array of size (NZ,NION), where
; NION=NZ+1 and usually NZ=30
; * note that in normal cases, total(IONFRAC,2)
; is an array of 1s.
; * if not given, assumed to be all neutral, IONFRAC[*,0]=1
; - unless DEM is given, see below
; * if non-zero scalar, assumed to be all neutral
; - DEM is ignored
; * if given as an array of size NZ, each element
; is assumed to have a total of that fraction,
; and all the ions have the same fraction
; - this is clearly an artificial case, to be used mainly
; for debugging and the like
; - DEM is studiously ignored
; * if given as a 1D array of size not matching NZ,
; is ignored with a warning and the default is used
; - unless DEM is given, see below
; * if given as a 2D array of size not matching (NZ,NZ+1),
; program quits with an error
; - unless DEM is given, see below
; * if DEM (see below) is given, is read in as a 3D array
; of size (NT,NZ+1,NZ1) from the database (using the PINTofALE
; function RD_IONEQ(), which calls the CHIANTI procedure
; READ_IONEQ), weights the fractions according to the DEM
; and is converted into a 2D array of size (NZ,NZ+1)
; vfkydir [INPUT] where to find the save files that contain the
; cross-sections calculated for each ionization state
; * default is '$ARDB'
; * NOTE: the cross-sections are not read in if the keyword
; ICRSTR is properly defined
; ikeV [INPUT] if set, assumes that W are in units of keV
; nconsec [INPUT] number of grid points that define the cross-sections
; that must fall within one bin of the user defined grid before
; the switch from averaging to interpolating occurs
; * default=1
; DEM [INPUT] the Differential Emission Measure that will be
; used to weight the different temperatures to compute the
; ionization fraction
; * looked at only if IONFRAC is not supplied as input
; * size must match LOGT, otherwise ignored
; logT [INPUT] the temperature grid over which DEM is defined
; * size must match DEM, otherwise DEM is ignored
; 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
; icrstr [I/O] the cross-sections calculated for each ionization
; state, stored in a structure of structures of the form
; ICRSTR.(Z) = {EVSPLAC, CROSSI}
; * if given on input, will use these numbers
; * if the structure does not contain enough tags, or if
; any one of the substructures does not contain enough
; tags, will be read in using the data in VFKYDIR
; verbose [INPUT] controls chatter
; _extra [INPUT ONLY] pass defined variables to subroutines
; -- RD_IONEQ : CHIDIR, EQFILE
;
;subroutines
; GETABUND()
; RD_IONEQ()
; READ_IONEQ
; INICON
; GETPOADEF()
;
;history
; written as IONTAU by Jeremy Drake (Apr07)
; further shoehorned into PoA format by Vinay Kashyap (Apr07)
; bug fixes and corrections and changed behavior of IONFRAC
; (VK; May09)
; updated references (JJD; Aug11)
; BUG: DEM was not being normalized when weighting IONFRAC (JJD; Nov11)
; added call to GETPOADEF (VK; Aug15)
;-
; 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='energy/wavelength grid is undefined' else $
if nw eq 1 then ok='energy/wavelength grid must be an array'
if ok ne 'ok' then begin
print,'Usage: sigabs=ionabs(w,abund=abund,ionfrac=ionfrac,vfkydir=vfkydir,$'
print,' /ikeV,nconsec=nconsec,logt=logt,verbose=verbose,noHeH=noHeH,$'
print,' icrstr=icrstr,chidir=chidir,eqfile=eqfile)'
print,' compute absorption cross-section given abundance and ion fractions'
print,' based on Verner et al. (1996) calculations'
if np ne 0 then message,ok,/informational
return,0.
endif
;message,'*** MEN AT WORK *** DO NOT USE ***'
; need this
inicon,atom=atom,aname=aname
; check keywords
;VERBOSE
vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1
;ABUND
abun=getabund('Grevesse & Sauval')
if n_elements(abund) gt 2 then abun=abund & nab=n_elements(abun)
atom=strlowcase(atom[0L:nab-1L]) & Z=indgen(nab)+1 & nZ=n_elements(Z)
;VFKYDIR
;vdir='/data/letg10/gorczyca/iontau'
;if size(vfkydir,/type) eq 7 then vdir=vfkydir[0]
;NCONSEC
mconsec=1 & if keyword_set(nconsec) then mconsec=fix(nconsec)
;IONFRAC, DEM(LOGT)
; if IONFRAC not given and DEM not given: neutral
; if IONFRAC not given and DEM is given: read in
; if IONFRAC is non-zero scalar: neutral, and always ignore DEM
; if IONFRAC is array of size NZ: flat across ions, and ignore DEM
; if IONFRAC is array of size !NZ: neutral, but read in if DEM given
; if IONFRAC is 2D array of size [NZ,NION]: use as is and ignore DEM
; if IONFRAC is 2D but ![NZ,NION]: quit with error, but read if DEM given
; anything else: read if DEM is given, quit with error otherwise
ok='neutral' & nion=n_elements(ionfrac) & szion=size(ionfrac)
aye='nay' & nDEM=n_elements(DEM) & nT=n_elements(logT)
if nDEM ne 0 and nT ne 0 and nDEM eq nT then aye='aye'
if not keyword_set(ionfrac) and aye eq 'aye' then ok='read'
if keyword_set(ionfrac) and nion eq 1 then ok='neutral'
if szion[0] eq 1 then begin
if nion eq nZ then ok='flat' else begin
if aye eq 'aye' then ok='read' else ok='neutral'
endelse
endif
if szion[0] eq 2 then begin
if szion[1] eq nZ and szion[2] eq nZ+1 then ok='ok' else begin
if aye eq 'aye' then ok='read' else ok='quit'
endelse
endif
if szion[0] gt 2 then begin
if aye eq 'aye' then ok='read' else ok='quit'
endif
;
case ok of
'ok': ifraction=ionfrac
'neutral': begin ;(default
ifraction=fltarr(nZ,nZ+1) & ifraction[*,0]=1.
end ;NEUTRAL)
'flat': begin ;(same for all ions
tmp=fltarr(nZ,nZ+1)
for iab=0,nZ-1 do tmp[i,0:i]=ionfrac[i]/(i+1.)
ifraction=tmp
end ;FLAT)
'read': begin ;(read in from database and fold in DEM
ionfrac=rd_ioneq(indgen(nZ)+1,logT, _extra=e)
tmp=ionfrac
for i=0L,nT-1L do tmp[i,*,*]=DEM[i]*ionfrac[i,*,*]/total(DEM) ;had to divide by total(DEM) to ensure proper fractional weighting
ifraction=total(tmp,1);/nT ;this division by nT seems to be a residue of previous attempts to normalize for the DEM, and so, given above correction by total(DEM), has now been commented out
; the output of RD_IONEQ is in [NT,NION,NZ],
; so the remaining axes must be transposed
ifraction=transpose(ifraction)
end ;READ)
else: begin ;(whatchewtakinaboot, willis?
message,'IONFRAC is not understandable; quitting',/informational
return,-1L
end ;QUIT)
endcase
ionfrac=ifraction
;ok='ok' & nion=n_elements(ionfrac) & szion=size(ionfrac)
;aye='flatDEM' & nDEM=n_elements(DEM) & nT=n_elements(logT)
;if nDEM eq 0 then aye='flatDEM' else $
; if nT eq 0 then aye='flatDEM' else $
; if nDEM ne nT then begin
; if nDEM eq 81 or nDEM eq 41 then aye='specialDEM'
; endif else aye='DEM'
;case szion[0] of
; 0: begin
; if nDEM gt 0 then ok='read' else ok='flat'
; mZ=nab
; end
; 1: begin
; if szion[1] eq nab then ok='flatZ' else $
; if szion[1] eq nab+1 then ok='flatION' else ok='flat'
; end
; 2: begin
; mZ=szion[1] & mION=szion[2]
; if mZ ne nab or mION ne nab+1 then ok='read'
; end
; 3: begin
; mT=szion[1] & mZ=szion[2] & mION=szion[3]
; if aye eq 'DEM' or aye eq 'specialDEM' then begin
; if mZ eq nab and mION eq nab+1 then ok='DEM' else $
; ok='read'
; endif
; end
; else: ok='read'
;endcase
;;
;ifraction=fltarr(nab,nab+1) & ifraction[*,0]=1.
;if ok eq 'ok' then ifraction=ionfrac else $
; if ok eq 'read' then ifraction=rd_ioneq(indgen(nZ)+1,logT, _extra=e)
; if ok eq 'flatZ' then for i=0L,nab-1L do ifraction[i,*]=ionfrac[i] else $
; if ok eq 'flatION' then for i=0L,nab do ifraction[*,i]=ionfrac[i] else $
; if ok eq 'DEM' or ok eq 'DEM' then begin
; xDEM=DEM
; if aye eq 'specialDEM' then begin
; if nDEM ne mT then begin
; xlogT=findgen(41)*0.1+4.
; if nDEM eq 81 then xlogT=findgen(81)*0.05+4.
; xDEM=interpol(DEM,xlogT,findgen(mT)*4./float(mt-1L)+4.)
; endif
; endif else xlogT=logT
; tmp=ionfrac
; for i=0L,mT-1L do tmp[i,*,*]=xDEM[i]*ionfrac[i,*,*]
; ifraction=total(tmp,1)/mT
; ionfrac=ifraction
; endif
;ionfrac=ifraction
;IKEV
if not keyword_set(ikev) then begin
ev=reverse(12398.5/w) ;input is in [Ang] -- convert to [eV]
endif else begin
ev=w*1e3 ;input is in [keV] -- convert to [eV]
endelse
;ICRSTR
nICR=n_tags(icrstr)
; cross-sections are stored in save files in the form of total
; cross-section for each ion for each element with Z=1-30
; energy grid is highly non-uniform so as to capture edge energies
; and structure accurately.
; structure of each array is (nion,nenergy)
;if not keyword_set(vfkydir) then vfkydir='/data/letg10/gorczyca/iontau/'
;zTOPDIR='/data/fubar/SCAR/'
;ivar=0 & defsysv,'!TOPDIR',exists=ivar ;if !TOPDIR exists
;if ivar ne 0 then setsysval,'TOPDIR',zTOPDIR,/getval else begin
; tmp=routine_info('IONABS',/source,/functions) & scardir=tmp.PATH
; jvar=strpos(scardir,'pro/ionabs.pro') ;we know where IONABS.PRO is
; zTOPDIR=strmid(scardir,0,jvar-1)
; if vv ge 10 then message,'PoA directory is: '+zTOPDIR,/info
;endelse
;zARDB=filepath('',root_dir=filepath('ardb',root_dir=zTOPDIR))
zTOPDIR=getpoadef()
zARDB=getpoadef('ARDB')
;if not keyword_set(vfkydir) then vfkydir=zARDB
if size(vfkydir,/type) ne 7 then vdir=zARDB else vdir=vfkydir[0]
elem=strlowcase(atom[0:nab-1])
elemnam=aname[0:nab-1]
; set up ev array bin boundaries and the binsize for rebin when user
; grid is more coarse than save file data (needed for conservation of
; cross-section over fine resonances)
dev=ev[1:*]-ev & dev=[dev[0],dev]
bev=[ev[0]-dev[0]/2.,ev+dev/2.]
; set up interpolation and rebin array to contain the cumulative cross-section
totcrossrb=fltarr(n_elements(ev))*0.
if vv gt 10000L then stop,'HALTing; type .CON to continue'
; main loop to read in data for each element, then loop over each
; element to sum up the cross-sections according to weighting
; by ion fractions. Weight total element cross-section by abundance,
; where abundance of H=1
for i=0,n_elements(elem)-1 do begin
splacfil=filepath(elem[i]+'splac.cross',root_dir=vdir)
if nICR eq 0 then begin
restore,splacfil
tmp=create_struct(elemnam[i],create_struct('EVSPLAC',EVSPLAC,'CROSSI',CROSSI))
if i eq 0 then icrstr=tmp else icrstr=create_struct(icrstr,tmp)
if vv gt 0 then print,splacfil
if vv gt 100000L then stop,'HALTing; type .CON to continue'
endif else begin
evsplac=icrstr.(i).EVSPLAC & crossi=icrstr.(i).CROSSI
endelse
; compute total cross-section:
; need total(ionfrac*cross-secton)*abun for each element
; then rebin onto supplied grid
; In order to sample properly when original data grid is finer than
; user requested, we need rebin
; To deal with regions in which user grid is finer than data grid,
; need to use interpolation.
; Here, both methods are used and then combined according to which of
; these regimes the grids are in.
totcross=abun[i]*reform(ifraction[i,0:i]#crossi)
totcrossi=interpol(totcross,evsplac,ev)
devsplac=evsplac[1:*]-evsplac & devsplac=[devsplac[0],devsplac]
devsplaci=interpol(devsplac,evsplac,ev)
oo=where(devsplaci lt dev,count)
if (count gt mconsec) then begin
doo=oo[1:*]-oo & doo=smooth([doo[0],doo],mconsec)
oo=where(doo lt 1.5,count)
if (count gt mconsec) then totcrossi[oo]=(rebinw(totcross,evsplac,bev))[oo]
endif
totcrossrb=totcrossrb+totcrossi
endfor
; convert to optical depth using column density. Recall that units
; from Verner routines are Mb (10^-18 cm^2)
;totcrossrb=nh*totcrossrb*1.e-18
totcrossrb=totcrossrb*1.e-18
if not keyword_set(ikev) then totcrossrb=reverse(totcrossrb)
return,totcrossrb
end