OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_sed_extractor.pro 10.9 KB
Edit Raw Normal View History
5c381fd9   Jean-Philippe Bernard   first commit 
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
 
FUNCTION dustem_sed_extractor,maps,index_roi,filters, $
					   		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, $
					   		_extra=_extra

;+
; NAME:
;    sed_extractor
; PURPOSE:
;    Extract SEDs from a set of maps 
; CATEGORY:
;    Dustemwrap
; CALLING SEQUENCE:
;	sed=sed_extractor(maps,index,filters[,comments=][,maps_order=][,method=][,mask_I_undef=][,/total_intensity_only])
; INPUTS:
;    maps 	   : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets]
;    index_roi : spatial index in Nx,Ny where the SED is to be extracted from (ROI=region of interest)
;    filters   : dustemwrap filters corresponding to each maps [Nmaps]
; 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')
; OUTPUTS:
;    SED  : a dustemwrap SED structure containing te extracted SED
; OPTIONAL OUTPUT PARAMETERS:
;    mask_I_undef = a map of saturated pixels found in Stokes I
;    comments   = a vector of strings to be writtena s comments into the SED.
; 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
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    None
; PROCEDURE:
;    None
; EXAMPLES
; MODIFICATION HISTORY:
;    Written by JPB 
;-

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

;===== do a map of saturated stokesI pixels
mask_I_undef=maps[*,*,0,0]
mask_I_undef[*]=0
FOR j=0L,Nfilt-1 DO BEGIN
	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
ENDFOR

used_maps=maps
;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

CASE use_method OF
	'MEAN':BEGIN
		FOR j=0L,Nfilt-1 DO BEGIN
			this_map=used_maps[*,*,indI,j]
			sed[j].stokesI=la_mean(this_map[index_roi])
			this_map=used_maps[*,*,indII,j]
			sed[j].sigmaII=la_div(la_mean(this_map[index_roi]),count_roi)
            IF not keyword_set(total_intensity_only) THEN BEGIN
				this_map=used_maps[*,*,indQ,j]
				sed[j].stokesQ=la_mean(this_map[index_roi])
				this_map=used_maps[*,*,indU,j]
				sed[j].stokesU=la_mean(this_map[index_roi])
				this_map=used_maps[*,*,indQQ,j]
				sed[j].sigmaQQ=la_div(la_mean(this_map[index_roi]),count_roi)
				this_map=used_maps[*,*,indUU,j]
				sed[j].sigmaUU=la_div(la_mean(this_map[index_roi]),count_roi)

				this_map=used_maps[*,*,indIQ,j]
				sed[j].sigmaIQ=la_div(la_mean(this_map[index_roi]),count_roi)
				this_map=used_maps[*,*,indIU,j]
				sed[j].sigmaIU=la_div(la_mean(this_map[index_roi]),count_roi)
				this_map=used_maps[*,*,indQU,j]
				sed[j].sigmaQU=la_div(la_mean(this_map[index_roi]),count_roi)
			ENDIF
		ENDFOR
	END
	'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
			this_map=used_maps[*,*,indI,j]
			sed[j].stokesI=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
			this_map=used_maps[*,*,indII,j]
			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
				this_map=used_maps[*,*,indQ,j]
				sed[j].stokesQ=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
				this_map=used_maps[*,*,indU,j]
				sed[j].stokesU=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
	            ;variances
				this_map=used_maps[*,*,indQQ,j]
				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)
				this_map=used_maps[*,*,indUU,j]
				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
				this_map=used_maps[*,*,indIQ,j]
				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)
				this_map=used_maps[*,*,indIU,j]
				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)
				this_map=used_maps[*,*,indQU,j]
				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
	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
        mapI_ref=used_maps[*,*,indI,indf]
        mapQ_ref=used_maps[*,*,indQ,indf]
        mapU_ref=used_maps[*,*,indU,indf]
        mapII_ref=used_maps[*,*,indII,indf]
        mapQQ_ref=used_maps[*,*,indQQ,indf]
        mapUU_ref=used_maps[*,*,indUU,indf]
		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]
			y=(used_maps[*,*,indI,j])[index_roi]
			sx=la_power(mapII_ref[index_roi],0.5)
			sy=la_power((used_maps[*,*,indI,j])[index_roi],0.5)
			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]]
			y=[(used_maps[*,*,indQ,j])[index_roi],(used_maps[*,*,indU,j])[index_roi]]
			sx=la_power([mapQQ_ref[index_roi],mapUU_ref[index_roi]],0.5)
			sy=la_power([(used_maps[*,*,indQ,j])[index_roi],(used_maps[*,*,indU,j])[index_roi]],0.5)
			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

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

;stop

the_end:
RETURN,sed

END