Blame view

src/idl/dustem_sed_extractor.pro 12.3 KB
88f84272   Jean-Philippe Bernard   changed
1
2
3
FUNCTION dustem_sed_extractor,maps, $
                              index_roi, $
                              filters, $
ecfb8fb3   Annie Hughes   updated comments ...
4
5
6
7
8
9
10
                              comments=comments, $
                              maps_order=maps_order, $
                              method=method, $
                              index_off=index_off, $
                              reference_filter=reference_filter, $
                              mask_I_undef=mask_I_undef, $
                              total_intensity_only=total_intensity_only, $
057f2b9c   Jean-Philippe Bernard   small changes
11
                              set_undefined=set_undefined, $
ecfb8fb3   Annie Hughes   updated comments ...
12
                              _extra=_extra
5c381fd9   Jean-Philippe Bernard   first commit
13
14
15

;+
; NAME:
cd49cdc0   Annie Hughes   updated help
16
;    dustem_sed_extractor
ecfb8fb3   Annie Hughes   updated comments ...
17
;
5c381fd9   Jean-Philippe Bernard   first commit
18
19
; PURPOSE:
;    Extract SEDs from a set of maps 
ecfb8fb3   Annie Hughes   updated comments ...
20
;
5c381fd9   Jean-Philippe Bernard   first commit
21
; CATEGORY:
ecfb8fb3   Annie Hughes   updated comments ...
22
23
;    DustEMWrap, Distributed, High-Level, Analysis Helper
;
5c381fd9   Jean-Philippe Bernard   first commit
24
25
; CALLING SEQUENCE:
;	sed=sed_extractor(maps,index,filters[,comments=][,maps_order=][,method=][,mask_I_undef=][,/total_intensity_only])
ecfb8fb3   Annie Hughes   updated comments ...
26
;
5c381fd9   Jean-Philippe Bernard   first commit
27
; INPUTS:
88f84272   Jean-Philippe Bernard   changed
28
;    maps      : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets]. Caution, maps are changed if set_undefined is set !
5c381fd9   Jean-Philippe Bernard   first commit
29
;    index_roi : spatial index in Nx,Ny where the SED is to be extracted from (ROI=region of interest)
ecfb8fb3   Annie Hughes   updated comments ...
30
31
;    filters   : DustEMWrap filters corresponding to each dataset [Ndata_sets]
;
5c381fd9   Jean-Philippe Bernard   first commit
32
33
34
35
36
37
38
; OPTIONAL INPUT PARAMETERS:
;    comments   = a vector of strings to be attached to the SED. Also an output.
;    maps_order = a vector of 9 integers giving the 0-based index for I,Q,U,II,QQ,UU,IQ,IU,QU in that order
;                 (default [0,1,2,3,4,5,6,7,8]). For total intensity only, use something like [0,-1,-1,1,-1,-1,-1,-1,-1].
;    method = extraction method (default = 'MEAN')/ Possible values : 'MEAN','MEAN_ONOFF','QU_CORREL'
;    index_off = same as index for the OFF region (needed only for method='MEAN_ONOFF')
;    reference_filter = reference filter to be correlated with each map (needed only for method='QU_CORREL')
ecfb8fb3   Annie Hughes   updated comments ...
39
;
5c381fd9   Jean-Philippe Bernard   first commit
40
41
; OUTPUTS:
;    SED  : a dustemwrap SED structure containing te extracted SED
ecfb8fb3   Annie Hughes   updated comments ...
42
;
5c381fd9   Jean-Philippe Bernard   first commit
43
44
; OPTIONAL OUTPUT PARAMETERS:
;    mask_I_undef = a map of saturated pixels found in Stokes I
ecfb8fb3   Annie Hughes   updated comments ...
45
46
;    comments   = a vector of strings to be written as comments into the SED.
;
5c381fd9   Jean-Philippe Bernard   first commit
47
48
49
; ACCEPTED KEY-WORDS:
;    total_intensity_only   = if set, data is assumed to contain only total intensity (no polarization)
;    comments               = comments to be added to the sed. sed_extractor adds comments to this
ecfb8fb3   Annie Hughes   updated comments ...
50
;
5c381fd9   Jean-Philippe Bernard   first commit
51
52
; COMMON BLOCKS:
;    None
ecfb8fb3   Annie Hughes   updated comments ...
53
;
5c381fd9   Jean-Philippe Bernard   first commit
54
55
; SIDE EFFECTS:
;    None
ecfb8fb3   Annie Hughes   updated comments ...
56
;
5c381fd9   Jean-Philippe Bernard   first commit
57
58
; RESTRICTIONS:
;    None
ecfb8fb3   Annie Hughes   updated comments ...
59
;
5c381fd9   Jean-Philippe Bernard   first commit
60
61
; PROCEDURE:
;    None
ecfb8fb3   Annie Hughes   updated comments ...
62
;
5c381fd9   Jean-Philippe Bernard   first commit
63
; EXAMPLES
ecfb8fb3   Annie Hughes   updated comments ...
64
;
5c381fd9   Jean-Philippe Bernard   first commit
65
; MODIFICATION HISTORY:
ecfb8fb3   Annie Hughes   updated comments ...
66
67
68
;    Written by JPB 2022
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
5c381fd9   Jean-Philippe Bernard   first commit
69
70
;-

