function flux_to_em,emis,flux=flux,logT=logT,wvl=wvl,Z=Z,NH=NH,defEM=defEM,$ noph=noph,thresh=thresh, _extra=e ;+ ;function flux_to_em ; return emission measures (-1 where undeterminable) at various ; temperatures consistent with given line fluxes and line emissivities. ; ; for each temperature, assume a delta-function emission measure, ; compute line fluxes seen through some instrument, account for ; interstellar absorption, and scale to match the observed flux. ; ;syntax ; em=flux_to_em(emis,flux=flux,logT=logT,wvl=wvl,Z=Z,NH=NH,$ ; defEM=defEM,noph=noph,thresh=thresh, abund=abund,/temp,/ikev,$ ; effar=effar,wvlar=wvlar,fh2=fh2,he1=he1,heII=heII,/fano) ; ;parameters ; emis [INPUT; required] 2D array of line emissivities, EMIS(LOGT,WVL) ; * if 1-D, assumed to be EMIS(LOGT) ; ;keywords ; flux [INPUT] observed fluxes ; * if not given, assumed to be 1 [(ph|erg)/s/cm^2] ; * size is forced to match 2nd dimension of EMIS ; if <, filled with 1s ; if >, extra elements are ignored ; * NOTE: the author realizes that there may be multiple IDs ; for the single observed line, and directs the inquirer ; to the function ID_TO_EMIS, which returns the appropriately ; concatenated emissivity table. ; logT [INPUT] log_10(Temperature [K]) at which EMIS is defined ; * if not given, assumed to go from 4 to 8 in even steps ; * interpolated if size does not match 1st dimension of EMIS ; wvl [INPUT] wavelengths at which EMIS is defined ; * if not given or if size does not match 2nd dimension of EMIS, ; (a) interstellar absorption correction is not made ; (b) all fluxes >gotta be< [ergs/s/cm^2] ; Z [INPUT] atomic numbers ; * if not given, assumed to be 1s (i.e., H. i.e., no ; abundance corrections are made) ; * size is forced to match 2nd dimension of EMIS ; if <, filled with 1s ; if >, extra elements are ignored ; NH [INPUT] H column density [cm^-2] ; * ignored if WVL is incorrect ; defEM [INPUT; default=1e14] default EM to use prior to scaling [cm^-5] ; noph [INPUT] if set, FLUX is assumed to be in [ergs/s/cm^2], and ; no attempts are made to convert things to [ph/...] ; * forcibly set if WVL is incorrect ; thresh [INPUT; default=1e-2] dynamic range in output curves, ; relative to maximum in each curve ; _extra [INPUT ONLY] use this to pass defined keywords to subroutines ; LINEFLX [ABUND,TEMP,IKEV,EFFAR,WVLAR] ; ISMTAU [FH2,HE1,HEII,FANO] ; ;restrictions ; * EFFAR and WVLAR must be correctly defined, or else the output ; will be garbage ; * requires subroutines ; LINEFLX ; WHEE ; ISMTAU ; ;history ; vinay kashyap (Oct98) ; bug correction re 0 EMIS (VK; Nov98) ; bug correction re 0 NH (VK; AugMM) ;- ; usage sze=size(emis) & nsze=n_elements(sze) & n=sze(nsze-2) if n eq 0 then begin print,'Usage: em=flux_to_em(emis,flux=flux,logT=logT,wvl=wvl,Z=Z,NH=NH,$' print,' defEM=defEM,noph=noph,thresh=thresh;' print,' LINEFLX:abund,temp,ikev,effar,wvlar; ISMTAU:fh2,he1,heII,fano)' return,-1L endif ; consistency checks edim=sze(0) ;dimension of EMIS if edim eq 0 then begin & mt=1L & mw=1L & endif ;scalar if edim eq 1 then begin & mt=sze(1) & mw=1L & endif ;EMIS=EMIS(LOGT) if edim eq 2 then begin & mt=sze(1) & mw=sze(2) & endif ;EMIS=EMIS(LOGT,WVL) if edim gt 2 then begin message,'emissivity array is too complex to understand',/info return,0*emis endif nf=n_elements(flux) & nt=n_elements(logT) & nw=n_elements(wvl) & nz=n_elements(Z) ; fx=fltarr(mw)+1. ;fluxes if nf gt 0 and nf lt mw then fx(0:nf-1)=flux if nf ge mw then fx=flux(0:mw-1) ; tlog=interpol([4.,8.],mt) ;temperatures if nt gt 0 and nt ne mt then tlog=interpol(logt,mt) if nt eq mt then tlog=logt ; if keyword_set(NH) then setabs=1 else setabs=0 ;include ISM absorption? ; ww=findgen(mw) ;wavelengths if nw eq 0 or nw ne mw then begin setnoph=1 & setabs=0 endif else ww=wvl ; zz=intarr(mw)+1 ;atomic numbers if nz gt 0 and nz lt mw then zz(0:nf-1)=Z if nz ge mw then zz=Z(0:mw-1) ; if keyword_set(noph) then setnoph=1 ;everything in [ergs/...] ; if not keyword_set(defEM) then dem=1d14 else dem=double(defEM) ;default EM [cm^-5] if not keyword_set(thresh) then thresh=1e-2 ;dynamic range ; get optical depths tau=fltarr(mw) if keyword_set(setabs) then begin tau=ismtau(ww,NH=NH,fH2=0,/Fano,_extra=e) < 69. endif ; get the predicted fluxes ct=dblarr(nt,nw) for iw=0,nw-1 do begin if edim eq 0 then line=[emis] if edim eq 1 then line=[emis(*)] if edim eq 2 then line=reform(emis(*,iw)) for it=0,nt-1 do begin tmp=lineflx(line(it),tlog(it),ww(iw),zz(iw),noph=setnoph,dem=dem, _extra=e) ct(it,iw)=tmp*exp(-tau(iw)) endfor endfor ; compare to observed and scale EM em=dblarr(nt,nw)-1 for iw=0,nw-1 do begin if edim eq 0 then begin pct=[ct] & line=[emis] endif if edim eq 1 then begin pct=[ct(*)] & line=[emis(*)] endif if edim eq 2 then begin pct=reform(ct(*,iw)) & line=reform(emis(*,iw)) endif oo=where(pct ge thresh*max(pct),moo) if max(pct) le 0 then moo=0L if moo gt 0 then begin frac=dem*fx(iw)/pct(oo) em(oo,iw)=frac(*) endif endfor return,em end