How to use Pint of Ale from the generation of a synthetic spectrum to the DEM creation
Kerry Neal


1. Initialize PINTofALE by .r initale

2. Call synthetic.pro, at first I just input a single ion so that I knew what I was looking for in the spectrum.

synthetic, 150., 250., 0.1, lambda, spectrum, list_wvl, list_ident, density=5.e8, /photons, sngl_ion=’Fe_10’

Chose mazzotta_etal_ext.ioneq for Ionization Equilibrium File, quiet_sun.dem for DEM input and sun_coronal_ext.abund for the Abundance File The spectrum needs to be an integer, not a double array, thus to change

spec=long64(spectrum)

3. Run lineid:

idstr=lineid(lambda, spec, stor=rd_line_linstr, n_e=!EDENS, dbdir=!LDBDIR, chifil=1, eqfile=!IONEQF, chidir=!CHIDIR, abund=!ABUND, effar=effar, wvlar=wvlar, verbose=!verbose)

Lineid takes awhile to run, it will spit out a GUI, which did not work for my computer according to the documentation. Right click and drag to select an area Left click and then right click to zoom in on the area selected Highlight a peak by using the left button, a line is drawn over the peak. If a line extends over multiple peaks then the highest one will be used in the next step. Right click then left click, another GUI pops up with a line list, usually the element that you want is the first entry, if the first entry has a negative value, then you want the second one. Click NEXT, then you are back to the graph, to zoom out click both right and left or the middle mouse button. If you have identified the last line instead of clicking NEXT click QUIT and it will save the line you just selected and all the ones previous. Back at the command line help, idstr, /struc and you will see a wvl tag and ID tags for each wavelength.

4. The next step is to use fitlines to fit Gaussians to the peaks to establish the flux and the flux errors. I have found this GUI to be rather touchy and of the two methods used in the fitlines.howto, I have used the second one, which is more manual. I will describe that one here. Unfortunately, the chi squared values have been extremely high, although the Gaussians that fit do not seem to fit well, in my case they are taller and skinnier than the actual peak, I can see that if you shortened the peak and made it wider, it would represent the peak pretty well. The integral under each curve, the actual and the fit, is in actuality pretty close to one another.

First, define the average position and the number of peaks in the idstr: pos=fltarr(3)+177 (! Watch out ! Do NOT put a period after the 177, this changes the position from an integer to a double and the program will not recognize double arrays, thus if you put a period the graph that comes up when you call fitlines will not appear because it cannot recognize the position of any of the lines)

type=’beta 1.8’
wdt=0.05
fstr=fitlines(lambda, spec, pos=pos, perrp=perrp, perrm=perrm, perrc=perrc,$ flx=flx, ferrp=ferrp, ferrm=ferrm, ferrc=ferrc,wdt=wdt, $ werrp=werrp, werrm=werrm, werrc=werrc, thaw=thaw, type=type, epithet=epithet, ties=ties,$ conlev=conlev,consig=consig,/dumb,verbose=!verbose)

A GUI will come up. First, click VIEW--> SHOW, makes sure that a graph comes up with the input spectrum, if it doesn’t see if you defined the pos input variable as a double array instead of an integer array, a tiny, but fatal detail.

Then define the components on the left hand side, in the first box entitled Component put 0 for the first component, 1 for the second, etc. Position is the wavelength of each component. Press enter after you enter the position. Then fill in the last box, Label, with the name of the line Fe10 a for example. Then press enter again. Repeat this process for each component.

Next, fit each line individually but going to the XRANGE --> KEYBOARD, enter a small range of 0.2-0.4 Angstroms wide to see just the peak. Input the range at the IDL prompt line, it will say X-Range: >.

Then, select YRANGE --> CURRENT, then VIEW --> SHOW.

If you are happy with this then click MODEL --> ADD COMPONENT.

On the graph, click with the left button on the top of the peak, then drag down to the right or left side of the peak to represent the width of the peak. A red line will outline the fit you just chose. Right click and then left click on the graph to exit.

