|
Giulio Del Zanna (G.Del-Zanna@damtp.cam.ac.uk) |
1 CDS data analysis
1.0.1 Reading a FITS file
file='s11484r06'
a=readcdsfits(file)
Searches the 's11484r06.fits' FITS file in the $CDS_FITS_DATA
directory.
IF found, it loads the data into the so-called quicklook data
structure (QLDS) a
Note that QLDS are 'special' entities CDS-specific, and normal
IDL commands MUST NOT be used. for example, to copy:
qlds = copy_qlds(a)
Once the data have been loaded, the QuickLook package can be used
to inspect (and analyse) the data:
dsp_menu, qlds
1.0.2 Subtraction of the CCD bias
vds_debias, qlds
1.0.3 Cosmic ray removal
An useful GUI is XCDS_COSMIC
which allows to use e.g. CDS_CLEAN_EXP (OK for pre-recovery
not Y-binned data) or CDS_CLEAN_SPIKE
Important NOTE: you will need to have defined
$CDS_EXTERNAL to point to a compiled version of external.f
if you are going to use CDS_CLEAN_EXP.
xcds_cosmic,qlds
1.0.4 Corrections for flat-fielding, burn-in and
nonlinear effects, using VDS_CALIB
vds_calib, qlds, exptime_corr=exptime_corr,THROUGHPUT=THROUGHPUT
; QLDS now has units PHOT-EVENTS/PIXEL/SECONDS
; exptime_corr,THROUGHPUT are outputs I added to the routine
;if you want COUNTS/PIXEL/S multiply for THROUGHPUT
;if you want PHOT-EVENTS/PIXEL multiply by exptime_corr
1.0.5 Radiometric Calibration
There are many ways of dealing with this issue.
I recommend to analyse the 'raw' data, and apply the
radiometric Calibration only at the end.
One possible way to do this is to create a 'dummy'
QLDS, fill the data with unity values, and then
apply the Radiometric Calibration.
In this way, you will have arrays with the correction
factors, that you can inspect and eventually modify.
qlds_calib=copy_qlds(qlds)
NWIN = N_ELEMENTS(QLDS_CALIB.DETDESC)
; fill the data with unity values:
;----------------------------------
FOR IWIN = 0,NWIN-1 DO BEGIN &$
ARRAY = DETDATA(QLDS_CALIB, IWIN+1) &$
ARRAY[*]=1. &$
ST_WINDATA,QLDS_CALIB,IWIN,ARRAY,/NO_COPY &$
ENDFOR
; apply the calibration:
;-----------------------
nis_calib,qlds_calib, /ERGS, /STERAD
; The default units for the returned data are
; PHOTONS/CM^2/SEC/ARCSEC^2, but
; since we used the /ERGS, /STERAD keywords, the
; unints are ERGS/CM^2/SEC/STERAD
;------------------------------------------------
1.0.6 Which data windows are available ?
lw = gt_wlimits(qlds, /short)
Window Label Band Xmin Xmax Wmin Wmax
0 O_4_554_51 NIS2 344 368 553.13 555.93
1 HE_1_584_33 NIS2 599 623 582.94 585.75
2 O_3_599_60 NIS2 729 753 598.17 600.99
3 MG_9_368_07 NIS1 850 874 367.22 368.91
4 MG_10_624_94 NIS2 945 969 623.52 626.34
5 O_5_629_73 NIS2 986 1010 628.34 631.16
The Xmin Xmax numbers give the detector start and stop pixels for each
data window.
The wavelengths are calculated using an average
pixel-wavelength conversion applicable to the date of observation.
The table above also shows that the data windows are 25 pixels wide.
LINES TO FIT:
--------------
Window=0
O_IV_553.5
O_IV_554.5
O_IV_555.5
Window=1 (not necessary)
He_I_584
Window=2
O_III_599
Window=3
Mg_VII_367
Mg_IX_368
Window=4
Mg_X_625
O_IV_625.8
Window=5
O_V_629
; If you are interested in comparing actual timings for the
; events as seen by different instruments, you need to
; get the starting time for each of the exposures:
time_start=gt_start(QLDS,/exposure,/truncate)
1.0.7 Put the data of each NIS windows into an array
;get the data for the first window:
;---------------------------------
window=0
data = gt_windata(qlds,window)
;multiply the data for the corrected exposure time to get PHOT-EVENTS/PIXEL:
; also, make sure that missing pixels are flagged.
print,min(data)
-100.000
missing=-100.
bad=where(data eq missing,nbad)
data =data * EXPTIME_CORR
if nbad gt 0 then data[bad]=missing
help,data
DATA FLOAT = Array[25, 10, 73]
The wavelength dimension has 25 pixels, 10 exposures were taken
by stepping in SolarX, and only 73 pixels along the slit (SolarY)
were telemetred.
In general, the variable data will then be a data array with from 2 to
4 dimensions. The
dimensions are always in the order (Wavelength , SolarX,
SolarY, Time).
The fourth dimension (Time) is only added where the observation sits at
the same SolarX location and takes repeated exposures. In
that case the SolarX dimension is retained, but only has
one element.
1.0.8 Get wavelengths
There are various ways. The simplest is to use:
dummy=gt_spectrum (qlds,WINDOW=window,xix=0,yix=0,lambda= w)
1.0.9 Get average spectra and define the fitting structure
; average over the Solar X,Y:
sp=average(average(data, 2, missing=-100),2, missing=-100)
window,0 & plot, w,sp,psym=10
delvarx,fit
xcfit, w,sp, fit
add components, then EXIT.
Adjust CFIT components. E.g. adjust min, max Width:
help,fit,/st
** Structure <883fb14>, 4 tags, length=1592, refs=1:
IGAUSS4 STRUCT -> COMPONENT_STC_3 Array[1]
IGAUSS3 STRUCT -> COMPONENT_STC_3 Array[1]
IGAUSS2 STRUCT -> COMPONENT_STC_3 Array[1]
BG STRUCT -> COMPONENT_STC_1 Array[1]
fit.igauss2.param[2].MIN_VAL=0.4
fit.igauss2.param[2].MAX_VAL=0.8
fit.igauss3.param[2].MIN_VAL=0.4
fit.igauss3.param[2].MAX_VAL=0.8
fit.igauss4.param[2].MIN_VAL=0.4
fit.igauss4.param[2].MAX_VAL=0.8
;fit.(2).name='O_IV_553.5'
;fit.(1).name='O_IV_554.5'
;fit.(0).name='O_IV_555.5'
Note: all the minimum, maximum values MUST be adjusted.
Also, all the names should be uniquely defined.
1.0.10 Do the fit and create image maps
I've written a simple routine that takes as
input the data, performs multiple fitting for each NIS window
at a time, and
outputs Dominic Zarro's maps.
; Do the fit:
;Input:
; file: file name.
; QLDS, QLDS_CALIB: QLDS as previously obtained.
; WINDOW: the integer of the NIS window
; EXPTIME_CORR: the corrected exposure time
; MISSING: the value of the missing pixel
cds_fit_win, file, QLDS, QLDS_CALIB, WINDOW, EXPTIME_CORR, MISSING
; at the end, DELETE the QLDS from memory:
delete_qlds,qlds
2 Obtain spatially-averaged intensity values and derive
EM, Te, Ne
;first restore the intensity maps
;-----------------------------------------
restgen, file='s11484r06_O_V_cds_map.genx',st=cds_map
intensity=cds_map.data
; units are ergs cm-2 s-1 sr-1
; select some pixels of interest:
;-----------------------------------------
; magnification factor:
mag=10
sz=size(intensity)
expand = rebin(intensity, mag*sz[1], mag*sz[2], /sample)
wdim = size(expand)
window,1,xpos=1275-wdim(1),ypos=5,xsize=wdim(1),ysize=wdim(2)
tvscl, expand
roi = defroi(sz[1], sz[2],zoom=mag)
ij = get_ij(roi,sz[1])
x = reform(ij[0, *])
y=reform(ij[1, *])
; GET averaged intensity values:
;-----------------------------------------
print, average(intensity[x,y],missing=missing)
; [Write down values -
;-----------------------------------------------------------
;-----------------------------------------------------------
; obtain the G(T) functions from the CHIANTI software, e.g.:
;-----------------------------------------------------------
; [default units are erg cm^+3 s^-1 sr-1 ]
; Here you have to select the G(T) corresponding to the lines
; for which you have intensities:
gofnt, 'o_3', 580, 635, t_o_3, g_o_3, d_o_3, $
press=1e15, abund_name=!xuvtop+'/abundance/sun_coronal.abund', $
ioneq_name=!xuvtop+'/ioneq/mazzotta_etal.ioneq'
gofnt, 'o_4', 550, 635, t_o_4, g_o_4, d_o_4, $
press=1e15, abund_name=!xuvtop+'/abundance/sun_coronal.abund', $
ioneq_name=!xuvtop+'/ioneq/mazzotta_etal.ioneq'
gofnt, 'o_5', 580, 635, t_o_5, g_o_5, d_o_5, $
press=1e15, abund_name=!xuvtop+'/abundance/sun_coronal.abund', $
ioneq_name=!xuvtop+'/ioneq/mazzotta_etal.ioneq'
gofnt, 'mg_9', 300, 380, t_mg_9, g_mg_9, d_mg_9, $
press=1e15, abund_name=!xuvtop+'/abundance/sun_coronal.abund', $
ioneq_name=!xuvtop+'/ioneq/mazzotta_etal.ioneq'
gofnt, 'mg_10', 580, 635, t_mg_10, g_mg_10, d_mg_10, $
press=1e15, abund_name=!xuvtop+'/abundance/sun_coronal.abund', $
ioneq_name=!xuvtop+'/ioneq/mazzotta_etal.ioneq'
;-----------------------------------------------------------
; Obtain emission measure loci curves (N_e*N_H*dh [cm^-5])
;-----------------------------------------------------------
F_o_3= 268./g_o_3 ; 599.59
F_o_4= 168./g_o_4 ; 553.37
F_o_5= 1812./g_o_5 ; 629.7
F_mg_9= 482./g_mg_9 ; 368.
F_mg_10= 97. /g_mg_10 ;625.
;-----------------------------------------------------------
; Plot them
;-----------------------------------------------------------
plot, alog10(t_o_3), alog10(F_o_3),$
xr=[4.5, 7.], yr=[25, 28], /yst, /xst, $
xtit='Log T [K]', ytit='Log F(T) [cm-5]', chars=1.3
oplot, alog10(t_o_4 ) , alog10(F_o_4 )
oplot, alog10(t_o_5 ) , alog10(F_o_5 )
oplot, alog10(t_mg_9 ) , alog10(F_mg_9)
oplot, alog10(t_mg_10 ) , alog10(F_mg_10 )
; x2ps, 'ft.ps'
;-----------------------------------------------------------
; Obtain a value of Ne by using
dens_plotter,'o_4'
; and the O IV intensities you have obtained.
;-----------------------------------------------------------
3 Overplot CDS, TRACE, MDI maps
;!p.background=0
;!p.color=255
loadct,0
set_line_color
black=0 & white=1 & yellow=2 & red=3 & green=4 & blue=5
orange=6 & purple=7 & magenta=8 & brown=9 & turquoise=10
; define the CDS FOV and Centre, e.g. (see help,CDS_MAP,/st) :
;-------------------------------------------------------------
fov =[0.8, 2.2] & center=[-218., 250.]
; plot the CDS map:
;------------------
window,0
plot_map, CDS_MAP, fov=fov,center=center, $ ; dmin=1,dmax=800,$
trans=[0,0], missing=-100.
; read a corresponding TRACE image, e.g.:
;----------------------------------------
read_trace, 'tri19980619.1200_0263.pl' , -1,index,data
;create TRACE map:
;-----------------
index2map, index,data, tracemap
tracemap = map2soho(tracemap)
; plot the TRACE map:
;--------------------
window,1
plot_map, tracemap, trans=[0,0],fov=fov,center=center
;rebin the TRACE map down to a resolution comparable to CDS:
;-----------------------------------------------------------
dummy=rebin(tracemap.data, 32, 80) ; factor of 8 binning
map=rep_tag_value(tracemap, dummy, 'data')
map.dx=8*map.dx & map.dy=8*map.dy
; plot the re-binned TRACE map:
;------------------------------
window,2
plot_map, map, trans=[0,0],fov=fov,center=center
; plot the re-binned TRACE map as overlay of the CDS map:
;--------------------------------------------------------
window,3
plot_map, CDS_MAP, fov=fov,center=center, $ ; dmin=1,dmax=800,$
trans=[0,0], missing=-100.
plot_map, map, trans=[0,0],/over, /cont
; NOTE: you will need to adjust the TRANS values.
;------------------------------------------------
; read a corresponding MDI magnetogram:
;--------------------------------------
rd_mdi,'hr_M_1024x500_01h.47892.0059.fits',$
index_mdi,data_mdi
; creat an MDI map:
;------------------
index2map,index_mdi,data_mdi,mdimap
; plot MDI map:
;---------------
window,4
plot_map, mdimap, trans=[0,0],fov=fov,center=center,$
dmin=-50,dmax=50
; plot MDI map overlayed on the TRACE map, with colour contours:
;---------------------------------------------------------------
window,5
plot_map, tracemap, charsize=1.4,margin=0.1, $
trans=[0,0],fov=fov,center=center, dmin=90, dmax=500
plot_map, mdimap, charsize=1.4,margin=0.1, $
trans=[0, 0], fov=fov,center=center,$
bottom=12,lev=[-20, -30,-40, -50, -100], lcol=green, /over, /cont
plot_map, mdimap, charsize=1.4,margin=0.1, $
trans=[0, 0], fov=fov,center=center,$
bottom=12,lev=[20, 30,40, 50, 100], lcol=red, /over, /cont
; NOTE: you might need to adjust the TRANS values.
;------------------------------------------------
; plot MDI map overlayed on the CDS map, with colour contours:
;-------------------------------------------------------------
window,6
plot_map, CDS_MAP, fov=fov,center=center, $ ; dmin=1,dmax=800,$
trans=[0,0], missing=-100.
plot_map, mdimap, charsize=1.4,margin=0.1, $
trans=[0, 0], fov=fov,center=center,$
bottom=12,lev=[-20, -30,-40, -50, -100], lcol=green, /over, /cont
plot_map, mdimap, charsize=1.4,margin=0.1, $
trans=[0, 0], fov=fov,center=center,$
bottom=12,lev=[20, 30,40, 50, 100], lcol=red, /over, /cont
; NOTE: you might need to adjust the TRANS values.
;------------------------------------------------
; End !
;******
File translated from
TEX
by
TTH,
version 3.33.
On 4 Nov 2003, 18:54.