function lsd,wrange,fluxes,wvls,elem=elem,edens=edens,DEM=DEM,tlog=tlog,$ dbdir=dbdir,ceiling=ceiling,flor=flor,ratmax=ratmax,flxmax=flxmax,$ outZ=outZ,outIon=outIon,outIDX=outIDX, _extra=e ;+ ;function lsd ; returns list of density sensitive lines in given wavelength range ; ;syntax ; ld=lsd(wrange,fluxes,wvls,elem=elem,edens=edens,DEM=DEM,tlog=tlog,$ ; dbdir=dbdir,ceiling=ceiling,flor=flor,ratmax=ratmax,flxmax=flxmax,$ ; outZ=outZ,outIon=outIon,outIDX=outIDX, tol=tol,chifil=chifil,$ ; chidir=chidir,eqfile=eqfile,/temp,abund=abund,/noph,effar=effar,$ ; wvlar=wvlar,/ikev) ; ;parameters ; wrange [INPUT; required] wavelength range in which to search for ; density sensitive lines ; * if only one element is given, set the range to be 10% of ; it on either side ; * if the single element is a string, apply RD_LIST rules ; to decipher it ("WVL +- dW", "WMIN,WMAX", "WMIN-WMAX" ; are all allowed) ; fluxes [OUTPUT] 2D array of predicted fluxes for lines that are ; found to be density sensitive, FLUXES(EDENS,WVLS) ; wvls [OUTPUT] array of wavelengths of lines that are found to ; be density sensitive ; ;keywords ; elem [INPUT] if set, confine attention to just those elements and ; ionic species specified. ; * default is to look at all available lines ; edens [I/O] array of electron densities [cm^-3] at which to ; check for density sensitivity. ; * default is [1e8,1e14] ; * if not specified, insufficiently specified, or otherwise ; meaningless, uses default ; DEM [INPUT] input array of differential emission measures [cm^-5] ; * default is to use 1e12 [cm^-5] at specified TLOG ; tlog [INPUT] log(temperatures) at which DEM is defined ; * if not specified, set to 6.0, unless DEM has more steps, ; in which case interpolated to the 4..8 range ; dbdir [INPUT] line database directory to use to search for the ; density sensitive lines ; * default is $CHIANTI ; ceiling [INPUT] maximum relative variation the fluxes of a given ; line must have before it is accepted as density sensitive ; * default is 2 ; * if 0 < CEILING < 1, it is assumed that the *reciprocal* ; of CEILING is given (e.g., CEILING=0.5 is translated to 2) ; * if CEILING < 0, assumed to be percentages (e.g., CEILING=-1 ; is translated to 1/0.01 == x100) ; flor [INPUT] if set, ignores any line which has a peak flux ; a factor FLOR below the flux of the strongest of the ; lines filtered through CEILING. ; * default is 0 ; * if FLOR > 1, assumed that the reciprocal of FLOR is given, ; i.e., FLOR=100 translates to 0.01 ; * if FLOR < 0, taken to be the absolute flux threshold, in ; whatever units it is being calculated in ; * NOTE on difference between CEILING and FLOR:- ; CEILING works on ratios, and FLOR works on the intensities ; ratmax [OUTPUT] the maximum deviations found for each of the ; density-sensitive lines found ; flxmax [OUTPUT] the peak fluxes for each of the density sensitive ; line found ; outZ [OUTPUT] atomic numbers ; outIon [OUTPUT] atomic numbers ; outIDX [OUTPUT] indices which passed the test ; * be careful in how you use this. a common mistake is to ; assume that it matches with whatever line emissivity ; database that has been read in. that might very well be ; true, but is definitely not guaranteed. ; ; _extra [INPUT] use this to pass defined keywords to subroutines:- ; FOLD_IONEQ: CHIFIL ; RD_IONEQ: CHIDIR, EQFILE ; LINEFLX: TEMP, ABUND, NOPH, EFFAR, WVLAR, IKEV ; ARRAYEQ: TOL ; ;example usage ; ls=lsd('10-120',dbdir='$CHIANTI',edens=[1e9,1e12],ceiling=2,ratmax=rmx) ; for i=0,n_elements(ls)-1 do print,ls(i),1./rmx(i),i ; ; ; eden=10.^(findgen(6)+8) ; ls=lsd('113.35?0.1',flx,elem='Fe20',dbdir='$CHIANTI',edens=eden) ; plot,eden,flx ; ;subroutines ; MK_DEM ; RD_LINE ; FOLD_IONEQ ; RD_IONEQ ; READ_IONEQ (CHIANTI subroutine) ; LINEFLX ; ARRAYEQ ; KABOOM ; INICON ; IS_KEYWORD_SET ; ;history ; vinay kashyap (Nov98) ; changed default EDENS to [1e8,1e14]; added keywords DEM,TLOG,OUTZ, ; OUTION; added call to MK_DEM; changed keyword FLOOR to FLOR; changed ; behavior of FLOR<0 (VK; Dec98) ; modified ion balance calcs (VK; 99May) ; changed call to INITSTUFF to INICON (VK; 99May) ; added keyword outIDX (VK; JanMMI) ; updated for IDL5.6 keyword_set([0]) behavior change for vectors ; (VK; 20Mar2006) ;- ; usage szw=size(wrange) & nszw=n_elements(szw) & nw=szw(nszw-1L) if szw(nszw-2) eq 0 or nw gt 2 then begin print,'Usage: ld=lsd(wrange,fluxes,wvls,elem=elem,edens=edens,DEM=DEM,$' print,' tlog=tlog,dbdir=dbdir,ceiling=ceiling,flor=flor,ratmax=rmx,$' print,' flxmax=fmx,outZ=outZ,outIon=outIon,outIDX=outIDX, tol=tol,$' print,' chifil=chifil,chidir=chidir,eqfile=eqfile,/temp,abund=abund,$' print,' /noph,effar=effar,wvlar=wvlar,/ikev)' print,' return list of density sensitive lines' return,'' endif ; check input ok='ok' if nw gt 2 then ok='cannot understand WRANGE -- too many elements?' if szw(nszw-2) eq 7 then begin ;(string input if nw gt 1 then message,'using only first element of WRANGE',/info c=strtrim(wrange(0),2) j=strpos(c,'[',0) & if j ge 0 then c=strmid(c,j+1,strlen(c)-j-1) j=strpos(c,']',0) & if j ge 0 then c=strmid(c,0,j) j=strpos(c,'(',0) & if j ge 0 then c=strmid(c,j+1,strlen(c)-j-1) j=strpos(c,')',0) & if j ge 0 then c=strmid(c,0,j) j=strpos(c,'{',0) & if j ge 0 then c=strmid(c,j+1,strlen(c)-j-1) j=strpos(c,'}',0) & if j ge 0 then c=strmid(c,0,j) j=0 ; W (exact match) if strpos(c,'+-',0) ge 0 then j=1 ; W +- dW if strpos(c,'-',0) ge 0 and j ne 1 then j=2 ; Wmin - Wmax if strpos(c,',',0) ge 0 then j=3 ; Wmin , Wmax if strpos(c,'?',0) ge 0 then j=4 ; W ? dW (same as W +- dW) case j of ;{for each case 1: begin cc=str_sep(c,'+-') & ncc=n_elements(cc) w=float(cc(0)) & if ncc eq 1 then dw=0.1*w else dw=float(cc(1)) w1=w-dw & w2=w+dw end 2: begin cc=str_sep(strtrim(c,2),'-') & ncc=n_elements(cc) w1=float(cc(0)) & if ncc eq 1 then w2=w1+1. else w2=float(cc(1)) end 3: begin cc=str_sep(strtrim(c,2),',') & ncc=n_elements(cc) w1=float(cc(0)) & if ncc eq 1 then w2=w1+1. else w2=float(cc(1)) end 4: begin cc=str_sep(strtrim(c,2),'?') & ncc=n_elements(cc) w=float(cc(0)) & if ncc eq 1 then dw=0.1*w else dw=float(cc(1)) w1=w-dw & w2=w+dw end else: begin w=float(c) & w1=0.9*w & w2=1.1*w endelse endcase ;J} endif ;WRANGE is a string) if szw(nszw-2) le 5 and szw(nszw-2) ge 2 then begin ;(numeric input if nw eq 1 then begin w=0.+wrange(0) & w1=0.9*w & w2=1.1*w endif else begin w1=0.+wrange(0) & w2=0.+wrange(1) endelse endif ;WRANGE=[Wmin,Wmax]) if w1 gt w2 then begin w=w1 & w1=w2 & w2=w ;switch if necessary endif if szw(nszw-2) gt 7 or szw(nszw-2) lt 2 then ok='WRANGE: Unknown data type' if ok ne 'ok' then begin message,ok,/info & return,'' endif ; check keywords if not keyword_set(elem) then elem=0 ; nden=n_elements(edens) if nden lt 2 then begin edens=[1e8,1e14] & nden=2L endif ; if not keyword_set(dbdir) then dbdir='$CHIANTI' ; cthr=0.5 if keyword_set(ceiling) then begin if ceiling(0) gt 1 then cthr=1./ceiling(0) if ceiling(0) gt 0 and ceiling(0) le 1 then cthr=ceiling(0) if ceiling(0) lt 0 then cthr=abs(ceiling(0))/100. endif if cthr gt 1 then cthr=1./cthr ; fthr=1e-6 if keyword_set(flor) then begin if flor(0) gt 0 then fthr=flor(0) if flor(0) lt 0 then fthr=1e-6 endif ;if fthr gt 1 then fthr=1./fthr ; need some initialization atom=1 & rom=1 & inicon,atom=atom,roman=rom nlin=0L ; extract the emissivities, fold in ion balance, get line fluxes for i=0L,nden-1L do begin ;{for each density ff=rd_line(elem,n_e=edens(i),wrange=[w1,w2],dbdir=dbdir,$ wvl=wvl,logT=logT,Z=Z,ion=ion,jon=jon, _extra=e) if ff(0) le -1 then goto,skip ;{no lines read.. wvl=abs(wvl) ff=fold_ioneq(ff,Z,jon,logT=logT, _extra=e) if keyword_set(DEM) then D_E_M=DEM if n_elements(D_E_M) ne n_elements(logT) then begin message,'interpolating DEM to fit the temperature grid',/info D_E_M=mk_dem('interpolate',logT=logT,indem=DEM,pardem=tlog) endif fx=lineflx(ff,logT,wvl,Z,DEM=D_E_M, _extra=e) if not keyword_set(ii) then begin ;(check consistency at EDENS(i) ii=1 ;first pass nlin=n_elements(fx) owvl=wvl & outZ=Z & oution=ion ;save them.. fmax=fx endif else begin ; subsequent passes nln=n_elements(fx) & aeq=arrayeq(wvl,owvl,_extra=e) if nln ne nlin or aeq ne 1 then begin ;(problem.. message,'holes exist in the line database',/info kaboom,/flash,/beep ffx=0*fx & k=0L if nln gt nlin then begin ;(more lines than before.. for j=0L,nlin-1L do begin if abs(owvl(j)-wvl(k)) lt 1e-3 and outZ(j) eq Z(k) and oution(j) eq ion(k) then begin ffx(j)=fx(k) & k=k+1L endif endfor endif else begin ;)(fewer lines than before.. for j=0L,nln-1L do begin if abs(owvl(k)-wvl(j)) lt 1e-3 and outZ(k) eq Z(j) and oution(k) eq ion(j) then begin ffx(k)=fx(j) & k=k+1L endif endfor endelse ;NLN v/s NLIN) fx=ffx endif ;..problem) fmax = fmax > fx ;locate the maxima endelse ;EDENS consistency check) ; form flux array.. if not is_keyword_set(flux) then flux=fx else flux=[flux,fx] if not is_keyword_set(DenE) then DenE=edens(i) else DenE=[DenE,edens(i)] skip: ;get here directly if no lines were read..} endfor ;I=0,NDEN-1} ; simple error checks if nlin eq 0 then return,'' ;nothing was found nfx=n_elements(flux) & if nfx lt nlin then return,'' ;nothing was found nd=n_elements(DenE) & if nd lt 2 then return,'' ;nothing was found ok=where(fmax gt 0,mok) & if mok eq 0 then return,'' ;nothing varied nd2=nfx/nlin & if nd ne nd2 then message,'bug' ; now find the density sensitive lines frat=fltarr(nlin)+1. for i=0L,nd-1L do begin ;{check all lines at each density k1=i*nlin & k2=k1+nlin-1L frat1=flux(k1:k2) & frat(ok)=frat1(ok)/fmax(ok) < (frat(ok)) endfor ;I=0,ND-1} ; ot=where(frat le cthr,mot) if mot eq 0 then begin message,'no density sensitive lines were found',/info return,'' endif if keyword_set(flor) then begin fmaxmax=max(fmax(ot)) if fthr lt 1 then ot=where(frat lt cthr and fmax ge fthr*fmaxmax,mot) else $ ot=where(frat lt cthr and fmax ge abs(flor(0)),mot) if mot eq 0 then message,'no lines were found for this flux threshold',/info endif ; the output sep=' ' ls=strarr(mot) for i=0L,mot-1L do ls(i)=atom(outZ(ot(i))-1)+' '+rom(oution(ot(i))-1)+sep+$ strtrim(string(owvl(ot(i)),'(g12.5)'),2)+sep+dbdir ratmax=frat(ot) flxmax=fmax(ot) outZ=outZ(ot) outIon=outIon(ot) outIDX=ot np=n_params() if np gt 1 then begin fluxes=fltarr(nd,mot) for i=0L,nd-1L do begin k1=i*nlin & k2=k1+nlin-1L fx=flux(k1:k2) & fluxes(i,*)=fx(ot) endfor endif if np gt 2 then wvls=owvl(ot) return,ls end