Example session using .par files and less explicit plotting etc

This annotated example walks a user through a PINTofALE session using predominantly .par files where available. It is similar to the explicit session, but is designed to be simpler, less wordy, less explicit plotting, etc. During this run, the user will

1. read in line and continuum emissivity databases,
2. compute a theoretical spectrum,
3. spectral browsing using SHOW_LINE
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

0. COPY .pars, INITIALIZE SETTINGS, AND LOAD DATA

# copy all the .pars from the ardb directory into working directory - 
# in practice one can do this as each .par is needed.

cp /foo/bar/PoA/ardb/*.par ./.

#	start IDL
idl

;	set up the PoA environment

.run initale

;NOTE:
;	if INITALE fails, run the script
;	@PoA_constructor

;	load the data and instrument files
restore,!ARDB+'/demo_v711Tau.save'
;	put data into standard variables for observed spectrum
;	used in .par files, and then plot them; first counts spectrum 

data_wvl=LAM_M1P_V711TAU & data_spec=SPC_M1P_V711TAU
wvlar=WVLAR_M1P & effar=EFFAR_M1P
!LOGT=LOGT_V711TAU & !DEM=DEM_V711TAU

plot,data_wvl,data_spec,psym=10,xrange=[5.,20.]
plot,wvlar,effar,xr=[5.,20.]		; effective area
plot,!LOGT,!DEM,psym=10,/ylog	; Rough DEM from EUVE Analysis

1. READ IN THE LINE AND CONTINUUM DATABASE

;	set the wavelength range of interest 

!wmin = min(data_wvl)
!wmax = max(data_wvl)

;	read in the line emissivity database; this .par writes, by default,
;	emissivities into out_linemis an id structure into out_linstr
;	Of course, these names can be edited at will in the .par
;	[1e-23 ergs cm^3 s^-1]

@rd_line.par

;	fold in line emissivities with ion-balance info
;	(needed for CHIANTI and SPEX emissivities that are stored
;	independently from ion balance)

@fold_ioneq.par

; put the corrected emissivities into the line emissivity structure
; note that rd_line_linemis will still have emissivites without ion
; balance folded in.

rd_line_linstr.LINE_INT = fold_ioneq_linemis

;	Read in two sets of densities and manually interpolate to
;	desired density.  Making interpolation weights is same as for
;	explicit example; for the rest, just remember to change !edens
;	to the values required, and remember to change back to desired default
;	value.  NOTE: in practise this should not be necessary unless 
;	synthesis of very density-sensitive features if required at a 
;	very specific density.  Even then, this interpolation is only 
;	linear in log10.

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)

!edens=n_e1
@rd_line.par
rd_line_linemis1=rd_line_linemis
!edens=n_e2
@rd_line.par
rd_line_linemis2=rd_line_linemis
!edens=n_e

rd_line_linemis=10.^(weight1*alog10(rd_line_linemis1>1e-30)+$
  weight2*alog10(rd_line_linemis2>1e-30))

; fold ion balance again and put into structure - again, note that
; rd_line_linemis will still have emissivites without ion balance folded
; in

@fold_ioneq.par
rd_line_linstr.LINE_INT = fold_ioneq_linemis

;	read in the continuum emissivity database; ion balance already
;	included in these; will be on wavelength grid out_cwgrid
;	[1e-23 ergs cm^3 s^-1 Ang^-1]

@rd_cont.par

2. COMPUTE THEORETICAL SPECTRUM

;NOTE:
;	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 carried out 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.  This scheme is shown below using some .par's,
;	though the wavelength grid first needs to be explicitly defined.

;	first define a wavelength grid to hold the spectrum.  Note
;	that one could simply use the grid of the observed spectrum,
;	though it is useful to then get a grid of the bin boundaries
;	since these are used by programs such as HASTOGRAM.
;	eg  
;  data_dwvl=data_wvl[1:*]-data_wvl
;  data_wvlb=[data_wvl-data_dwvl/2. , data_wvl[n_elements(data_wvl)-1] $
;             + data_dwvl[n_elements(data_dwvl)-1]/2]

sp_NBIN=20000L	;number of bins
sp_DW=(!WMAX-!WMIN)/sp_NBIN	;bin width
sp_wgrid=findgen(sp_NBIN+1L)*sp_DW+!WMIN

;	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]
;NOTE:  could use lineflx+ismtau.par here to do flux and ism
;       absorption all at once.  Note, in this .par example we multiply the
;	combined line+continuum spectrum by the exposure time right at
;	the end, unlike in the explicit example.

@lineflx.par

;	bin the line fluxes into a spectrum.  hastogram is an improved
;	general histogram routine, much faster than the IDL histogram

linspec=hastogram(abs(rd_line_WVL),sp_wgrid,wts=lineflx_flx)

;	compute continuum fluxes (is actually a flux density), then multiply
;	by bin width. contflx.par actually calls 
;	lineflx, but works consistently for continuum too.  Output 
;	in contflx_flx

@contflx.par
rd_cont_CDW=rd_cont_wvl[1:*]-rd_cont_wvl
contflx_flx=contflx_flx*rd_cont_CDW

;	rebin to defined grid sp_grid - another completely general PoA utility

conspec=rebinw(contflx_flx,rd_cont_wvl,sp_wgrid,/perbin)

;	add the line and continuum spectra

spec=linspec+conspec

;	correct for ISM absorption.  Remember ISM parameters are set
;	up in variables !NH,!FH2,!HE1,!HEII that are used in
;	ismtau.par and output is written to out_tau
;	Need to set in_wvl variable to wavelength grid

in_wvl=sp_wgrid
@ismtau.par
spec=spec*exp(-ismtau_tau)

;	Now convert to total counts and plot the result - without axis titles!

spec=spec*EXPTIM_V711TAU

plot,sp_wgrid,spec,psym=10,xr=[5.,30.],/xs,yr=[1e-2,5e3],/ys,/ylog
setkolor,'red',1 & setkolor,'yellow',2
oplot,sp_wgrid,linspec,color=1,psym=10
oplot,sp_wgrid,conspec,color=2,psym=10

;	convolve with lorentzian of half-width=5*bin-width to simulate LRF
;	This is an example of a "general utility" PoA routine that is
;	not coronal plasma specific.  Use "beta model" with exponent 1.8
 
lrfspec=voorsmooth(spec,5.,type='beta=1.8')

;	plot the result and unconvolved spectrum

plot,sp_wgrid,spec,psym=10,xr=[5,8],/xs,/ylog
oplot,sp_wgrid,lrfspec,color=1,psym=10,thick=2

;NOTE:
;	these variables have been saved in !ARDB/example_1.save:
;save,file=!ARDB+'/example_1.save',sp_wgrid,spec,lrfspec,linspec,conspec

; The above will give a "complete" spectrum with all the lines in the
; database.  To compute, eg, just the spectrum of Ne with no continuum 
; is easy by selecting only the lines from elements with atomic number
; = 10 (note that using the full wavelength range of sp_wgrid is a bit
; inefficient here though).

iNe=where(rd_line_linstr.z eq 10)
linespec_Ne=hastogram(abs(rd_line_WVL[iNe]),sp_wgrid,wts=lineflx_flx[iNe])
lrfspec_Ne=voorsmooth(linespec_Ne,5.,type='beta=1.8')

plot,sp_wgrid,linespec_Ne,psym=10,xr=[9,14],/xs
oplot,sp_wgrid,lrfspec_Ne,color=1,psym=10,thick=2
; zoom in on "triplet"
plot,sp_wgrid,linespec_Ne,psym=10,xr=[13.3,13.8],/xs
oplot,sp_wgrid,lrfspec_Ne,color=1,psym=10,thick=2


; Yet another way to generate a line spectrum from scratch is to use 
; LINESPEC - a routine built up from the individual routines called
; above to arrive at a spectrum.  See linespec.par for the I/O
; used and assumed (eg, !DEM, !ABUND, !WMIN,!WMAX, effar, wvlar etc).
; Output of LINSPEC is a flux density, so need to x by binsize to get 
; flux or count/bin

io_ww=sp_wgrid
@linespec.par

; This should create an output linespec_spec that, provided none of the
; system or other standard variables have been modified, should be the
; same as linspec created above after multiplying by bin width:

plot,linspec,linespec_spec*sp_dw,psym=3,/xlog,/ylog

3. SPECTRAL BROWSING

; find locations of lines of a given element and ionisation stage,
; say, all lines of Ne. - recall iNe defined above by 
; iNe=where(rd_line_linstr.z eq 10)
; First set inputs to showline in_x,in_z,in_ion (latter two define the
; labels and are not required.  Can also specify in_y to set y level
; for labels.  If in_y is not set, labels appear at 2/3 y range.

in_x=abs(rd_line_linstr.wvl[iNe])
in_z=rd_line_linstr.z[iNe]
in_ion=rd_line_linstr.ion[iNe]
@show_line.par

; In show_line.par, a single middle mouse button click will place the
; line labels to their default positions (will also toggle log/linear 
; axis).  Can test out features of
; pickrange by 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 colour, encapsulated
; 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: 

istrong=where(lineflx_flx/max(lineflx_flx) gt 0.1)
in_x=abs(rd_line_linstr.wvl[istrong])
in_z=rd_line_linstr.z[istrong]
in_ion=rd_line_linstr.ion[istrong]
@show_line.par

; still in show_line (actually within pickrange), try, eg, browsing
; through the spectrum by zooming in on a narrow wavelength range to
; the short wavelength end (right button click & 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 from smoothing
; the spectrum

ife17=where(rd_line_linstr.z eq 26 and rd_line_linstr.ion eq 17)
ifexvii=where(rd_line_linstr.z eq 26 and rd_line_linstr.ion eq 17 and $
  lineflx_flx/max(lineflx_flx[ife17]) gt 0.05)
in_x=abs(rd_line_linstr.wvl[ifexvii])
in_z=rd_line_linstr.z[ifexvii]
in_ion=rd_line_linstr.ion[ifexvii]
in_y=200.*lineflx_flx[ifexvii]/max(lineflx_flx[ifexvii]) $
 + interpol(smooth(data_spec,1000),data_wvl,in_x)
@show_line.par

; 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.  This example picks lines whose fluxes
; change by a factor of 5 or more, and that are 1/100 or stronger
; than the strongest density sensitive lines at any of the densities
; investigated.  Do this using SPEX emissivities (but many fewer density
; sensitive lines than CHIANTI)

in_ceiling=5.
in_flor=1.e-2
in_edens=[1.e9,1.e13]
!LDBDIR='$SPEX'
@lsd.par
delvar,in_y
in_x=lsd_wvl
in_z=lsd_z
in_ion=lsd_ion
@show_line.par
!LDBDIR='$CHIANTI'		; re-set !LDBDIR to CHIANTI for later

; find lines that are relatively unblended - less than 20%.  First set the
; "resolution" - the +/- delta wavelength from line centre within
; which to consider neighbouring lines as blended.  Output in blend_frac

in_dwvl = 0.015		; in AA, reasonable rough number for MEG
@bland.par
iblendfree=where(blend_frac gt 0.8 and lineflx_flx/max(lineflx_flx) gt 0.1)
in_x=abs(rd_line_linstr.wvl[iblendfree])
in_z=rd_line_linstr.z[iblendfree]
in_ion=rd_line_linstr.ion[iblendfree]
@show_line.par

; Now run smudge to find instrument features and use in showline too,
; first selecting the instrument config (this is accomplished by
; logical contructions using in_keys, in_only and in_nisi in
; smudge.par - yes, it is not so obvious but was done that way to 
; avoid a huge and repetitive smudge.lst file

in_keys='HRMA MEG'
@smudge.par
@show_line.par

4. IDENTIFY LINE FEATURES

;	See 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.
;	As mentioned in explicit demo, rd_line can of course work with
;	emissivities already read in.  This is controlled by the STOR
;	keyword.  In lineid.par, stor=out_linstr and if out_linstr is
;	a valid structure from rd_line then rd_line and fold_ioneq will not be
;	called.  Otherwise, these will be called and output will be put in 
;	out_linstr.  In this case, calls to rd_line and fold_ioneq are
;	controlled by their keywords in lineid.par

@lineid.par

; Try playing with the features of lineid - eg zoom in on a region
; using the standard PICKRANGE cursor or keyboar commands; select a
; line using a click and drag over the line peak; play with changing
; the parameters in the pop-up GUI listing potential line ID's

;NOTE:
;	The sample IDs are in an IDL save file, at !ARDB/example_2.save:
;save,file=!ARDB+'/example_2.save',out_idstr

5. DERIVE LINE FLUXES

;	see 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 examples below.

;	extract some useful information out of the ID structure
;
; To start from here without going through rest of example session,
; just do
;
;.run initale
;restore,!ARDB+'/demo_v711Tau.save'
;restore,!ARDB+'/example_2.save'

;	set the widths to some reasonable number and 
;	choose a function type for the line profile

fitlines_wdt=0.01
fitlines_type='beta=1.8'

;	fit -  fitlines is very powerful in terms of ability to tie
;	together variables during fitting.  See FITLINES.howto in the 
;	doc directory for more information.  One very useful feature
;	is the ability to set up ties between parameters on calling,
;	so that, eg all the line widths are tied together.
;	Here we are calling lineid_fitlines.par that sets up the line
;	parameters from the ID structure output from LINEID

@lineid_fitlines.par

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

6. COMPUTE EMISSION MEASURES

;.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 line emissivities to derive upper limits
;	to the emission-measure distribution at each temperature.
;	("Carole Jordan type" plot!)

in_flx=fitlines_fstr.flx
@potem_tool.par

;	plot the output EM curves

@plotem.par

;NOTE:
;	the emission measures are in the IDL save file !ARDB+'example_5.save'
;save,file=!ARDB+'example_5.save',out_LOGT,out_EM,out_multid

7. CREATE A LaTeX TABLE OF THE IDs

;	first, update the locations of the features and the
;	line fluxes according to the measured (fitted) values
;	output from updateid is updateid_idstr

@fitlines_updateid.par

;	now write out the updated information to a LaTeX file (this
;	part is identical to the explicit demo except for idstr name)

tekku,updateid_idstr,texfil=!ARDB+'/example'

;	the resulting LaTeX file may be compiled and displayed
cd,'/tmp'
spawn,'latex /example'
spawn,'xdvi example'