ecfb8fb3   Annie Hughes   updated comments ...
71
72


5c381fd9   Jean-Philippe Bernard   first commit
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
IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_sed_extractor'
  sed=0.
  goto,the_end
ENDIF

Nfilt=n_elements(filters)
sed=dustem_initialize_sed(Nfilt,comments=comments)
sed.filter=filters
sed.instru=dustem_filter2instru(sed.filter)
sed.wave=dustem_filter2wav(filters)
comments=[comments,'\Computed using sed_extractor on '+systime(0)]
order=sort(sed.wave)  ;wavelength order

count_roi=n_elements(index_roi)

use_method='MEAN'  ;This is the default extraction method
IF keyword_set(method) THEN BEGIN
	use_method=method
ENDIF
comments=[comments,'SED extracted using method '+use_method]
comments=[comments,'ROI has '+strtrim(count_roi,2)+' map pixels']
N1=(size(maps))[1]
N2=(size(maps))[2]
Nplanes=(size(maps))[3]

IF keyword_set(total_intensity_only) THEN BEGIN
	indI=0L & indII=1L
	IF keyword_set(maps_order) THEN BEGIN
		indI=maps_order[0]
		indII=maps_order[1]
	ENDIF
ENDIF ELSE BEGIN
	indI=0L & indQ=1L & indU=2L & indII=3L & indQQ=4L & indUU=5L &indIQ=6L & indIU=7L &indQU=8L
	IF keyword_set(maps_order) THEN BEGIN
		indI=maps_order[0]
		indQ=maps_order[1]
		indU=maps_order[2]
		indII=maps_order[3]
		indQQ=maps_order[4]
		indUU=maps_order[5]
		indIQ=maps_order[6]
		indIU=maps_order[7]
		indQU=maps_order[8]
	ENDIF
ENDELSE

88f84272   Jean-Philippe Bernard   changed
120
;used_maps=maps
057f2b9c   Jean-Philippe Bernard   small changes
121
122

IF keyword_set(set_undefined) THEN BEGIN
88f84272   Jean-Philippe Bernard   changed
123
	used_maps=maps
057f2b9c   Jean-Philippe Bernard   small changes
124
125
126
	;===== do a map of saturated stokesI pixels
	mask_I_undef=maps[*,*,0,0]
	mask_I_undef[*]=0
5c381fd9   Jean-Philippe Bernard   first commit
127
	FOR j=0L,Nfilt-1 DO BEGIN
057f2b9c   Jean-Philippe Bernard   small changes
128
129
130
131
132
133
134
		this_map=maps[*,*,indI,j]
		ind=where(this_map EQ la_undef(),count)
		IF count NE 0 AND count NE N1*N2 THEN BEGIN
			message,'Found '+strtrim(count,2)+' undefined values for map of '+filters[j],/continue
			;stop
			mask_I_undef[ind]=1
		ENDIF
5c381fd9   Jean-Philippe Bernard   first commit
135
	ENDFOR
