-
Notifications
You must be signed in to change notification settings - Fork 3
/
aia_teem_map.old2.pro
362 lines (312 loc) · 12.3 KB
/
aia_teem_map.old2.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
pro aia_teem_map,fileset = fileset,FOV = fov,wave_,npix,teem_table,teem_map, FILElist = filelist, VERBOSE = verbose, filename_extra = FILENAME_extra, SAVE_DIR = save_dir, DEBUG = debug
;+
; Project : AIA/SDO
;
; Name : AIA_TEMPMAP
;
; Category : Data analysis
;
; Explanation : calculates EM and Te temperature maps
; based on single-Gaussian fits in each macropixel
;
; Syntax : IDL>aia_teem_map,fileset,fov,wave_,npix,teem_table,teem_map
;
; Inputs : fileset = common filename part of 6 wavelength FITS images
; wave_ = strarr(6) with wavelengths in Angstroem
; fov(4) = [i1,j1,i2,j2] pixel ranges of field-of-view
; npix = size of macropixel (spatial resolution)
; teem_table = filename of DEM lookup table
; (calculated previously with AIA_TEEM_TABLE.PRO)
; teem_map = savefile containing EM and Te maps
;
; Outputs : postscript file <plotname>_col.ps (if io=2)
;
; History : 3-Mar-2011, Version 1 written by Markus J. Aschwanden
; 16-May-2011, made fov optional. If not given, uses the
; whole images.
; made fileset optional. Added filelist keyword to just give a list of files.
; : 5-Oct-2011 - Added search for 1 sigma contour in chi^2 space to calculate errors in telog and tsig - A. Inglis
;
; Contact : [email protected]
;-
default, save_dir, ''
;_________________________________________________________________________
t1 = systime(0,/seconds)
nwave = n_elements(wave_)
default, filename_extra, ''
IF keyword_set(fileset) THEN BEGIN
files = strarr(6)
FOR iw = 0, nwave-1 DO BEGIN
;searchstring=fileset+'*'+wave_(iw)+'.fits'
searchstring=fileset+'*'+wave_[iw]+'_.fts'
file_iw = file_search(searchstring, count=nfiles)
files[iw]=file_iw[u]
ENDFOR
ENDIF ELSE BEGIN
;need to sort the files in the proper order
files = filelist
;s = fltarr(nwave)
;FOR i = 0, n_elements(wave_)-1 DO s[i] = where(strmatch(files,'*' + num2str(wave_[i]) + '_.fts') EQ 1)
;files = files[s]
ENDELSE
;_____________________REBINNING IMAGES___________________________________
texp_ =fltarr(nwave)
for iw=0,nwave-1 do begin
;read_sdo,files(iw),index,data
fits2map, files[iw], map
IF keyword_set(FOV) THEN BEGIN
sub_map, map, smap, xrange = fov[0:1], yrange = fov[2:3]
map = smap
ENDIF
dim = size(map.data)
; schriste (16-may-2011)
; If Keyword FOV not given then determine the FOV
; of the whole image automatically
;IF keyword_set(fov) THEN BEGIN
; i1=fov[0] & j1=fov[1] & i2=fov[2] & j2=fov[3]
;ENDIF ELSE BEGIN
;s = size(data)
s = size(map.data)
i1=0 & j1=0 & i2=s[1]-1 & j2=s[2]-1
;ENDELSE
;dim =size(data)
nx0 =dim[1]
ny0 =dim[2]
;IF (i1 ge nx0) or (i2 ge nx0) or (j1 ge ny0) or (j2 ge ny0) then begin
; print,'FOV=',i1,i2,j1,j2
; print,'image size=',nx0,ny0
; stop,'ERROR in subimage range i1,i2,j1,j2 for image with size nx0,ny0
;ENDIF
;image = data[i1:i2,j1:j2]
image = map.data
;texp = index.exptime
texp = map.dur
;dateobs = index.date_obs
dateobs = anytim(map.time, /CCSDS)
IF (iw eq 0) THEN BEGIN
dim =size(image)
nx =dim[1]
ny =dim[2]
nxx =(nx/npix)
nyy =(ny/npix)
i3 =nxx*npix-1
j3 =nyy*npix-1
x =i1+(npix*findgen(nxx)+0.5)
y =j1+(npix*findgen(nyy)+0.5)
images=fltarr(nxx,nyy,nwave)
ENDIF
if (npix eq 1) then images[*,*,iw]=float(image)/texp
if (npix gt 1) then images[*,*,iw]=rebin(float(image[0:i3,0:j3]),nxx,nyy)/texp
texp_(iw)=texp
;PLOT the image as a MAP
id = map.id
;id = strmid(index.TELESCOP,0,3) + '/' + strmid(index.instrume,0,3)
xy_cen = get_map_center(map)
xcen = xy_cen[0]
ycen = xy_cen[1]
aia_lct, rr, gg, bb, wavelnth=wave_[iw], /load
;a = index.CROTA2
;xcen = index.CRVAL1 + index.CDELT1*cos(a)*(0.5*(index.NAXIS1+1)- index.CRPIX1)-index.CDELT2*sin(a)*((index.NAXIS2+1)/2-index.CRPIX2)
;ycen = index.CRVAL2 + index.CDELT1*sin(a)*((index.NAXIS1+1)*0.5-index.CRPIX1) + index.CDELT2*cos(a)*((index.NAXIS2+1)*0.5-index.CRPIX2)
dx = map.dx*npix
dy = map.dy*npix
map = make_map(images[*,*,iw], time = dateobs, id = id, dur = texp, xc = xcen, yc = ycen, dx = dx, dy = dy)
aia_lct, rr, gg, bb, wavelnth=wave_[iw], /load
plot_map, map, /limb, /log
tvscl, images[*,*,iw]
endfor
;________________________TEMPERATURE MAP_________________________________
restore,save_dir + teem_table ;-->wave_,q94,area,resp_corr,telog,dte,tsig,flux
dim =size(flux)
nte =dim[1]
nsig =dim[2]
nwave =dim[3]
ntot =nte*nsig
te_map =fltarr(nxx,nyy)
em_map =fltarr(nxx,nyy)
sig_map =fltarr(nxx,nyy)
flux_dem_map = fltarr(nxx, nyy, nwave)
chi_map = fltarr(nxx,nyy)
em_errmap = fltarr(nxx,nyy)
telog_errmap=fltarr(nxx,nyy,2)
tsig_errmap=fltarr(nxx,nyy,2)
telog_errmap_symmetric=fltarr(nxx,nyy)
tsig_errmap_symmetric=fltarr(nxx,nyy)
chi_map6=fltarr(nxx,nyy,6)
te_best =0.
em_best =0.
sig_best=0.
r0 =0.95*(4096/2)
x0 =4096/2
y0 =4096/2
nfree =3
FOR j = 0, nyy-1 DO BEGIN
if (nx eq 4096) then ind=where(sqrt((x-x0)^2+(y(j)-y0)^2) le r0,nind)
if (nx lt 4096) then ind=findgen(nxx)
i1 = min(ind) > 0
i2 = max(ind) < (nxx-1)
FOR i = i1, i2 DO BEGIN
flux_obs=reform(images[i,j,*])
counts=flux_obs*texp_
noise=sqrt(counts)/texp_
chimin=9999.
chi6min=9999.
chi2d = fltarr(nte,nsig)
FOR k=0,nte-1 DO BEGIN
FOR l=0,nsig-1 DO BEGIN
flux_dem1=reform(flux(k,l,*))
em1 =total(flux_obs)/total(flux_dem1)
em1_err = sqrt(total(noise^2))/total(flux_dem1)
flux_dem=flux_dem1*em1
chi =sqrt(total((flux_obs-flux_dem)^2/noise^2)/(nwave-nfree))
chi6=abs(flux_obs-flux_dem)/noise
chi2d[k,l] = chi
IF (chi le chimin) THEN BEGIN
chimin =chi
chi6min =chi6
em_best =alog10(em1)
em_best_err = em1_err/em_best
te_best =telog(k)
sig_best =tsig(l)
flux_dem_best = flux_dem
ENDIF
ENDFOR
ENDFOR
; find errors in tsig and telog by finding the 1 sigma contour in chi^2 space for each pixel - A. Inglis 5-Oct-2011
IF (chimin ne 9999.) THEN BEGIN
;chimin defaults to 9999. if bad pixel (i.e. no data), so skip those
;find the 1 sigma contour in chi^2-space to get errors in tsig and telog
contour, chi2d, telog, tsig, levels = [chimin +2.3],path_xy=path_xy,path_info=path_info,/path_data_coords,closed=0
IF exist(path_xy) THEN BEGIN
testxmax=MAX(path_xy(0,*))
testxmin=MIN(path_xy(0,*))
testymax=MAX(path_xy(1,*))
testymin=MIN(path_xy(1,*))
telog_err=[testxmin,testxmax] - te_best
;error bars are asymmetric - symmetrise!
telog_errsymmetric=sqrt(telog_err(0)^2 + telog_err(1)^2)
tsig_err=[sig_best - testymin,testymax - sig_best]
;error bars are asymmetric - symmetrise!
;tsig_errsymmetric=sqrt(tsig_err(0)^2 + tsig_err(1)^2)
tsig_errsymmetric = max(tsig_err)
;print,telog_errsymmetric
;print,tsig_errsymmetric
;print,telog_err
;print,tsig_err
ENDIF ELSE BEGIN
;contingency statements in case bad pixel
telog_err=[0,0]
tsig_err=[0,0]
tsig_errsymmetric=0.
telog_errsymmetric=0.
ENDELSE
ENDIF ELSE BEGIN
;contingency statements in case bad pixel
telog_err=[0,0]
tsig_err=[0,0]
tsig_errsymmetric=0.
telog_errsymmetric=0.
ENDELSE
IF keyword_set(VERBOSE) THEN BEGIN
print, i, j
print, 'Using wavelengths:', wave_
print, 'Obs flux:', flux_obs
print, 'Fit flux:', flux_dem_best
print, 'Fit/Obs', flux_dem_best/flux_obs
print, 'min Chisq = ', chimin
print, 'log(T) = ', te_best
print, 'telog err = ', telog_err + te_best
print, 'T sigma = ', sig_best
print, 'T sig err = ', tsig_err + sig_best
print, 'log(EM) = ', em_best
print, ''
loadct, 0
hsi_linecolors
nlevels = 20
levels = chimin * (findgen(nlevels)*0.1 + 1.0)
anot = strarr(20)
FOR y = 0, nlevels-1 DO anot[y] = num2str(levels[y])
contour, chi2d, telog, tsig, levels = levels, c_annotation = anot, xtitle = 'log[Temp]', ytitle = 't_sig'
oplot, [te_best], [sig_best], psym = symcat(16), color = 6
leg = num2str(chimin) + '[log(T) = ' + num2str(te_best) + ', T_sig = ' + num2str(sig_best) + ']'
legend, leg, psym = 4
contour,chi2d,telog,tsig,levels=[chimin + 2.3],/over,thick=2, color = 6
IF keyword_set(DEBUG) THEN stop
ENDIF
flux_dem_map[i,j,*] = flux_dem_best
em_map[i,j]=em_best
te_map[i,j]=te_best
sig_map[i,j]=sig_best
chi_map[i,j]=chimin
chi_map6[i,j,*]=chi6min
telog_errmap[i,j,0:1]=telog_err[0:1]
tsig_errmap[i,j,0:1]=tsig_err[0:1]
telog_errmap_symmetric[i,j]=telog_errsymmetric
tsig_errmap_symmetric[i,j]=tsig_errsymmetric
;telog_errmap(i,j,*)=telog_err
;tsig_errmap(i,j,*)=tsig_err
ENDFOR
IF (j mod 10) eq 0 THEN print,j,nyy
ENDFOR
help,telog_errmap_symmetric
help,tsig_errmap_symmetric
;print,telog_errmap_symmetric(50:70,50:70)
;print,tsig_errmap_symmetric(50:70,50:70)
aia_map = make_map(images[*,*,0], time = dateobs, id = id, dur = texp, xc = xcen, yc = ycen, dx = dx, dy = dy)
aia_map_cube = replicate(aia_map, 6)
aia_simul_map_cube = aia_map_cube
FOR i = 0, n_elements(wave_)-1 DO BEGIN
aia_map_cube[i].data = images[*,*,i]
aia_map_cube[i].id = 'SDO/AIA ' + num2str(wave_[i]) + ' A'
ENDFOR
FOR i = 0, n_elements(wave_)-1 DO BEGIN
aia_simul_map_cube[i].data = flux_dem_map[*,*,i]
aia_simul_map_cube[i].id = 'Simulated SDO/AIA ' + num2str(wave_[i]) + ' A'
ENDFOR
temperature_map = aia_map
temperature_map.data = te_map
temperature_map.id = 'log(Temperature [MK])'
emission_map = aia_map
emission_map.data = em_map
emission_map.id = 'log(Emission Measure [cm!U-3!N]'
sigma_map = aia_map
sigma_map.data = sig_map
sigma_map.id = 'Sigma'
chisq_map = aia_map
chisq_map.data = chi_map
chisq_map.id = textoidl('\chi^2')
tsig_error_map_minus = aia_map
tsig_error_map_minus.data = tsig_errmap[*,*,0]
tsig_error_map_minus.id = 'Tsig error minus'
telog_error_map_minus = aia_map
telog_error_map_minus.data = tsig_errmap[*,*,1]
telog_error_map_minus.id = 'Telog error minus'
tsig_error_map_plus = aia_map
tsig_error_map_plus.data = tsig_errmap[*,*,0]
tsig_error_map_plus.id = 'Tsig error plus'
telog_error_map_plus = aia_map
telog_error_map_plus.data = tsig_errmap[*,*,1]
telog_error_map_plus.id = 'Telog error plus'
;temperature_map = make_map( te_map, xc = xcen, yc = ycen, dx = index.cdelt1*npix, dy = index.cdelt2*npix, id = 'log(Temperature [MK])', time = dateobs)
;emission_map = make_map( em_map, xc = xcen, yc = ycen, dx = index.cdelt1*npix, dy = index.cdelt2*npix, id = 'log(Emission Measure [cm!U-3!N]', time = dateobs)
;sigma_map = make_map( sig_map, xc = xcen, yc = ycen, dx = index.cdelt1*npix, dy = index.cdelt2*npix, id = 'Sigma', time = dateobs)
;chisq_map = make_map( chi_map, xc = xcen, yc = ycen, dx = index.cdelt1*npix, dy = index.cdelt2*npix, id = textoidl('\chi^2'), time = dateobs)
;tsig_error_map_minus = make_map( tsig_errmap(*,*,0), xc = xcen, yc = ycen, dx = index.cdelt1*npix, dy = index.cdelt2*npix, id = textoidl('\chi^2'), time = dateobs)
;telog_error_map_minus = make_map( telog_errmap(*,*,0), xc = xcen, yc = ycen, dx = index.cdelt1*npix, dy = index.cdelt2*npix, id = textoidl('\chi^2'), time = dateobs)
;tsig_error_map_plus = make_map( tsig_errmap(*,*,1), xc = xcen, yc = ycen, dx = index.cdelt1*npix, dy = index.cdelt2*npix, id = textoidl('\chi^2'), time = dateobs)
;telog_error_map_plus = make_map( telog_errmap(*,*,1), xc = xcen, yc = ycen, dx = index.cdelt1*npix, dy = index.cdelt2*npix, id = textoidl('\chi^2'), time = dateobs)
;________________________STATISTICS OF CHI-2 FITS_________________
print,'Statistics of chi2='
statistic,chi_map
print,'Statistics of chi2 (94 A)'
statistic,chi_map6(*,*,5)
print,'log(EM)-range = ',minmax(em_map)
;________________________SAVE MAPS________________________________
;save,filename=teem_map,te_map,em_map,sig_map,chi_map,dateobs
;save,filename=save_dir + teem_map,te_map,em_map,sig_map,chi_map,tsig_errmap_symmetric,telog_errmap_symmetric,temperature_map,emission_map,sigma_map,chisq_map,aia_map_cube, dateobs , aia_simul_map_cube
save,temperature_map,emission_map,sigma_map,chisq_map,tsig_error_map_minus,tsig_error_map_plus,telog_error_map_minus,telog_error_map_plus,aia_simul_map_cube, aia_map_cube, filename=save_dir + 'tempmap' + filename_extra + '.sav'
print,'TE+EM maps saved in file : ',teem_map
t2 =systime(0,/seconds)
print,'Computation time = ',(t2-t1)/60.0,' min'
IF keyword_set(DEBUG) THEN stop
END