Go back to the MODEL and choose FREEZE UNIMPORTANT, which will freeze the other peaks. Click on the RENORM button, just above the Component box. Then VIEW.

Next define the continuum.

CONTINUUM --> PIECEWISE Go over to the graph and left click below the x-axis to get rid of Poisson noise. Then click with the left button somewhere a little above zero and drag vertically downwards for just a bit, maybe 4-5 notches on the y-axis.

To exit, once again right click and then left click.

To accept the continuum, go to, CONTINUUM -->ACCEPT.

Then RENORM SHOW FIT, the chi squared value appears after you fit.

VIEW

Repeat this from the top, by redefining the XRANGE. If you want to fit another component just go back to the MODEL, ADD COMPONENT and continue from there. Fit each line.

When you are finished SAVE, then DUMP and at the IDL prompt line create a .sav file. Once this is done QUIT.

 Save the fstr once done with the GUIs.

Now it’s time to prepare the data for the DEM analysis I found these next few steps to be somewhat difficult because you need to figure out which components are the ones that you want, if you fitted a curve more than once then there are multiple read outs for the positions, fluxes and flux errors.

In the Pint of Ale example it says: flxerr=(fstr.ferrp+fstr.ferrm)/2 idstr=updatid(idstr, fstr.flx, fstr.pos, flxerr=flxerr)
I found it necessary though to print out the arrays of fstr.ferrp and fstr.ferrm to compare with the fstr.pos, to find which fits were the ones that I wanted to keep. Then I just redefined these arrays such as:
Flx=[fstr.flx[0], fstr.flx[3], fstr.flx[5]]
Then I put those arrays into the commands above.
Next step is to consolidate multiple IDs:
idstr=squishem(idstr)
Define wavelengths, fluxes, flux errors and emissivities:
wavelengths=idstr.wvl
fluxes=[idstr.id1.flux, idstr.id2.flux, …]
fluxerr=[idstr.id1.fluxerr, …]
 z=[idstr.id1.z, …]
emissivity=dblarr(81, ncomp of wavelengths)
emissivity[*,0]=idstr.id1.emis[*]
emissivity[*,1]=idstr.id2.emis[*]
.
.
.
newlogt=(findgen(num of temp elements you have)+blah)/10

Line Ratios!
Ratio1=fluxes[2]/fluxes[4]
Ratio2=(fluxes[5]-fluxes[3])/(fluxes[2]+fluxes[3]+fluxes[5])
onlyrat=[‘’,’’, ‘+N1, +D2’, ‘-N2’, ‘+D1,+D2’, ‘+N2,+D2’]

Find where there is no line, just noise for this next part:
plot, lambda, spec, xr=[noise1, noise2]
ll=where(lambda ge noise1 and lambda le noise2]
wgrid=mid2bound(lambda)
lrf=0.02 & type=’beta=2.5’
licospec, wgrid[ll], lspec, cspec, verbose=10, lstr=lstr, cstr=cstr, n_e=1e10, lrf=lrf, type=type, lcthr=0.05, constr=constr
plot, wgrid, lspec/cspec, /ylog
oo=where((lspec/cspec) lt 0.05)
tmp=spec[ll]
contflx=total(tmp[oo])
conemiso=constr.cont_int*(max(wgrid)-min(wgrid))
conemis=rebinx(conemiso, !logt, newlogt)
fluxes=[fluxes, contflx]
fluxerrs=[flxerr, 1+sqrt(0.75+contflx)]
wavelengths=[wavelengths, middle of noise1 and noise2]
emissivity=dblarr(81, # wavelength components now)
emissivity[*,0]=idstr.id1.emis[*]
.
.
.
newemis=rebinx(emissivity, !logt, newlogt)
newemis[*,last elem of array]=conemis[*]

Ready!
adem=mcmc_dem(wavelengths, fluxes, newemis, z=z, logt=newlogt, fsigma=flxerr, onlyrat=onlyrat, simdem=simdem, demerr=demerr, simprb=simprb, savfil=’file.sav’)