057f2b9c   Jean-Philippe Bernard   small changes
136
137
138
139
140
141
142
143
144
145
146
147
	;set saturated stokesI pixels to la_undef in 
	indbad=where(mask_I_undef EQ 1,count)
	IF count NE 0 THEN BEGIN
		message,'Setting Stokes I saturated values to undef in the whole cube',/continue
		FOR j=0L,Nfilt-1 DO BEGIN
			FOR k=0L,Nplanes-1 DO BEGIN
				this_map=used_maps[*,*,k,j]
				this_map[indbad]=la_undef()
		  		used_maps[*,*,k,j]=this_map
		 	ENDFOR
		ENDFOR
	ENDIF
88f84272   Jean-Philippe Bernard   changed
148
149
	maps=used_maps
	used_maps=0
5c381fd9   Jean-Philippe Bernard   first commit
150
151
152
153
ENDIF

CASE use_method OF
	'MEAN':BEGIN
88f84272   Jean-Philippe Bernard   changed
154
155
		   FOR j=0L,Nfilt-1 DO BEGIN
                   this_map=maps[*,*,indI,j]
ecfb8fb3   Annie Hughes   updated comments ...
156
                   sed[j].stokesI=la_mean(this_map[index_roi])
88f84272   Jean-Philippe Bernard   changed
157
                   this_map=maps[*,*,indII,j]
ecfb8fb3   Annie Hughes   updated comments ...
158
159
                   sed[j].sigmaII=la_div(la_mean(this_map[index_roi]),count_roi)
                   IF not keyword_set(total_intensity_only) THEN BEGIN
88f84272   Jean-Philippe Bernard   changed
160
                      this_map=maps[*,*,indQ,j]
ecfb8fb3   Annie Hughes   updated comments ...
161
                      sed[j].stokesQ=la_mean(this_map[index_roi])
88f84272   Jean-Philippe Bernard   changed
162
                      this_map=maps[*,*,indU,j]
ecfb8fb3   Annie Hughes   updated comments ...
163
                      sed[j].stokesU=la_mean(this_map[index_roi])
88f84272   Jean-Philippe Bernard   changed
164
                      this_map=maps[*,*,indQQ,j]
ecfb8fb3   Annie Hughes   updated comments ...
165
                      sed[j].sigmaQQ=la_div(la_mean(this_map[index_roi]),count_roi)
88f84272   Jean-Philippe Bernard   changed
166
                      this_map=maps[*,*,indUU,j]
ecfb8fb3   Annie Hughes   updated comments ...
167
168
                      sed[j].sigmaUU=la_div(la_mean(this_map[index_roi]),count_roi)
                      
88f84272   Jean-Philippe Bernard   changed
169
                      this_map=maps[*,*,indIQ,j]
ecfb8fb3   Annie Hughes   updated comments ...
170
                      sed[j].sigmaIQ=la_div(la_mean(this_map[index_roi]),count_roi)
88f84272   Jean-Philippe Bernard   changed
171
                      this_map=maps[*,*,indIU,j]
ecfb8fb3   Annie Hughes   updated comments ...
172
                      sed[j].sigmaIU=la_div(la_mean(this_map[index_roi]),count_roi)
88f84272   Jean-Philippe Bernard   changed
173
                      this_map=maps[*,*,indQU,j]
ecfb8fb3   Annie Hughes   updated comments ...
174
175
                      sed[j].sigmaQU=la_div(la_mean(this_map[index_roi]),count_roi)
                   ENDIF
88f84272   Jean-Philippe Bernard   changed
176
177
            ENDFOR
         END
ecfb8fb3   Annie Hughes   updated comments ...
178
179
180
181
182
183
184
185
        'MEAN_ONOFF': BEGIN     ;This is mean ON - mean MEAN_ONOFF
           IF not keyword_set(index_off) THEN BEGIN
              message,'index_off keyword must be set with method '+use_method,/continue
              stop
           ENDIF
           count_off=n_elements(index_off)
           comments=[comments,'OFF has '+strtrim(count_off,2)+' map pixels']
           FOR j=0L,Nfilt-1 DO BEGIN
88f84272   Jean-Philippe Bernard   changed
186
              this_map=maps[*,*,indI,j]
ecfb8fb3   Annie Hughes   updated comments ...
187
              sed[j].stokesI=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
88f84272   Jean-Philippe Bernard   changed
188
              this_map=maps[*,*,indII,j]
ecfb8fb3   Annie Hughes   updated comments ...
189
190
191
192
              sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
              sig_off=la_div(la_mean(this_map[index_off]),count_off)
              sed[j].sigmaII=la_add(sig_on,sig_off)
              IF not keyword_set(total_intensity_only) THEN BEGIN
