pro did2emis2em,idstr,edens,fluxes=fluxes,deblend=deblend,DEM=DEM,logT=logT,$ ldir=ldir,NH=NH,EMs=EMs,FXs=FXs,xtitle=xtitle,ytitle=ytitle,$ title=title,verbose=verbose,n_e=n_e, _extra=e ;+ ;procedure did2emis2em ; for a specified line, read in emissivities at various densities, ; compute fluxes, and invert to get emission measure estimates at ; each density. ; ;syntax ; did2emis2em,idstr,edens,fluxes=fluxes,/deblend,DEM=DEM,logT=logT,$ ; ldir=ldir,NH=NH,EMs=EMs,FXs=FXs,xtitle=xtit,ytitle=ytit,title=tit,$ ; verbose=v,dWVL=dWVL,eps=eps,/incieq,mapping=mapping,chifil=chifil,$ ; chidir=chidir,eqfile=eqfile,abund=abund,/noph,effar=effar,wvlar=wvlar,$ ; /kev,tol=tol,fH2=fH2,He1=He1,HeII=HeII,/Fano,/wam,/bam,/mam,/Zbeda,$ ; /Ibeda,Wbeda=Wbeda,/Lbeku,wform=wform,wstyle=wstyle,ziform=ziform,$ ; /nuthin,/noday,/notime,/nouser,/nopack,stacol=stacol,stasiz=stasiz,$ ; stathk=stathk,charsize=charsize,align=align ; ;parameters ; idstr [INPUT; required] structure containing ID information ; (see LINEID for description) ; edens [I/O] array of electron densities [cm^-3] ; * default is [1e8,1e14] ; * if not specified, insufficiently specified, or otherwise ; meaningless, uses default and overwrites input ; ;keywords ; fluxes [I/O] observed fluxes ; * if size does not match the number of components ; or the number of IDs, will instead use IDSTR.(#).FLUX ; and overwrite the input ; * see UPDATID() for details ; deblend [INPUT] if _not_ set (which is the default), squishes multiple ; IDs of a single feature into one single flux and calculates ; only a "composite" ID. on the other hand, if set, leaves ; the multiple IDs alone. ; DEM [INPUT] differential emission measure to use in computing ; the fluxes ; * default is to use 1e12 [cm^-5] at each specified LOGT ; logT [INPUT] temperatures at which DEM is defined ; * if not specified, set to 6.5, unless DEM has more elements, ; in which case is interpolated into the 4..8 range ; ldir [INPUT] string array of database directories to use to ; search for the emissivities ; * default is '$CHIANTI' ; * if size does not match either number of IDSTR components ; or IDs, only the first element is used ; NH [INPUT] H column density [cm^-2] ; EMs [OUTPUT] array of emission measures EMS[EDEN,WVL] that ; are plotted. ; * at any given point, EMS = \int DEM[LOGT]*dLOGT, ; unless DEM has one or two elements only, in which ; case EMS = SUM(DEM) ; * what this implies is that the input DEM is scaled ; at each density to match the observed flux at each ; wavelength and this scaled DEM is represented in the ; plot by the summed measure. This ensures that all ; wavelengths considered are compared over the same ; temperature range, and the ambiguity regarding the ; temperature of peak emissivity, etc. has been removed. ; FXs [OUTPUT] the predicted fluxes FXS[EDEN,WVL] which are used ; to scale the EMs above. ; xtitle [INPUT] passed to PLOT ; ytitle [INPUT] passed to PLOT ; title [INPUT] passed to PLOT ; verbose [INPUT] controls chatter ; _extra [INPUT] pass defined keyword values to subroutines: ; ID2EMIS2ID: DWVL ; RD_LIST: EPS, INCIEQ, MAPPING ; FOLD_IONEQ: CHIFIL, EQFILE ; RD_IONEQ: CHIDIR ; SQUISHEM: ABUND ; LINEFLX: ABUND, NOPH, EFFAR, WVLAR, KEV ; ARRAYEQ: TOL ; ISMTAU: fH2, He1, HeII, FANO, WAM, BAM, MAM ; IDLABEL: ZBEDA, IBEDA, WBEDA, LBEKU, WFORM, WSTYLE, ZIFORM ; STAMPLE: NUTHIN,NODAY,NOTIME,NOUSER,NOPACK,STACOL,STASIZ,STATHK ; XYOUTS: CHARSIZE,ALIGN ; n_e [IGNORE] here only to trap boxing gloves ; ;restrictions ; requires IDL5.3+ because of use of STRMATCH, etc. ; requires subroutines ; BAMABS ; CAT_LN ; FOLD_IONEQ ; GETABUND ; ID2EMIS2ID ; IDLABEL ; INICON ; ISMTAU ; LAT2ARAB ; LINEFLX ; RDABUND ; RD_IONEQ ; RD_LINE ; RD_LIST ; READ_IONEQ ; SQUISHEM ; SYMB2ZION ; SYZE ; UPDATID ; WHEE ; ZION2SYMB ; ;side-effects ; generates a plot which erases any plot that exists beforehand ; ;history ; vinay kashyap (JanMMI; based on DID2EM) ; improved color-scale setting for 24-bit consoles (VK; FebMMI) ;- ; usage nid=n_tags(idstr) if nid eq 0 then begin print,'Usage: did2emis2em,idstr,edens,fluxes=fluxes,/deblend,DEM=DEM,$' print,' logT=logT,ldir=ldir,NH=NH,EMs=EMs,FXs=FXs,xtitle=xt,ytitle=yt,$' print,' title=tt,verbose=v, dWVL=dWVL,/incieq,mapping=mapping,$' print,' chifil=chifil,chidir=chidir,eqfile=eqfile,abund=abund,/noph,$' print,' effar=effar,wvlar=wvlar,/kev,tol=tol,fH2=fH2,He1=He1,$' print,' HeII=HeII,/Fano,/wam,/bam,/mam,/Zbeda,/Ibeda,Wbeda=Wbeda,$' print,' /Lbeku,wform=wform,wstyle=wstyle,ziform=ziform,/nuthin,$' print,' /noday,/notime,/nouser,/nopack,stacol=stacol,stasiz=stasiz,$' print,' stathk=stathk,charsize=charsize,align=align' print,' generates emission measure estimates at a variety of electron' print,' density values and plots them' return endif if float(strmid(!version.RELEASE,0,3)) lt 5.3 then begin message,'Requires IDL v5.3 or higher. Returning',/info return endif ; verify input ok='ok' & tagid=tag_names(idstr) if tagid[0] ne 'WVL' then ok='ID Structure not in right format' else $ if tagid[0] eq 'WVL_COMMENT' then ok='ID Structure contains no data' else $ if tagid[0] eq 'Z' then ok='ID structure appears to have been stripped' if ok ne 'ok' then begin message,ok,/info & return endif ; ; how many components? how many IDs? which wavelengths? obswvl=idstr.WVL & nwvl=n_elements(obswvl) wvl=abs(idstr.(1).WVL) & for i=1L,nwvl-1L do wvl=[wvl,abs(idstr.(i+1).WVL)] mwvl=n_elements(wvl) ; check other inputs nden=n_elements(edens) if nden lt 2 then begin message,'Setting density array to default values',/info edens=[1e8,1e14] & nden=2L endif ; check and set keywords ; verbosity v=0 & if keyword_set(verbose) then v=long(verbose[0]) > 1 ; deblend or no? if not keyword_set(deblend) then begin ;(squish multiple IDs if v gt 3 then message,'squishing together multiple IDs into one flux',/info nflux=nwvl obsflx=fltarr(nwvl) & for i=0L,nwvl-1L do obsflx[i]=total(idstr.(i+1L).FLUX) endif else begin ;)(leave them deblended if v gt 3 then message,'multiple IDs are deblended',/info nflux=mwvl obsflx=idstr.(1).FLUX & for i=1L,nwvl-1L do obsflx=[obsflx,idstr.(i+1).FLUX] endelse ;DEBLEND) ; fluxes nf=n_elements(fluxes) if nf eq 0 then fluxes=obsflx if nf gt 0 and nf ne nflux then begin message,'Ignoring input fluxes; using IDSTR.(#).FLUX',/info fluxes=obsflx endif ;fold in new fluxes strid=updatid(idstr,fluxes,verbose=v) ;DEM and logT mD=n_elements(DEM) & mT=n_elements(logT) if mT eq 0 then begin ;(set LOGT if mD le 1 then begin if mD eq 0 then begin DEM=1d12 if v ge 5 then message,'DEM undefined: setting to default',/info endif logT=6.5 & mT=1 & mD=1 if v ge 5 then message,'LOGT undefined; defaulting to '+$ strtrim(logT,2),/info endif else begin mT=mD & logT=findgen(mT)*(4./float(mD-1.))+4. if v ge 5 then message,'LOGT undefined; interpolating to ['+$ strtrim(min(logT),2)+','+strtrim(max(logT),2)+']',/info endelse endif ;LOGT) if mD eq 0 then begin ;(set DEM DEM=dblarr(mT)+1d12 & mD=mT if v ge 5 then message,'DEM undefined: setting to default',/info endif ;DEM) if mD gt 0 then begin ;(convert log(DEM) to DEM if max(DEM) lt 100. then begin DEM=10.D^(DEM) message,'converting input DEM from log10',/info unlog=1 ;check this keyword later to convert back if needed endif endif ;log(DEM)->DEM ;database directories ldbdir=strarr(mwvl)+'$CHIANTI' ;set line database directory for each ID mdir=n_elements(ldir) & szd=size(ldir) & nszd=n_elements(szd) if mdir gt 0 and szd[nszd-2] ne 7 then begin ;(LDIR is not a string mdir=1L if fix(ldir[0]) eq 1 then ldir='$SPEX' else $ if fix(ldir[0]) eq 2 then ldir='$CHIANTI' else $ ldir='$SCAR' if v gt 1 then message,'Setting emissivity database to '+ldir[0],/info endif ;decoding LDIR) if mdir gt 0 then begin ;(set LDIR ldbdir[*]=ldir[0] if mdir eq mwvl then ldbdir=ldir if mdir eq nwvl then begin k=0L for i=0L,nwvl-1L do begin mw=n_elements(strid.(i+1).WVL) ldbdir[k:k+mw-1L]=ldir[[i]] k=k+mw endfor endif endif ;LDIR) ; NH etau=fltarr(nflux)+1. if keyword_set(NH) then begin ;(set absorption if v ge 5 then message,'Computing ISM absorption',/info if keyword_set(deblend) then tau=ismtau(abs([wvl]),NH=NH, _extra=e) else $ tau=ismtau(abs([obswvl]),NH=NH, _extra=e) oo=where(tau lt 69.,moo) if moo gt 0 then etau[oo]=exp(-tau[oo]) endif ;ETAU) ; initializations atom=1 & rom=1 ;need this for IDLv<5 inicon,atom=atom,roman=rom atom=['X',atom] & rom=['',rom] ;atomic symbols EMs=dblarr(nden,nflux) ;output FXs=fltarr(nden,nflux) ;output if keyword_set(deblend) then $ labels=idlabel(idstr,idlabel_idx,/Wbeda,Lbeku=0, _extra=e) else $ labels=idlabel(squishem(idstr,_extra=e),idlabel_idx,/Wbeda,Lbeku=0,_extra=e) ; figure out the baseline, default EM if mD gt 2 then begin mdT=alog(10.)*median(logT[1:*]-logT) tmp=DEM & tmp[[0L,mD-1L]]=0.5*tmp[[0L,mD-1L]] defEM=total(tmp*mdT) endif else defEM=total(DEM) ; the basic strategy is to use ID2EMIS2ID to read in ; new emissivities at each density for each of the lines ; in the input ID structure, and to compute fluxes in ; each case. for i=0L,nden-1L do begin ;{for each density n_e=edens[i] newid=id2emis2id(strid,ldbdir,verbose=v,n_e=n_e, _extra=e) ; predflx=fltarr(nflux) & k=0L for j=0L,nwvl-1L do begin ;{for all components tmpID=newid.(j+1L) & mw=n_elements(tmpID.wvl) tmpfx=lineflx(tmpID.EMIS,tmpID.LOGT,tmpID.WVL,tmpID.Z, _extra=e) if keyword_set(deblend) then begin ;(deblend predflx[k:k+mw-1L]=tmpfx & k=k+mw endif else predflx[j]=total(tmpfx) ;squished) endfor ;J=0,NFLUX-1} predflx=predflx*etau ;ISM correction ; oo=where(predflx gt 0,moo) if moo gt 0 then EMs[i,oo]=defEM*fluxes[oo]/predflx[oo] FXs[i,*]=predflx endfor ;I=0,NDEN-1} ; plot ; first figure out plot ranges oo=where(EMs gt 0,moo) if moo eq 0 then begin message,'No lines found',/info & stop & return endif emmin=min(EMs[oo],max=emmax) if emmax/emmin gt 10 then begin ylog=1 ;plot Y-axis in log scale emmin=emmin/5. & emmax=emmax*5. ;stretch the window drop=10.^((alog10(emmax)-alog10(emmin))/nflux) ;label drops endif else begin ylog=0 ;plot Y-axis in normal scale emmin=emmin/2. > 0 emmax=emmax*2. ;stretch the window drop=(emmax-emmin)/nflux ;label drops endelse ; denmin=min(edens,max=denmax) if denmax/denmin gt 100 then begin xlog=1 ;plot X-axis in log scale denmin=denmin/10. & denmax=denmax*10. ;stretch the window endif else begin xlog=0 ;plot X-axis in normal scale denmin=denmin-10. > 0. denmax=denmax+10. ;stretch the window endelse ; ; set up the plot plot,[1],/nodata,xrange=[denmin,denmax],yrange=[emmin,emmax],$ xlog=xlog,ylog=ylog,xstyle=1,ystyle=1,$ xtitle=xtitle,ytitle=ytitle,title=title stample, _extra=e ; ; and plot dncolors=256. > !D.N_COLORS ;24-bit color screen temporary fix curvcol=fix(200*(!d.n_colors)/dncolors+1) lablcol=fix(150*(!d.n_colors)/dncolors+1) lablpos=fltarr(nflux) for i=0L,nflux-1L do begin y=reform(EMs[*,i]) & oo=where(y gt 0,moo) if moo gt 0 then begin oplot,edens[oo],y[oo],color=curvcol lpos=max(y,ipos) & ok=where(lpos eq lablpos,mok) if mok eq 0 then lablpos[i]=lpos else lablpos[i]=lablpos[ok[mok-1L]]-drop if ipos eq 0 then align=1 else align=0 xyouts,edens[ipos],lpos,labels[i],align=align,col=lablcol, _extra=e if mok ne 0 then oplot,[max(edens),denmax],[max(y),lpos],col=lablcol endif endfor if keyword_set(unlog) then DEM=alog10(DEM) return end