function fold_ioneq,emis,z,ion,logt=logt,tmax=tmax,trng=trng,level=level,$
userarr=userarr,eqfile=eqfile,chifil=chifil,verbose=verbose, _extra=e
;+
;function fold_ioneq
; folds in ionization equilibrium fractions of the various ions
; of a given element with line emissivities (which include only
; level population info; see RD_LINE()). returns
; emis(T;Z,ion)*ioneq(T;Z,ion)
;
;syntax
; lflx=fold_ioneq(emis,Z,ion,logt=logt,tmax=tmax,trng=trng,$
; level=level,userarr=userarr,chifil=chifil,chidir=chidir,$
; eqfile=eqfile,verbose=v)
;
;parameters
; emis [INPUT; required] 1- or 2-D array of line emissivities as
; a function of temperature
; * if 2D, assumed to be EMIS(Temperature,Wavelength)
; * if 1D, assumed to be EMIS(Temperature), unless the
; number of elements matches that of Z
; Z [INPUT; required] atomic number(s) corresponding to
; each wavelength
; * if size does not match the 2nd dimension of EMIS, then
; if >, ignore extra elements
; if <, but >0, use available ones, and use Z[0] for rest
; ion [INPUT; default: Z+1] ionic state(s) corresponding to
; each wavelength
; * if size does not match the 2nd dimension of EMIS, then
; if >, ignore extra elements
; if <, but >0, use available ones, and use ION[0] for rest
; * NOTE: this is the ionic state >>that matters<<, i.e.,
; the ionic state of the species that populates the
; upper level. See keyword JON of RD_LINE, which is
; what is needed here.
;
;keywords
; logt [INPUT] 1D array of log10(temperature [K]) at which EMIS
; is given. (default: findgen(81)*0.05+4)
; * if size does not match that of EMIS, use default
; level [INPUT; default: 0.5] level at which to determine the
; temperature range on line flux (in other words, default
; will be to return the FWHM)
; * if <1, TRNG=LOGT(where(FLUX>LEVEL*MAX(FLUX))
; * if >=1, TRNG=LOGT(where(FLUX>(LEVEL/100.)*MAX(FLUX))
; * if <=0 or >=100, set to 0.5
; tmax [OUTPUT] log10(temperature) at which contribution is maximum
; trng [OUTPUT] range in LOGT (LEVEL of MAX)
; -- TRNG(0,*) is lower bound, TRNG(1,*) is upper bound
; userarr [INPUT] an array of ionization fractions in the same
; format as is returned by RD_IONEQ(), and has the size
; [n(T),max(Z)+1,n(Z)]
; * if this is given and is legit (i.e., matches the
; supplied EMIS and LOGT), then EQFILE and CHIFIL are
; ignored
; eqfile [I/O] file from which to read in ion-balance data
; * default is to read in the CHIANTI file, !IONEQF
; * hard-coded default is
; ioneq/chianti.ioneq
; * if not given, but CHIFIL is given, uses that.
; chifil [I/O] if set, reads in from CHIANTI-type database; must be
; set to the name of the file containing the ionization
; equilibrium info (i.e., CHIDIR/CHIFIL)
; * may be overridden with EQFILE
; * default is ioneq/chianti.ioneq
; * currently, set automatically if not given.
; verbose [INPUT] higher the number, greater the chatter
; _extra [INPUT] use to transmit defined keywords to called functions
; * RD_IONEQ [CHIDIR]
;
;subroutines
; -- WHEE
; -- SETSYSVAL
; -- RD_IONEQ [READ_IONEQ (%)]
; (%) CHIANTI subroutine, used as is
;
;history
; vinay kashyap (Dec96)
; removed call to KILROY and added call to WHEE, hacked default
; use of READ_IONEQ until more options become available; corrected
; bug that was reading in ionstate-1 rather than ionstate; added
; code to handle Z=0 (VK; Feb97)
; bug: if Z is scalar take em all to be same element (VK; Apr97)
; added keyword VERBOSE (VK; MarMM)
; streamlined meanings of EQFILE, CHIFIL (VK; DecMM)
; changed default for EQFILE; now pass VERBOSE to RD_IONEQ; added
; call to SETSYSVAL (VK; Jul01)
; allowed EQFILE to be "none" (VK; Feb04)
; changed default behavior of Z and ION when sizes don't match EMIS;
; now uses first element of input if given (VK; Nov04)
; bug correction, ION is offset by 1 from array index (LL/VK; Dec04)
; added keyword USERARR, slightly modified behavior when EMIS
; was input as 1D (VK; Aug08)
; changed default of EQFILE and CHIFIL to ioneq/chianti.ioneq (VK; Aug13)
;-
; usage
ok='ok' & np=n_params() & nem=n_elements(emis) & szl=size(emis)
nz=n_elements(z) & ni=n_elements(ion)
if np lt 2 then ok='Insufficient parameters' else $
if nem eq 0 then ok='EMIS is not defined' else $
if nz eq 0 then ok='Z is not defined' else $
if szl(0) gt 2 then ok='EMIS is not in recognizable format'
if ok ne 'ok' then begin
print,'Usage: lflx=fold_ioneq(emissivity,Z,ion,logt=logt,tmax=tmax,trng=trng,$'
print,' level=level,userarr=userarr,eqfile=eqfile,chifil=chifil,$'
print,' chidir=chidir,verbose=v)'
print,' fold line emissivities read in with RD_LINE with ion balance'
print,' also accepts defined keywords CHIDIR and CHIFIL for RD_IONEQ'
if np ne 0 then message,ok,/informational
return,-1L
endif
; check input
if nz eq 0 then begin & z=1 & nz=1 & endif ;default element is H
if ni eq 0 then begin & ion=z+1 & ni=nz & endif ;default ionic state is Z+1
case szl(0) of
1: begin ;emis=emis(LOGT) or emis(WVL)
if szl(1) eq nz then begin
nt=1 & nw=nz & atno=[z(*)] & jon=[ion(*)]-1 ;the "-1" coz of IDL
endif else begin
nt=szl(1) & nw=1 & atno=[z(0)] & jon=[ion(0)]-1 ;the "-1" coz of IDL
endelse
end
2: begin ;emis=emis(LOGT,WVL)
nt=szl(1) & nw=szl(2) & atno=intarr(nw)+1 & jon=atno
if nz gt 0 then atno(*)=z(0)
if ni gt 0 then jon(*)=ion(0)
if nz gt 0 and nz le nw then atno(0L:nz-1L)=z
if ni gt 0 and ni le nw then jon(0L:nz-1L)=ion-1L
;if nz le nw then atno(0:nz-1)=([z])(*) else atno(*)=z(0:nw-1)
;if ni le nw then jon(0:ni-1)=([ion])(*)-1 else jon(*)=ion(0:nw-1)-1
end
else: return,-1L ;nolle comprendo
endcase
;
if not keyword_set(logt) then begin
tlog=findgen(nt)*(4./((nt>2)-1.))+4.
endif else tlog=logt
if n_elements(tlog) ne nt then begin
message,'Line emissivity does not match temperature grid',/info
tlog=findgen(nt)*(4./((nt>2)-1.))+4.
endif
;
if not keyword_set(level) then lev=0.5 else lev=level
if lev le 0. or lev ge 100. then lev=0.5
if lev ge 1. then lev=lev/100.
;
; HACK ALERT!!! until there are other ways to compute the
; ion balance, use CHIANTI!!! to avoid using CHIANTI, explicitly
; set CHIFIL=0
if n_elements(chifil) eq 0 then chifil=1
if keyword_set(userarr) then begin
ok='ok' & szu=size(userarr)
if szu(0) lt 2 then ok='USERARR must be 3D (Temp,Z,ION)' else $
if szu(1) ne nt then ok='USERARR does not match LOGT' else $
if szu(2) ne nz then ok='USERARR does not match Z' else $
if szu(3) ne max(Z)+1 then ok='USERARR does not match ION'
if ok eq 'ok' then useuserarr=1 else $ ;no need to call CHIANTI!
useuserarr=0
endif
;
zEQFILE='ioneq/chianti.ioneq'
ivar=0 & defsysv,'!IONEQF',exists=ivar ;if !IONEQF exists
if ivar ne 0 then setsysval,'IONEQF',zEQFILE,/getval
if not keyword_set(eqfile) then eqfile=zEQFILE
;def_eqfil='ioneq/arnaud_rothenflug.ioneq'
sze=size(eqfile) & nsze=n_elements(sze)
szc=size(chifil) & nszc=n_elements(szc)
if not keyword_set(eqfile) or sze(nsze-2) ne 7 then begin ;no EQFILE
if szc(nszc-2) eq 7 then eqfile=chifil else eqfile=zEQFILE ;def_eqfil
endif
; chatter
vv=0 & if keyword_set(verbose) then vv=long(verbose(0)) > 1
; initialize
lflx=0*emis & tmax=fltarr(nw) & trng=fltarr(2,nw)
; find all the unique elements
zz=atno(uniq(atno,sort(atno))) & mz=n_elements(zz)
; for each element, read in ionization equilibrium and multiply
for iz=0L,mz-1L do begin
oz=where(atno eq zz(iz),noz)
if zz(iz) le 0 then begin ;{error! error!
ioneq=fltarr(nt)+1.
if noz gt 0 then jon(oz)=0
endif else begin ;}{all OK
ioneq=fltarr(nt,zz(iz)+1)+1. ;default -- multiply by 1
;
whee,spoke=spoke,/moveit ;we're working!
if keyword_set(useuserarr) then ioneq=userarr(*,*,iz) else begin
if keyword_set(chifil) then begin ;use CHIANTI
szc=size(chifil) & nszc=n_elements(szc)
if szc(nszc-2) eq 7 then eqfile=chifil else begin
if not keyword_set(eqfile) then eqfile=zEQFILE ;'ioneq/chianti.ioneq'
endelse
;eqfile=chifil
;if szc(nszc-2) ne 7 then eqfile=zEQFILE ;'ioneq/chianti.ioneq'
if strlowcase(eqfile) ne 'none' then begin
if vv gt 0 then message,'Using ion balances from: '+eqfile,/info
ioneq=rd_ioneq(zz(iz),tlog,eqfile=eqfile,verbose=vv,_extra=e)
endif
endif
endelse
endelse ;}
;
for ii=0L,noz-1L do begin
if vv gt 1 then whee,spoke=spoke,/moveit ;we're working!
; fold in ionization equilibrium fractions
lflx(*,oz(ii))=emis(*,oz(ii))*ioneq(*,jon(oz(ii)))
; and here figure out TMAX ...
tmp=reform(lflx(*,oz(ii))) ;& imax=where(tmp eq max(tmp),moo)
tmpmax=max(tmp,imax)
tmax(oz(ii))=total(tlog(imax));/moo
; ... and TRNG
oo=where(tmp gt lev*tmpmax,moo)
;oo=where(tmp gt lev*tmp(imax(0)),moo)
if moo ne 0 then begin
trng(0,oz(ii))=tlog(oo(0)) & trng(1,oz(ii))=tlog(oo(moo-1))
endif
endfor
if vv gt 1 then whee,spoke=spoke,/moveit ;we're working!
;
endfor
return,lflx
end