88f84272   Jean-Philippe Bernard   changed
193
                 this_map=maps[*,*,indQ,j]
ecfb8fb3   Annie Hughes   updated comments ...
194
                 sed[j].stokesQ=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
88f84272   Jean-Philippe Bernard   changed
195
                 this_map=maps[*,*,indU,j]
ecfb8fb3   Annie Hughes   updated comments ...
196
197
                 sed[j].stokesU=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
                                ;variances
88f84272   Jean-Philippe Bernard   changed
198
                 this_map=maps[*,*,indQQ,j]
ecfb8fb3   Annie Hughes   updated comments ...
199
200
201
                 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
                 sig_off=la_div(la_mean(this_map[index_off]),count_off)
                 sed[j].sigmaQQ=la_add(sig_on,sig_off)
88f84272   Jean-Philippe Bernard   changed
202
                 this_map=maps[*,*,indUU,j]
ecfb8fb3   Annie Hughes   updated comments ...
203
204
205
206
                 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
                 sig_off=la_div(la_mean(this_map[index_off]),count_off)			
                 sed[j].sigmaUU=la_add(sig_on,sig_off)
                                ;co-variances
88f84272   Jean-Philippe Bernard   changed
207
                 this_map=maps[*,*,indIQ,j]
ecfb8fb3   Annie Hughes   updated comments ...
208
209
210
                 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
                 sig_off=la_div(la_mean(this_map[index_off]),count_off)			
                 sed[j].sigmaIQ=la_add(sig_on,sig_off)
88f84272   Jean-Philippe Bernard   changed
211
                 this_map=maps[*,*,indIU,j]
ecfb8fb3   Annie Hughes   updated comments ...
212
213
214
                 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
                 sig_off=la_div(la_mean(this_map[index_off]),count_off)			
                 sed[j].sigmaIU=la_add(sig_on,sig_off)
88f84272   Jean-Philippe Bernard   changed
215
                 this_map=maps[*,*,indQU,j]
ecfb8fb3   Annie Hughes   updated comments ...
216
217
218
219
220
221
                 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
                 sig_off=la_div(la_mean(this_map[index_off]),count_off)			
                 sed[j].sigmaQU=la_add(sig_on,sig_off)
              ENDIF
           ENDFOR	
                                ;stop
5c381fd9   Jean-Philippe Bernard   first commit
222
223
224
225
226
227
228
229
230
231
232
233
234
235
	END
	'QU_CORREL': BEGIN    ;This is to compute SED of polarization correlated to a reference
        IF keyword_set(total_intensity_only) THEN BEGIN
        	message,'mode QU_CORREL can only be used with polarization data'
        ENDIF
	    IF not keyword_set(reference_filter) THEN BEGIN
	    	message,'reference_filter keyword must be set with method '+use_method,/continue
	    	stop
	    ENDIF
	    indf=where(filters EQ reference_filter,countf)
	    IF countf EQ 0 THEN BEGIN
	    	message,'reference filter '+reference_filter+' not found in provided filter list',/continue
	    	stop
	    ENDIF
88f84272   Jean-Philippe Bernard   changed
236
237
238
239
240
241
        mapI_ref=maps[*,*,indI,indf]
        mapQ_ref=maps[*,*,indQ,indf]
        mapU_ref=maps[*,*,indU,indf]
        mapII_ref=maps[*,*,indII,indf]
        mapQQ_ref=maps[*,*,indQQ,indf]
        mapUU_ref=maps[*,*,indUU,indf]
5c381fd9   Jean-Philippe Bernard   first commit
242
243
244
245
246
247
		Iref=la_mean(mapI_ref[index_roi])
		Qref=la_mean(mapQ_ref[index_roi])
		Uref=la_mean(mapU_ref[index_roi])
		FOR j=0L,Nfilt-1 DO BEGIN
			;==== do the correlation in Stokes I
			x=mapI_ref[index_roi]
88f84272   Jean-Philippe Bernard   changed
248
			y=(maps[*,*,indI,j])[index_roi]
5c381fd9   Jean-Philippe Bernard   first commit
249
			sx=la_power(mapII_ref[index_roi],0.5)
88f84272   Jean-Philippe Bernard   changed
250
			sy=la_power((maps[*,*,indI,j])[index_roi],0.5)
