function lineflx,line,logT,wvl,Z,DEM=DEM,temp=temp,abund=abund,noph=noph,$ nhne=nhne,effar=effar,wvlar=wvlar,ikev=ikev,noabund=noabund,$ regrid=regrid, _extra=e ;+ ;function lineflx ; computes observeable counts [ct/s] or flux [ergs/s] from spectral ; line for specified DEM by simple trapezoidal integration over ; temperatures. ; if effective area is not specified, output will be [ph/s/cm^2] ; or [ergs/s/cm^2] ; ; input line emissivities are in [1e-23 ergs cm^3/s]. these include ; ionization balance, but not (necessarily) abundances. the line ; emissivities are multiplied by the differential emission measure ; and element abundance to get fluxes, then divided by the line energy ; to get photon fluxes, and multiplied by the supplied effective areas ; to get count rates. no effort is made to determine line profiles ; or to include spectral responses (beyond the effective areas). ; ;syntax ; fx=lineflx(line,logT,wvl,Z,DEM=DEM,/temp,abund=abund,/noph,$ ; nhne=nhne,effar=effar,wvlar=wvlar,/ikev,/noabund) ; ;parameters ; line [INPUT; required] array of line cooling emissivities ; in units of 1e-23 erg cm^3/s ; * if 2D array, LINE==LINE(logT,WVL) ; ; * WARNING: will return garbage if given 1D array LINE(WVL) ; (or even 2D array LINE(1,WVL)) ; use a for-loop to handle such a case ; ; * WARNING: will get converted to 1-element vector if input ; as scalar ; logT [INPUT; required] array of log10(Temperature [K]) at ; which emissivities are given. ; * array size MUST match that of LINE ; * if not regularly gridded, remember to set REGRID! ; wvl [INPUT; required] wavelength(s) [Angstrom] of lines at ; which LINE is given ; * array size MUST match LINE ; Z [INPUT] atomic number of element that generates each WVL ; * if Z has less elements than WVL, Z is ignored ; * if Z is to be ignored, all lines are assumed to be from ; same element, and abundance values are ignored ; ;keywords ; DEM [INPUT] Differential Emission Measure at each T [cm^-5/logK] ; * default is a constant=1e14!! ; * if array size does not match that of LOGT, then DEM is ; assumed to cover the same range and is linear interpolated ; * if defined at only 1 temperature, assumed to be EM [cm^-5] ; temp [INPUT] if set, assumes that logT is actually in T [K], and ; that DEM is given in units of [cm^-5/K] ; * ignored if logT is defined at only one bin ; abund [INPUT] abundances relative to H (abund(0)=1) ; * abund(Z-1) contains the abundance for element Z ; * if array size is smaller than the largest present Z, ; the last element is replicated to fill the gap ; * default: Anders & Grevesse (1989) ; noph [INPUT] if set, does not do the conversion to [ph] from [ergs] ; nhne [INPUT] the ratio of N(H)/n_e needed to complete the ; intensity calculation ; * if not set, assumes 0.83, which is the limiting value ; for high-temperature plasma with approximately cosmic ; abundances ; * if array and size does not match LOGT, then uses only ; the first element ; * if Z is not defined, or is 1, then it is **IGNORED** ; * see discussion in the documentation of POPSOL() ; effar [INPUT] effective area [cm^2] ; * if not set, assumed to be 1 cm^2 ; wvlar [INPUT] wavelengths at which effective area is defined. ; * if not set, assumed to be WVL ; * array size MUST match that of EFFAR ; ikev [INPUT] if set, assumes that WVL and WVLAR are in [keV], ; NOT [Angstrom] ; noabund [INPUT] if set, ABUND values are ignored in the calculations, ; same as though they (or Z) had been set to 1 ; * but NHNE gets applied ; regrid [INPUT] if set, assumes that LOGT is irregularly sampled ; and regrids it (and LINE) to a regular grid ; * NOT IMPLEMENTED YET ; _extra [INPUT] junk -- here only to prevent crashing the program ; ;subroutines ; GETABUND ; WHEE ; ;usage examples ; * f=lineflx(rd_line(logP=20,wvl=w,logT=T,/allah),T,w,dem=10.D^(5*t)) ; * fl=rd_line(wr=16.776,wvl=w,logT=T,/allah) & t=10.^(t) ; f=lineflx(fl,T,w,dem=10.D^(5*alog10(t)),/temp) ; ;history ; vinay kashyap (Nov96) ; added keywords EFFAR and WVLAR (VK; Jan96) ; added call to WHEE, removed call to KILROY, ; ensured DEM stretches to fit logT, added case of EM masquerading ; as DEM, removed the [/sr] from the output (VK; Feb97) ; changed integration style; restructured to ignore effective area ; if not supplied; added keyword REGRID (VK; Apr97) ; corrected ln v/s log10 bug (VK; Jun97) ; added keyword NOPH; bug correction -- dlogT should in in base 10 ; (VK; Jan98) ; changed keyword KEV to IKEV (VK; DecMM) ; EFFAR and WVLAR ignored if set to 0 (VK; JanMMI) ; added keyword NHNE (VK; Jun02) ; added keyword NOABUND (VK; Apr03) ; changed bad input check from LINE(0)=-1 to total(LINE)=-N(LINE) ; (VK; Feb06) ; changed back integration style to simple from trapezoidal (VK; Mar06) ; changed check for whether WVL are given (VK; Jun07) ; bug fix: was including 0.83 factor even when continuua and broad-band ; emissivities were being used [thx P.Testa]; now anything with Z=1 ; will exclude the NHNE correction -- for H lines, fold the NHNE factors ; into the emissivities manually (VK; Dec10) ;- ; check input ;is the cooling function supplied? if n_elements(line) eq 0 then begin print,'Usage: fx=lineflx(line,logT,wvl,Z,DEM=DEM,/temp,abund=abund,$' print,' /noph,nhne=nhne,effar=effar,wvlar=wvlar,/ikev,/noabund)' print,' compute observed flux ([ph/s/cm^2] or [ph/s]) from spectral line' return,0. endif szc=size(line) & nt=szc(1) & nw=szc(2) if szc(0) eq 0 and szc(1) ne 0 then begin ;LINE is input as a scalar -- convert to vector line=[line] & szc=size(line) & nt=szc(1) & nw=1L endif if szc(0) eq 1 then nw=1L ;was the cooling function properly read? ;if line(0) eq -1. then begin -- this is needlessly restrictive if total(line) eq -1L*n_elements(line) then begin ;(means filled with -1s message,' No lines input; returning',/info & return,line endif ;LINE==-1) ;are the temperatures given? if n_elements(logT) ne nt then begin ;size mismatch message,'Cooling Function and Temperature arrays do not match',/info return,line endif tt=[logT] ;are the wavelengths given? ang=[0.] & if keyword_set(wvl) then ang=[wvl] if szc(0) ne 1 then begin if n_elements(wvl) ne nw then begin ;size mismatch message,'Cooling Function and wavelength arrays do not match',/info return,line endif ang=[wvl] endif ang=abs(ang) ;the CHIANTI case ;if wavelengths are given, so also should atomic numbers be atom=intarr(nw)+1 ;default is H ;if ang(0) ne 0 then begin ;wavelengths are given if total(ang) ne 0 then begin ;wavelengths are given nz=n_elements(z) if nz gt 0 and nz le nw then atom(0)=z if nz gt nw then atom(*)=z(0:nw-1) endif oz=where(atom eq 1,moz) ;is the DEM given? diffem=dblarr(nt)+1d14 ;default if keyword_set(dem) then begin ndem=n_elements(dem) if ndem eq 1 then diffem(*)=dem ;scalar if ndem eq nt then diffem=[dem] ;what the doc ordered if ndem gt 1 and ndem ne nt then begin ;interpolate message,'Interpolating DEM onto a default grid',/info t0=logT(0) & t1=logT(nt-1) & ti=findgen(ndem)*(t1-t0)/(ndem-1.)+t0 tmp=alog10(dem>1e-20) diffem=interpol(tmp,ti,logT) & diffem=10.^(diffem) oo=where(diffem lt 1e-19) & if oo(0) ne -1 then diffem(oo)=0. endif endif ;are the temperatures in T or logT? if keyword_set(temp) then begin if nt gt 1 then diffem=diffem*tt ;[cm^-5/K]->[cm^-5/logK] ;(if nt eq 1, then DEM is actually EM) tt=alog10(abs(logT)) ;convert to log10(T) endif ;is the N(H)/N(e) fraction given? nh_ne=fltarr(nt)+0.83 & nhe=n_elements(nhne) if nhe gt 0 then begin if nhe eq 1 or nhe ne nt then begin if nhne(0) le 0 then nh_ne(*)=1. else nh_ne(*)=nhne(0) endif else nh_ne=float(nhne)>0. endif if nz eq 0 then nh_ne(*)=1. ;ignore if Z is not given if moz ne 0 then nh_ne(oz)=1. ;or is "H" because this is how ;continuum and solar calcs are done ;check if temperature gridding is uniform if nt gt 1 then begin dlogt=tt(1:*)-tt & mdt=total(dlogt)/(nt-1.) sdt=sqrt(total((dlogt-mdt)^2)/(nt-1.)) if sdt/mdt gt 0.1 and not keyword_set(regrid) then $ message,'Temperature gridding is not regular, but ignoring!',/info endif ;regrid if required if keyword_set(regrid) then begin message,' sorry, REGRID not implemented yet!',/info endif ;convert dlog10 to dl ;if nt gt 1 then mdt=alog(10.^(tt(1)-tt(0))) ;THIS IS WRONG!! if nt gt 1 then mdt=alog(10.)*mdt ;are the abundances given? if n_elements(abund) ne 30 then abund=getabund('anders & grevesse') ;are the effective areas and wavelengths given? nar=n_elements(effar) & nwv=n_elements(wvlar) if nar eq 1 and not keyword_set(effar) then nar=0 if nwv eq 1 and not keyword_set(wvlar) then nwv=0 if nar eq 0 and nwv eq 0 then area=0*ang+1. else begin if nar ne 0 then areff=[effar] else areff=[1.,1.] if nwv ne 0 then arwvl=[wvlar] else arwvl=[min(abs(ang)),max(abs(ang))] if nar ne nwv then begin message,'Effective area v/s wavelength mismatch; doing nothing',/info area=0*ang+1. endif else begin ow=sort(arwvl) & arwvl=arwvl(ow) & areff=areff(ow) ;sort area=interpol(areff,arwvl,abs(ang))>0 ;effective areas at WVLs endelse endelse ; integrate flx=dblarr(nw) & avdem=diffem*1d-23 for iw=0L,nw-1L do begin if keyword_set(noabund) then zab=1. else $ zab=(abund([atom(iw)-1]))(0) ;abundance if nt eq 1 then begin ;(EM, not DEM cfnc=[nh_ne(0)*line(0)] flx(iw)=zab*avdem*cfnc ;[ergs/s/cm^2]=[1e23 cm^-5]*[1e-23 ergs cm^3/s] endif else begin ;)(DEM, not EM cfnc=reform(nh_ne(*)*line(*,iw)) tmp=zab*avdem*cfnc ;tmp([0,nt-1])=0.5*tmp([0,nt-1]) ;; the 0.5 weighting at the ends is due to trapezoidal integration ;which has been commented out by VK because it is just very annoying ;to have those 2x dips at the temperature grid extremes (Mar2006) flx(iw)=mdt*total(tmp,/nan) ;[ergs/s/cm^2]=[1e23 cm^-5/logK]*[logK]*[1e-23 ergs cm^3/s] ;NOTE: divide by (4.D*!dpi) to include [/sr] -- not done here ;;;;flx(iw)=int_tabulated(tt,zab*avdem*cfnc) ;;;; uses 5-pt Newton-Cotes, but slows LINEFLX to a ;;;; crawl even on a SPARC Ultra endelse ;NT.ne.1) endfor ; convert from ergs/s/cm^2 to ph/s if not keyword_set(noph) then begin h=6.626176e-27 & c=2.9979e10 ;[ergs s] & [cm s^-1] nrg=h*c*1e8/abs(ang) ;[ergs] if keyword_set(ikev) then nrg=1.6021892e-9*abs(ang) ;[ergs] flx=flx/nrg ;[ph/s/cm^2] endif flx=flx*area ;[(erg|ph)/s] return,flx end