A detailed example session

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

  1. read in the line and continuum emissivity databases,
  2. compute a theoretical spectrum,
  3. browse spectra with special attention to lines of specific species, density sensitive lines, and unblended lines, etc.
  4. identify line features in the observed spectrum,
  5. derive line fluxes by fitting line profiles,
  6. compute emission measures appropriate for the derived fluxes, and
  7. create a LaTeX table of the identifications


#	start IDL

;	set up the PINTofALE environment
.run initale

;	if INITALE fails, run the script
;	@PoA_constructor

;	load the data and instrument files

;	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.],$
for i=0L,n_elements(data_wvl)-1L do oplot,$
  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
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, $

;	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)
;	and do not forget to include ion balance again..
;	.. 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)


There are many different programs in PINTofALE (e.g., LINESPEC, LINESPEC_EM, MAKE_SPECTRUM, etc.) that will allow one to generate a theoretical spectrum. The basic scheme however is to first compute observed line fluxes from the line emissivities and bin them into a spectrum, and add the continuum spectrum, derived by computing the fluxes and rebinning to the same wavelength grid as the line spectrum. Below we will illustrate this scheme in detail, first making a complete spectrum, and second exploring smaller parts of the spectrum. Note that this section requires variables defined in the previous section.
  1. Compute Spectrum
  2. Subsets of Spectrum


;	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, $

;	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]

;	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

;	these variables have been saved in !ARDB/example_1.save thusly:

;	restore stored values to continue with current demo


The above will give a "complete" spectrum with all the lines in the database. To compute, e.g., just the spectrum of Ne with no continuum is easy by selecting only the lines from elements with atomic number = 10
;	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
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


This section illustrates some interactive ways to locate ``interesting'' regions and features in an observed spectrum by linking it to the database read in above via theoretically biased filters.
  1. Locate Specific Lines
  2. Explore Strong Lines
  3. Strong Lines of Particular Species
  4. Density sensitive Lines
  5. Unblended Lines
  6. Instrument Features


;	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'

;	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'

;	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_y=200.*lineflx_flx[ae]/max(lineflx_flx[ae]) + $
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, $

;	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], $

;	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.


See http://hea-www.harvard.edu/PINTofALE/doc/LINEID.howto for a detailed description on how to use LINEID(). Here we only present a sample call and a sample set of IDs, for use in examples below. The following calls require that the variables defined in the initialization section are present.
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, $

;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:

;	restore stored values to continue with current demo


See http://hea-www.harvard.edu/PINTofALE/doc/FITLINES.howto for a demo of fitting line profiles. Here, as above, we only present a sample call proceeding from previously derived line identifications, for use in further examples below.
;	need do only the following to restore all necessary variables
;	(uncomment if cold start)
;.run initale

;	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 --
for i=0L,ncomp-1L do fitlines_flx[i]=total(lineid_idstr.(i+1).FLUX)

;	set the widths to some reasonable number

;	choose a function type for the line profile

;	make labels for each feature, so that it will be easy to keep track
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)

;	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

;	restore stored values to continue with current demo


;	need do only the following to restore all necessary variables
;	(uncomment if cold start)
;.run initale

;	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

;	the emission measures are in the IDL save file !ARDB+'example_5.save'

;	restore stored values to continue with current demo


;	requires the outputs of feature identifications and
;	line profile fitting.  (uncomment if cold start)
;.run initale

;	first, update the locations of the features and the
;	line fluxes according to the measured (fitted) values
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'

pintofale()head.cfa.harvard.edu (VK)