5c381fd9   Jean-Philippe Bernard   first commit
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
			ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good)
			comments=[comments,'I corelation has '+strtrim(count_good,2)+' map pixels']
			IF count_good NE 0 THEN BEGIN
	            x=x[ind_good]
	            y=y[ind_good]
	            sx=sx[ind_good]
	            sy=sy[ind_good]
	            start=[1.,1.]
		    	try=linear_mpfit(x,y,sx,sy,start,status=status $
	                      ,perror=perror,bestnorm=bestnorm,_extra=_extra)
		    	;stop
	    	ENDIF ELSE BEGIN
	    		status=0
	    		try=[la_undef(),la_undef()]
	    	ENDELSE
	    	IF status NE 0 AND status NE -16 THEN BEGIN
				sed[j].stokesI=try[0]*Iref
				sed[j].sigmaII=la_power(la_mul(perror[0],Iref),2.)
			ENDIF
			;Do the correlation between Q maps
			x=[mapQ_ref[index_roi],mapU_ref[index_roi]]
88f84272   Jean-Philippe Bernard   changed
272
			y=[(maps[*,*,indQ,j])[index_roi],(maps[*,*,indU,j])[index_roi]]
5c381fd9   Jean-Philippe Bernard   first commit
273
			sx=la_power([mapQQ_ref[index_roi],mapUU_ref[index_roi]],0.5)
88f84272   Jean-Philippe Bernard   changed
274
			sy=la_power([(maps[*,*,indQ,j])[index_roi],(maps[*,*,indU,j])[index_roi]],0.5)
5c381fd9   Jean-Philippe Bernard   first commit
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
			ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good)
			comments=[comments,'QU corelation has '+strtrim(count_good,2)+' map pixels']
			IF count_good NE 0 THEN BEGIN
	            x=x[ind_good]
	            y=y[ind_good]
	            sx=sx[ind_good]
	            sy=sy[ind_good]
	            start=[1.,1.]
		    	try=linear_mpfit(x,y,sx,sy,start,status=status $
                      ,perror=perror,bestnorm=bestnorm,_extra=_extra)
		    	;stop
	    	ENDIF ELSE BEGIN
	    		status=0
	    		try=[la_undef(),la_undef()]
	    	ENDELSE
	    	IF status NE 0 AND status NE -16 THEN BEGIN
				sed[j].stokesQ=la_mul(try[0],Qref)
				sed[j].stokesU=la_mul(try[0],Uref)
				sed[j].sigmaQQ=la_power(la_mul(perror[0],Qref),2.)
				sed[j].sigmaUU=la_power(la_mul(perror[0],Uref),2.)
	            ;for now assume no co-variances. But we could aclculate this, right ?
				sed[j].sigmaIQ=0.
				sed[j].sigmaIU=0.
				sed[j].sigmaQU=0.
				;IF filters[j] EQ reference_filter THEN stop
			ENDIF
		ENDFOR
		comments=[comments,'Reference filter for QU correlation '+reference_filter]     
	END
    ELSE: BEGIN
    	message,'Method '+use_method+' not implemented',/continue
    	sed=0L
    	goto,the_end
    END
ENDCASE

057f2b9c   Jean-Philippe Bernard   small changes
311
312
313
314
315
316
317
318
319
320
321
322
IF not keyword_set(total_intensity_only) THEN BEGIN
	polar_iqu2ippsi,sed.stokesI,sed.stokesQ,sed.stokesU,p_values,psi_values,/lac,largep=largep
	sed.smallp=p_values
	sed.psi=psi_values
	sed.largep=largep
	polar_variance_iqu2ippsi,sed.stokesI,sed.stokesQ,sed.stokesU,sed.sigmaII,sed.sigmaQQ,sed.sigmaUU,sed.sigmaIQ,sed.sigmaIU,sed.sigmaQU,variances_p,variances_psi,variances_largep,/lac
	sed.sigma_smallp=variances_p
	sed.sigma_largep=variances_largep
	sed.sigma_psi=variances_psi
	ssed=sed[order]
	sed=ssed
ENDIF
5c381fd9   Jean-Philippe Bernard   first commit
323
324
325
326
327
328

;stop

the_end:
RETURN,sed

ecfb8fb3   Annie Hughes   updated comments ...
329
END