|
Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://hea-www.harvard.edu/PINTofALE/doc/EXAMPLE_VLK.html
Дата изменения: Unknown Дата индексирования: Tue Oct 2 01:41:39 2012 Кодировка: Поисковые слова: п п п п р п р п р п р п р п р п р п р п р п |
This annotated example walks a user through a PINTofALE session, sparing no details however gory, unlike the similar thread that uses .par files. During this run, the user will
# start IDL idl ; set up the PINTofALE environment ; .run initale ;NOTE: ; if INITALE fails, run the script ; @PoA_constructor ; load the data and instrument files ; restore,!ARDB+'/demo_v711Tau.save',/verbose ; copy the relevant variables to more normal looking names ; data_spec = SPC_M1P_V711TAU ;observed spectrum [ct/bin] data_wvl = LAM_M1P_V711TAU ;wavlengths for observed spectrum [Ang] data_spec_err = SPE_M1P_V711TAU ;errors on observed spectrum [ct/bin] effar = EFFAR_M1P ;effective area [cm^2] wvlar = WVLAR_M1P ;wavelengths for effective area [Ang] ; and plot them ; ; counts spectrum and errors ; plot, data_wvl,data_spec,psym=10,xtitle='Wavelength ['+!AA+']',$ ytitle='Counts/bin',title='HR 1099: HETG/MEG +1',xrange=[5.,20.],$ yr=[0.1,1200],/ys for i=0L,n_elements(data_wvl)-1L do oplot,$ data_wvl[i]*[1,1],data_spec[i]+data_spec_err[i]*[-1,1],$ color=1 & setkolor,'red',1 ; ; instrument effective area ; plot, wvlar,effar,xtitle='Wavelength ['+!AA+']',$ ytitle='Effective Area [cm!u2!n]',title='HETG/MEG +1',xr=[5.,20.] ; ; normalized spectrum ; tmp_effar=interpol(effar,wvlar,data_wvl) plot, data_wvl,data_spec/tmp_effar/EXPTIM_V711TAU,$ psym=10,xtitle='Wavelength ['+!AA+']',ytitle='Counts [cm!u-2!n s!u-1!n]',$ title='HR 1099: HETG/MEG +1',xr=[5.,20.] ; ; DEM ; plot, LOGT_V711TAU,DEM_V711TAU,psym=10,/ylog,xtitle='log!d10!n (T [K])',$ ytitle='DEM [cm!u-5!n]',title='EUVE based approximate DEM for HR1099'
; set the wavelength range of interest ; !WMIN = min(data_wvl) !WMAX = max(data_wvl) ; read in the line emissivity database ; [1e-23 ergs cm^3 s^-1] ; rd_line_linemis = rd_line(atom, n_e=!EDENS, wrange=[!WMIN,!WMAX], $ dbdir=!LDBDIR, wvl=rd_line_wvl, logT=rd_line_logT, z=rd_line_z, $ ion=rd_line_ion, jon=rd_line_jon, src=rd_line_src, fstr=rd_line_linstr, $ verbose=!VERBOSE) ; fold in line emissivities with ion-balance info ; (this is needed only if the input emissivities do not already ; contain the ion-balance info. For example, this is needed while ; using our implementation of the SPEX and CHIANTI line emissivity ; databases. on the other hand, APED, as well as the continuum ; emissivity databases, come with ion-balance already folded in.) ; rd_line_linemis = fold_ioneq(rd_line_linemis, rd_line_Z, rd_line_JON, $ logT=rd_line_LOGT, tmax=fold_ioneq_tmax, trng=fold_ioneq_trng, level=0.1, $ chifil=chifil, chidir=!CHIDIR, eqfile=!IONEQF, verbose=!VERBOSE) ;NOTE 1: ; it may be advantageous (and less of a headache later) to overwrite ; the corrected emissivities into the line emissivity structure. ; the main output, rd_line_linemis, will still have emissivites ; without ion balance folded in. ; rd_line_linstr.LINE_INT = rd_line_linemis ;NOTE 2: ; the line emissivities are read in from an existing database ; computed for a density nearest to that specified. In order ; to interpolate from emissivities computed at two densities, ; it is necessary to read in each set and to manually interpolate. ; here is one way that it can be accomplished. ; ; In practice, this should not be necessary unless synthesis of very ; density-sensitive features is required at a very specific density. ; Remember that this interpolation is only linear in log. ; n_e=5e10 & n_e1=1e10 & n_e2=1e11 lne=alog10(n_e) & lne1=alog10(n_e1) & lne2=alog10(n_e2) weight1=(lne-lne1)/(lne2-lne1) & weight2=(lne2-lne)/(lne2-lne1) ; rd_line_linemis1=rd_line(atom,n_e=n_e1,wrange=[!WMIN,!WMAX],dbdir=!LDBDIR,$ wvl=rd_line_wvl,logT=rd_line_logT,z=rd_line_z,ion=rd_line_ion,$ jon=rd_line_jon,src=rd_line_src,verbose=!VERBOSE) ; rd_line_linemis2=rd_line(atom,n_e=n_e2,wrange=[!WMIN,!WMAX],dbdir=!LDBDIR,$ wvl=rd_line_wvl,logT=rd_line_logT,z=rd_line_z,ion=rd_line_ion,$ jon=rd_line_jon,src=rd_line_src,verbose=!VERBOSE) ; rd_line_linemis=10.^(weight1*alog10(rd_line_linemis1>1e-30)+$ weight2*alog10(rd_line_linemis2>1e-30)) ; ; and do not forget to include ion balance again.. ; rd_line_linemis=fold_ioneq(rd_line_linemis,rd_line_Z,rd_line_ION,$ logT=rd_line_LOGT,tmax=fold_ioneq_tmax,trng=fold_ioneq_trng,level=0.1,$ chifil=chifil,chidir=!CHIDIR,eqfile=!IONEQF,verbose=!VERBOSE) ; ; .. and push the emissivities into the structure ; rd_line_linstr.LINE_INT = rd_line_linemis ; read in the continuum emissivity database ; [1e-23 ergs cm^3 s^-1 Ang^-1] ; (the continuum emissivities are very weakly dependent on ; changes in ion balance, and so it is already included in ; the databases) ; rd_cont_conemis = rd_cont(!CEROOT, n_e=!EDENS, wrange=[!WMIN,!WMAX], $ dbdir=!CDBDIR, abund=!ABUND, metal=1., wvl=rd_cont_cwgrid, $ logT=rd_cont_LOGT, fcstr=rd_cont_constr, verbose=!VERBOSE)
; first make sure that the DEM is defined on the same temperature
; grid as are the emissivities
;
!DEM = mk_DEM('interpolate', logT=!LOGT, indem=DEM_V711TAU, $
pardem=LOGT_V711TAU)
;NOTE:
; we require that the DEM always be defined in the form
; n^2*dV/dlnT
; modifications such as n^2 dh/dlnT, or division by (4 pi d^2),
; or multiplication by (R^2/d^2) are OK too -- as long as
; the user keeps track of things. For the sake of definiteness,
; in what follows, we will assume that the units are [cm^-5]
; compute the line fluxes
; emissivities are [ergs cm^3/s], DEM is [/cm^5], effective
; area is [cm^2]. the [ergs] are converted internally to [ph],
; so the output of LINEFLX has the units [ct/s]
;
lineflx_flx = lineflx(rd_line_linemis, rd_line_LOGT, rd_line_WVL, rd_line_Z, $
DEM=!DEM, abund=!ABUND, effar=effar, wvlar=wvlar) * EXPTIM_V711TAU
; units at this stage: [ct]
; correct for ISM absorption
;
lineflx_flx = lineflx_flx*$
exp(-ismtau(abs(rd_line_WVL), NH=!NH, fH2=!FH2, He1=!HE1, HeII=!HEII, $
/bam, abund=!ABUND, verbose=!VERBOSE))
; define a wavelength grid to hold the spectrum
;
sp_NBIN = 20000L ;number of bins
sp_DW = (!WMAX-!WMIN)/sp_NBIN ;bin width
sp_wgrid = findgen(sp_NBIN+1L)*sp_DW+!WMIN ;bin boundaries
; bin the line fluxes into a spectrum
;
linspec = hastogram(abs(rd_line_WVL), sp_wgrid, wts=lineflx_flx)
; compute continuum fluxes
; emissivities are [ergs cm^3/s/A], DEM is [/cm^5], effective
; area is [cm^2]. the [ergs] are converted internally to [ph],
; so the output of LINEFLX has the units [ct/s/A], and is
; multiplied by the bin widths [A] and the exposure time [s].
;
rd_cont_CDW = rd_cont_cwgrid[1:*]-rd_cont_cwgrid
;
contflx_flx = lineflx(rd_cont_conemis, rd_cont_LOGT, rd_cont_constr.MIDWVL, $
DEM=!DEM, effar=effar, wvlar=wvlar) * rd_cont_CDW * EXPTIM_V711TAU
; units at this stage: [ct]
;NOTE:
; LINEFLX, despite its name, is a fairly general routine, which
; integrates the emissivity function over the temperature grid
; weighted by the DEM at each wavelength. The fact that we might
; be dealing with continuum emissivities rather than line emissivities
; is immaterial to the program. Also, when the atomic numbers are
; omitted, all the wavelengths are assumed to belong to H, with
; an abundance of 1.0, and hence the output will be perfectly fine.
; correct for ISM absorption
;
contflx_flx = contflx_flx*$
exp(-ismtau(rd_cont_constr.MIDWVL, NH=!NH, fH2=!FH2, He1=!HE1, $
HeII=!HEII, /bam, abund=!ABUND, verbose=!VERBOSE))
; rebin
;
conspec = rebinw(contflx_flx, rd_cont_cwgrid, sp_wgrid, /perbin)
; add the line and continuum spectra
;
spec = linspec + conspec
; plot the result
;
plot, sp_wgrid,spec,psym=10,xr=[5.,30.],/xs,yr=[1e-2,5e3],/ys,/ylog,$
xtitle='Wavelength ['+!AA+']',ytitle='Counts [ct]',$
title='Theoretical spectrum'
oplot, sp_wgrid,linspec,color=1,psym=10
oplot, sp_wgrid,conspec,color=2,psym=10
setkolor, 'red',1 & setkolor, 'yellow',2
; convolve with lorentzian of half-width=5*bin-width to simulate LRF
; This is an example of a "general utility" PINTofALE routine that
; is not coronal plasma specific. Use a "Lorentzian beta model" with
; exponent 1.8 .
;
lrfspec = voorsmooth(spec, 5., type='beta=1.8')
; plot the result
;
plot, sp_wgrid,spec,psym=10,xr=[5,8],/xs,/ylog,$
xtitle='Wavelength ['+!AA+']',ytitle='Counts [ct s!u-1!n]',$
title='Predicted spectrum'
oplot, sp_wgrid,lrfspec,color=1,psym=10,thick=2
;NOTE:
; these variables have been saved in !ARDB/example_1.save thusly:
;save,file=!ARDB+'/example_1.save',sp_wgrid,spec,lrfspec,linspec,conspec
; restore stored values to continue with current demo
;
restore,!ARDB+'/example_1.save',/verbose
; select the lines of interest ; ae=where(rd_line_Z eq 10) ; make the spectrum (note that using the full wavelength range ; of sp_wgrid is a bit inefficient here though) ; linspec_Ne = hastogram(abs(rd_line_WVL[ae]), sp_wgrid, wts=lineflx_flx[ae]) ; convolve with the "LRF" ; lrfspec_Ne = voorsmooth(linspec_Ne, 5., type='beta=1.8') ; plot the result ; !p.multi=[0,1,2] plot, sp_wgrid,linspec_Ne,psym=10,xr=[9,14],/xs,$ xtitle='Wavelength ['+!AA+']',ytitle='Counts [ct]',$ title='Predicted Ne spectrum' oplot, sp_wgrid,lrfspec_Ne,color=1,psym=10,thick=2 ; zoom in on Ne IX triplet plot, sp_wgrid,linspec_Ne,psym=10,xr=[13.3,13.8],/xs,$ xtitle='Wavelength ['+!AA+']',ytitle='Counts [ct]',$ title='Ne IX triplet' oplot, sp_wgrid,lrfspec_Ne,color=1,psym=10,thick=2 !p.multi=0
; find locations of lines of a given element and ionisation stage, ; say, all lines of Ne. ae=where(rd_line_linstr.z eq 10) show_line, abs(rd_line_wvl[ae]), Z=rd_line_Z[ae], ion=rd_line_ion[ae], $ lambda=data_wvl, spec=data_spec, order='1,3', sep=',', xo=xo, yo=yo, $ oxr=oxr, oyr=oyr, title='Looking for Neon' ;NOTE: ; SHOW_LINE calls the "general utility" routine PICKRANGE, ; which is described elsewhere. To test out features of ; PICKRANGE, try zooming in on a rectangular region drawn with ; the click-and-drag of the right mouse button, followed by ; left button click to enact the zoom; when an interesting ; feature has been found, get cursor control with a left button ; click, then dump a postscript hardcopy by typing h (can ; toggle color, encapsulated postscript, or full landscape ; by typing c, e, and l), then = to dump the file.
; To get a rough first idea of line ID's for prominent lines, ; just select lines that are stronger than 1/10 the flux of ; the strongest line ae=where(lineflx_flx/max(lineflx_flx) gt 0.1) show_line, abs(rd_line_wvl[ae]), Z=rd_line_Z[ae], ion=rd_line_ion[ae], $ lambda=data_wvl, spec=data_spec, order='1,2,3', sep=',', xo=xo, yo=yo, $ oxr=oxr, oyr=oyr, title='Expected strong lines' ;NOTE: ; continuing to test PICKRANGE: try, e.g., browsing through ; the spectrum by zooming in on a narrow wavelength range to ; the short wavelength end (right button click and drag, then ; click with left button to select), then click again with ; left button to get keyboard control; now > and < scroll ; through the spectrum to right and left.
; One more: Fe XVII lines, but only the stronger ones again, ; and scaling the y position of the labels according to the ; theoretical fluxes from lineflx, using ~200 as peak of ; Fe XVII resonance line at 15.015 AA, and adding to very ; rough continuum derived above (section 2). aei=where(rd_line_Z eq 26 and rd_line_ION eq 17) ae=where(rd_line_Z eq 26 and rd_line_ION eq 17 and $ lineflx_flx/max(lineflx_flx[aei]) gt 0.05) in_x=abs(rd_line_wvl[ae]) in_y=200.*lineflx_flx[ae]/max(lineflx_flx[ae]) + $ interpol(conspec,0.5*(sp_wgrid[1:*]+sp_wgrid),in_x) ; show_line, in_x, in_y, Z=rd_line_Z[ae], ion=rd_line_ion[ae], $ lambda=data_wvl, spec=data_spec, order='1,2,3', sep=',', $ xo=xo, yo=yo, oxr=oxr, oyr=oyr, title='Fe XVII' ; Note that zooming in here does not change y label ; location because y positions are fixed by in_y
; find density sensitive lines using LSD. This routine needs ; to read in the emissivities for different densities in order ; to find which lines have changed flux. Set the threshold for ; density sensitivity to factor of 1.5 lsd_list = lsd([!WMIN,!WMAX], lsd_flx, lsd_wvl, edens=[1.e9,1.e13], $ DEM=!DEM, tlog=!LOGT, dbdir='$SPEX', ceiling=5., flor=1.e-2, $ ratmax=lsd_ratmax, flxmax=lsd_flxmax, outZ=lsd_Z, outIon=lsd_Ion, $ chifil=1, chidir=!CHIDIR, eqfile=!IONEQF, abund=!ABUND, $ effar=effar, wvlar=wvlar) show_line, lsd_wvl, Z=lsd_Z, ion=lsd_ION, lambda=data_wvl, spec=data_spec, $ order='1,2,3', sep=',', xo=xo, yo=yo, oxr=oxr, oyr=oyr, $ title='Density Sensitive Lines'
; find lines that are relatively unblended, say less than 20%. ; First set the "resolution" - the +/- delta wavelength from ; line centre within which to consider neighbouring lines as ; blended. ; in_dwvl=0.015 ;in [AA], reasonable rough number for MEG ; figure out the blending factor for all the lines ; bland_frac = bland(rd_line_linstr, in_dwvl, flx=lineflx_flx, $ wflx=abs(rd_line_wvl)) ; choose the "least blended" ; ae=where(bland_frac gt 0.8 and lineflx_flx/max(lineflx_flx) gt 0.1) ; and display them show_line, abs(rd_line_wvl[ae]), Z=rd_line_Z[ae], ion=rd_line_ION[ae], $ lambda=data_wvl, spec=data_spec, order='1,2,3', sep=',', xo=xo, yo=yo, $ oxr=oxr, oyr=oyr, title='Unblended Lines'
; Now run SMUDGE to find instrument features and use in ; SHOW_LINE too, first selecting the instrument configuration ; (this is accomplished by logical contructions using keywords ; keys (OR), only (AND) and nisi (NOT) in SMUDGE - yes, it is ; not so obvious but was done that way to avoid a huge and ; repetitive instrument information database file ; smudge, smudge_loc, smudge_label, keys='HRMA HEG', range=[!WMIN,!WMAX], $ infile=!ARDB+'/smudge.lst' ; and display them _along with_ the unblended lines show_line, abs(rd_line_wvl[ae]), Z=rd_line_Z[ae], ion=rd_line_ION[ae], $ lambda=data_wvl, spec=data_spec, wmark=smudge_loc, lmark=smudge_label, $ order='1,2,3', sep=',', cmark=150, markc=222, xo=xo, yo=yo, oxr=oxr, $ oyr=oyr, title='Unblended Lines + Instrument Features' ;NOTE: uncorrected bug -- first two unblended lines seem to get ; label locations forcibly set to 0.
lineid_idstr = lineid(data_wvl, data_spec, stor=rd_line_linstr, $ n_e=!EDENS, dbdir=!LDBDIR, chifil=1, eqfile=!IONEQF, chidir=!CHIDIR, $ DEM=!DEM*EXPTIM_V711TAU, abund=!ABUND, effar=effar, wvlar=wvlar, $ verbose=!VERBOSE) ;NOTE 1: ; since RD_LINE_LINSTR already exists with the correct emissivities, ; there is no need to read them all in all over again, which LINEID ; does automatically otherwise. The RD_LINE()/FOLD_IONEQ() keywords ; such as N_E, DBDIR, EQFILE, etc. are used only if STOR is not ; properly defined. ;NOTE 2: ; The sample IDs are in an IDL save file, at !ARDB/example_2.save: ;save,file=!ARDB+'/example_2.save',lineid_idstr ; restore stored values to continue with current demo ; restore,!ARDB+'/example_2.save',/verbose
; need do only the following to restore all necessary variables ; (uncomment if cold start) ; ;.run initale ;restore,!ARDB+'/demo_v711Tau.save' ;restore,!ARDB+'/example_2.save' ; extract some useful information out of the ID structure ; first, the number of independent line features -- ; fitlines_pos=lineid_idstr.WVL & ncomp=n_elements(fitlines_pos) ; then, the expected fluxes -- ; fitlines_flx=fltarr(ncomp) for i=0L,ncomp-1L do fitlines_flx[i]=total(lineid_idstr.(i+1).FLUX) ; set the widths to some reasonable number ; fitlines_wdt=0.01 ; choose a function type for the line profile ; fitlines_type='beta=1.8' ; make labels for each feature, so that it will be easy to keep track ; fitlines_epithet=strarr(ncomp) all_labels = idlabel(lineid_idstr, j, /Wbeda) for i=0L,n_elements(j)-1L do $ fitlines_epithet[j[i]]=fitlines_epithet[j[i]]+all_labels[i]+' ' ; fit ; fitlines_fitstr = fitlines(data_wvl, data_spec, ysig=data_spec_err, $ pos=fitlines_pos, flx=fitlines_flx, wdt=fitlines_wdt, $ perrp=fitlines_perrp, perrm=fitlines_perrm, perrc=fitlines_perrc, $ ferrp=fitlines_ferrp, ferrm=fitlines_ferrm, ferrc=fitlines_ferrc, $ werrp=fitlines_werrp, werrm=fitlines_werrm, werrc=fitlines_werrc, $ thaw=fitlines_thaw, type=fitlines_type, epithet=fitlines_epithet, $ ties=fitlines_ties, conlev=fitlines_conlev, consig=fitlines_consig, $ /dumb, verbose=!verbose) ;NOTE: ; the state of the fit is in the IDL save file (dumped from ; within the GUI) !ARDB+'example_3.save' ; the output is in an IDL save file ;save,file=!ARDB+'/example_4.save',fitlines_fitstr ; restore stored values to continue with current demo ; restore,!ARDB+'/example_4.save',/verbose
; need do only the following to restore all necessary variables ; (uncomment if cold start) ; ;.run initale ;restore,!ARDB+'/demo_v711Tau.save' ;restore,!ARDB+'/example_2.save' ;restore,!ARDB+'/example_4.save' ; here we will use the measured fluxes for the identified lines ; in conjunction with the emissivity profiles to derive upper limits ; to the emission-measure distribution at each temperature. potem_tool_EM = potem_tool(lineid_idstr, fitlines_fitstr.flx, logT=!LOGT, $ abund=!ABUND, multid=potem_tool_multid, effar=effar, wvlar=wvlar, $ NH=!NH, He1=!He1, HeII=!HeII, verbose=!verbose) ; plot the output EM curves ; plotem, !LOGT, potem_tool_EM, idstr=lineid_idstr, multid=potem_tool_multid, $ xr=[6,8], lsiz=1 ;NOTE: ; the emission measures are in the IDL save file !ARDB+'example_5.save' ;save,file=!ARDB+'example_5.save',potem_tool_EM,potem_tool_multid ; restore stored values to continue with current demo ; restore,!ARDB+'/example_5.save',/verbose
; requires the outputs of feature identifications and ; line profile fitting. (uncomment if cold start) ; ;.run initale ;restore,!ARDB+'example_2.save' ;restore,!ARDB+'example_4.save' ; first, update the locations of the features and the ; line fluxes according to the measured (fitted) values ; fitlines_flxerr=0.5*(fitlines_fitstr.ferrp+fitlines_fitstr.ferrm) lineid_idstr = updatid(lineid_idstr, fitlines_fitstr.flx, $ fitlines_fitstr.pos, flxerr=fitlines_flxerr) help, cat_id(lineid_idstr) ; now write out the updated information to a LaTeX file ; tekku, lineid_idstr, texfil=!ARDB+'/example' ; the resulting LaTeX file may be compiled and displayed ; cd, '/tmp' spawn, 'latex '+!ARDB+'/example' spawn, 'xdvi example'