dustem_check_data.pro
18.7 KB
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
FUNCTION dustem_check_data,data_sed=data_sed,data_ext=data_ext,sed=sed,ext=ext,polext=polext,polsed=polsed,polfrac=polfrac,psi_em=psi_em,qsed=qsed,used=used,qext=qext,uext=uext,psi_ext=psi_ext
;WRITE HELP SECTION . URGENT
;NOTA BENE: LargeP has not been debiasied yet.
;CHECKING THE FORMAT OF THE SED STRUCTURE
IF KEYWORD_SET(data_sed) THEN BEGIN
ans=''
wavs=data_sed.wave
;=== Impose central wavelengths for photometric channels
ind=where(data_sed.filter NE 'SPECTRUM',count)
IF count NE 0 THEN BEGIN
wavs(ind)=dustem_filter2wav(data_sed(ind).filter)
ENDIF
;=== Initializing observation structure
obs_EM={instru_names:data_sed.instru,filt_names:data_sed.filter,wav:wavs, values:data_sed.StokesI,sigma:sqrt(abs(data_sed.sigmaII))}
;INITIALIZING OUTPUT STRUCTURES
sed = obs_EM
qsed = sed
used = sed
;Test if there is StokesI data. If not move to polsed testing.
ind = where(data_sed.StokesI NE la_undef() and finite(data_sed.StokesI) EQ 1 and data_sed.sigmaII NE la_undef() and finite(data_sed.sigmaII) EQ 1,N_sed) ;Testing if the data and errors are set
IF N_sed EQ 0 THEN BEGIN
sed = ptr_new()
goto, qzone
ENDIF ELSE BEGIN
ind = where(finite(data_sed.StokesI) EQ 0 or finite(data_sed.sigmaII) EQ 0, countnoki) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
IF countnoki NE 0 THEN BEGIN
message, "SED format isn't met. NaN value(s) detected in StokesI data.",/continue
message, 'Without automatic format fitering you cannot proceed.',/continue
read, ans, prompt='Would you like automatic format filtering? (Y/N):'
IF strupcase(ans) EQ 'N' THEN BEGIN
message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
message, 'Aborting ...'
ENDIF
ENDIF
IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
ELSE message, 'Data format is met. Removal of undefined data points...',/continue
tdaerri=where((data_sed.sigmaII EQ la_undef() and data_sed.StokesI NE la_undef()) or (data_sed.sigmaII NE la_undef() and data_sed.StokesI EQ la_undef()),ctdaerri)
IF ctdaerri NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
testi=where(data_sed.sigmaII NE 0. and finite(data_sed.sigmaII) EQ 1 and data_sed.sigmaII NE la_undef() and data_sed.StokesI NE la_undef() and finite(data_sed.StokesI) EQ 1,ctesti)
IF ctesti NE 0 THEN BEGIN
instru_names = sed.instru_names(testi)
filt_names = sed.filt_names(testi)
wav = sed.wav(testi)
values = sed.values(testi)
sigma = sed.sigma(testi)
sed={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
ENDIF
ENDELSE
qzone:
ans=''
qsed.values=data_sed.StokesQ
qsed.sigma=sqrt(abs(data_sed.sigmaQQ))
;Testig if there is StokesQ data. If not move to used data
ind = where(data_sed.StokesQ NE la_undef() and finite(data_sed.StokesQ) EQ 1 and data_sed.sigmaQQ NE la_undef() and finite(data_sed.sigmaQQ) EQ 1,NQ)
IF NQ EQ 0 THEN BEGIN
qsed = ptr_new()
goto, uzone
ENDIF ELSE BEGIN
ind = where(finite(data_sed.StokesQ) EQ 0 or finite(data_sed.sigmaQQ) EQ 0, countnokq) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
IF countnokq NE 0 THEN BEGIN
message, "SED format isn't met. NaN value(s) detected in StokesQ data.",/continue
message, 'Without automatic format fitering you cannot proceed.',/continue
read, ans, prompt='Would you like automatic format filtering? (Y/N):'
IF strupcase(ans) EQ 'N' THEN BEGIN
message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
message, 'Aborting ...'
ENDIF
ENDIF
IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
ELSE message, 'Data format is met. Removal of undefined data points...',/continue
tdaerrq=where((data_sed.sigmaQQ EQ la_undef() and data_sed.StokesQ NE la_undef()) or (data_sed.sigmaQQ NE la_undef() and data_sed.StokesQ EQ la_undef()),ctdaerrq)
IF ctdaerrq NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
testq=where(data_sed.sigmaQQ NE 0. and finite(data_sed.sigmaQQ) EQ 1 and data_sed.sigmaQQ NE la_undef() and data_sed.StokesQ NE la_undef() and finite(data_sed.StokesQ) EQ 1,ctestq)
IF ctestq NE 0 THEN BEGIN
instru_names = qsed.instru_names(testq)
filt_names = qsed.filt_names(testq)
wav = qsed.wav(testq)
values = qsed.values(testq)
sigma = qsed.sigma(testq)
qsed={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
ENDIF
ENDELSE
uzone:
ans=''
used.values=data_sed.StokesU
used.sigma=sqrt(abs(data_sed.sigmaUU))
;Testig if there is StokesU data. If not return empty value and end loop.
ind = where(data_sed.StokesU NE la_undef() and finite(data_sed.StokesU) EQ 1 and data_sed.sigmaUU NE la_undef() and finite(data_sed.sigmaUU) EQ 1,NU)
IF NU EQ 0 THEN BEGIN
used = ptr_new()
goto, psiemzone
ENDIF ELSE BEGIN
ind = where(finite(data_sed.StokesU) EQ 0 or finite(data_sed.sigmaUU) EQ 0, countnoku) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
IF countnoku NE 0 THEN BEGIN
message, "SED format isn't met. NaN value(s) detected in StokesU data.",/continue
message, 'Without automatic format fitering you cannot proceed.',/continue
read, ans, prompt='Would you like automatic format filtering? (Y/N):'
IF strupcase(ans) EQ 'N' THEN BEGIN
message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
message, 'Aborting ...'
ENDIF
ENDIF
IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
ELSE message, 'Data format is met. Removal of undefined data points...',/continue
tdaerru=where((data_sed.sigmaUU EQ la_undef() and data_sed.StokesU NE la_undef()) or (data_sed.sigmaUU NE la_undef() and data_sed.StokesU EQ la_undef()),ctdaerru)
IF ctdaerru NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
testu=where(data_sed.sigmaUU NE 0. and finite(data_sed.sigmaUU) EQ 1 and data_sed.sigmaUU NE la_undef() and data_sed.StokesU NE la_undef() and finite(data_sed.StokesU) EQ 1,ctestu)
IF ctestu NE 0 THEN BEGIN
instru_names = used.instru_names(testu)
filt_names = used.filt_names(testu)
wav = used.wav(testu)
values = used.values(testu)
sigma = used.sigma(testu)
used={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
ENDIF
ENDELSE
psiemzone:
if nu eq 0 and nq eq 0 then goto, end_sed
degtorad = !pi/180
;Generation of POLFRAC,POLSED and PSI data by default
;Testing on the presence of covariances otherwise an error is returned
;The testing happens here because covariances are needed for the error estimation of P,p and PSI.
tstcovu = where((data_sed(testu).sigmaIU) EQ la_undef() or (data_sed(testu).sigmaQU) EQ la_undef(), ctcovu)
if ctcovu NE 0 then message, 'Covariance data is missing in StokesU data. Data, errors and covariances do not match. Aborting...'
tstcovq = where((data_sed(testq).sigmaIQ) EQ la_undef() or (data_sed(testq).sigmaQU) EQ la_undef(), ctcovq)
if ctcovq NE 0 then message, 'Covariance data is missing in StokesQ data. Data, errors and covariances do not match. Aborting...'
;Matching both filter and spectrum points
match2,used.wav,qsed.wav,qpts,upts
;removing -1 values from the indices (qpts & upts)
indnotq = where(qpts eq -1, cnotq)
if cnotq ne 0 then remove, indnotq, qpts
indnotu = where(upts eq -1, cnotu)
if cnotu ne 0 then remove, indnotu, upts
match2, sed.wav, qsed.wav(qpts), qtmp, sedpts
indnoti = where(sedpts eq -1, cnoti)
if cnoti ne 0 then remove, indnoti, sedpts
instru_names = used.instru_names(upts)
filt_names = used.filt_names(upts)
wav = used.wav(upts)
values = 0.5*atan(used.values(upts),qsed.values(qpts))/degtorad
polar_variance_iqu2ippsi,sed.values(sedpts),qsed.values(qpts),used.values(upts),(data_sed(testi).sigmaII)(sedpts),(data_sed(testq).sigmaQQ)(qpts),(data_sed(testu).sigmaUU)(upts),(data_sed(testi).sigmaIQ)(sedpts),(data_sed(testi).sigmaIU)(sedpts),(data_sed(testi).sigmaQU)(sedpts),variance_smallp,variance_psi,variance_largep
sigma = sqrt(variance_psi)
psi_em={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
goto, pzone ; the variances for polfrac and polsed were already computed.
pzone:
polsed = psi_em
polsed.values = sqrt(used.values(upts)^2+qsed.values(qpts)^2)
polsed.sigma = sqrt(variance_largep)
goto, fpzone
fpzone:
polfrac = psi_em
polfrac.values = polsed.values/sed.values(sedpts)
polfrac.sigma = sqrt(variance_smallp)
goto, end_sed
end_sed:
ENDIF
;EXTINCTION HASN'T BEEN TREATED YET
IF KEYWORD_SET(data_ext) THEN begin
ans=''
wavs=data_ext.wave
;=== Impose central wavelengths for photometric channels
ind=where(data_ext.filter NE 'SPECTRUM',count)
IF count NE 0 THEN BEGIN
wavs(ind)=dustem_filter2wav(data_ext(ind).filter)
ENDIF
;=== define observations structure
obs_EXT={instru_names:data_ext.instru,filt_names:data_ext.filter,wav:wavs, values:data_ext.EXT_I,sigma:sqrt(abs(data_ext.sigextII))}
;INITIALIZING OUTPUT STRUCTURES
ext = obs_EXT
qext = ext
uext = ext
ind = where(data_ext.EXT_I NE la_undef() and finite(data_ext.EXT_I) EQ 1 and data_ext.sigextII NE la_undef() and finite(data_ext.sigextII) EQ 1,N_ext) ;Testing if the data and errors are set
IF N_ext EQ 0 THEN BEGIN
ext = ptr_new()
goto, qexzone
ENDIF ELSE BEGIN
ind = where(finite(data_ext.EXT_I) EQ 0 or finite(data_ext.sigextII) EQ 0, countnokex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
IF countnokex NE 0 THEN BEGIN
message, "SED format isn't met. NaN value(s) detected in EXT_I data.",/continue
message, 'Without automatic format fitering you cannot proceed.',/continue
ans=''
read, ans, prompt='Would you like automatic format filtering? (Y/N):'
IF strupcase(ans) EQ 'N' THEN BEGIN
message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
message, 'Aborting ...'
ENDIF
ENDIF
IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
ELSE message, 'Data format is met. Removal of undefined data points...',/continue
tdaerrex=where((data_ext.sigextII EQ la_undef() and data_ext.EXT_I NE la_undef()) or (data_ext.sigextII NE la_undef() and data_ext.EXT_I EQ la_undef()),ctdaerrex)
IF ctdaerrex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
testex=where(data_ext.sigextII NE 0. and finite(data_ext.sigextII) EQ 1 and data_ext.sigextII NE la_undef() and data_ext.EXT_I NE la_undef() and finite(data_ext.EXT_I) EQ 1,ctestex)
IF ctestex NE 0 THEN BEGIN
instru_names = ext.instru_names(testex)
filt_names = ext.filt_names(testex)
wav = ext.wav(testex)
values = ext.values(testex)
sigma = ext.sigma(testex)
ext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
ENDIF
ENDELSE
qexzone:
ans=''
qext.values=data_ext.EXT_Q
qext.sigma=sqrt(abs(data_ext.sigextQQ))
;Testig if there is EXT_P data. If not move to EXT_U testing
ind = where(data_ext.EXT_Q NE la_undef() and finite(data_ext.EXT_Q) EQ 1 and data_ext.sigextQQ NE la_undef() and finite(data_ext.sigextQQ) EQ 1,NQ_ext) ;Testing if the data and errors are set
;dustem_set_data.prostop
IF NQ_ext EQ 0 THEN BEGIN
qext = ptr_new()
goto, uexzone
ENDIF ELSE BEGIN
ind = where(finite(data_ext.EXT_Q) EQ 0 or finite(data_ext.sigextQQ) EQ 0, countnokqex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
IF countnokqex NE 0 THEN BEGIN
message, "SED format isn't met. NaN value(s) detected in EXT_Q data.",/continue
message, 'Without automatic format fitering you cannot proceed.',/continue
ans=''
read, ans, prompt='Would you like automatic format filtering? (Y/N):'
IF strupcase(ans) EQ 'N' THEN BEGIN
message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
message, 'Aborting ...'
ENDIF
ENDIF
IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
ELSE message, 'Data format is met. Removal of undefined data points...',/continue
tdaerrqex=where((data_ext.sigextQQ EQ la_undef() and data_ext.EXT_Q NE la_undef()) or (data_ext.sigextQQ NE la_undef() and data_ext.EXT_Q EQ la_undef()),ctdaerrqex)
IF ctdaerrqex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
testqex=where(data_ext.sigextQQ NE 0. and finite(data_ext.sigextQQ) EQ 1 and data_ext.sigextQQ NE la_undef() and data_ext.EXT_Q NE la_undef() and finite(data_ext.EXT_Q) EQ 1,ctestqex)
IF ctestqex NE 0 THEN BEGIN
instru_names = qext.instru_names(testqex)
filt_names = qext.filt_names(testqex)
wav = qext.wav(testqex)
values = qext.values(testqex)
sigma = qext.sigma(testqex)
qext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
ENDIF
ENDELSE
uexzone:
ans=''
uext.values=data_ext.EXT_U
uext.sigma=sqrt(abs(data_ext.sigextUU))
;Testig if there is EXT_P data. If not return empty value and end loop.
ind = where(data_ext.EXT_U NE la_undef() and finite(data_ext.EXT_U) EQ 1 and data_ext.sigextUU NE la_undef() and finite(data_ext.sigextUU) EQ 1,NU_ext) ;Testing if the data and errors are set
IF NU_ext EQ 0 THEN BEGIN
uext = ptr_new()
goto, psiexzone
ENDIF ELSE BEGIN
ind = where(finite(data_ext.EXT_U) EQ 0 or finite(data_ext.sigextUU) EQ 0, countnokuex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value
IF countnokuex NE 0 THEN BEGIN
message, "SED format isn't met. NaN value(s) detected in EXT_U data.",/continue
message, 'Without automatic format fitering you cannot proceed.',/continue
ans=''
read, ans, prompt='Would you like automatic format filtering? (Y/N):'
IF strupcase(ans) EQ 'N' THEN BEGIN
message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue
message, 'Aborting ...'
ENDIF
ENDIF
IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $
ELSE message, 'Data format is met. Removal of undefined data points...',/continue
tdaerruex=where((data_ext.sigextUU EQ la_undef() and data_ext.EXT_U NE la_undef()) or (data_ext.sigextUU NE la_undef() and data_ext.EXT_U EQ la_undef()),ctdaerruex)
IF ctdaerruex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...'
testuex=where(data_ext.sigextUU NE 0. and finite(data_ext.sigextUU) EQ 1 and data_ext.sigextUU NE la_undef() and data_ext.EXT_U NE la_undef() and finite(data_ext.EXT_U) EQ 1,ctestuex)
IF ctestuex NE 0 THEN BEGIN
instru_names = uext.instru_names(testuex)
filt_names = uext.filt_names(testuex)
wav = uext.wav(testuex)
values = uext.values(testuex)
sigma = uext.sigma(testuex)
uext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
ENDIF
ENDELSE
psiexzone: ;Starting the analogy here.
if NU_EXT eq 0 and NQ_EXT eq 0 then goto, pexzone
degtorad = !pi/180
;Generation of POLEXT and PSI_EXT data by default
;Testing on the presence of covariances otherwise an error is returned
;The testing happens here because covariances are needed for the error estimation of EXT_P,PSI_EXT.
tstcovuex = where((data_ext(testuex).sigextIU) EQ la_undef() or (data_ext(testuex).sigmaQU) EQ la_undef(), ctcovuex)
if ctcovuex NE 0 then message, 'Covariance data is missing in the extinction data pertaining to the StokesU parameter. Data, errors and covariances do not match. Aborting...'
tstcovqex = where((data_ext(testqex).sigextIQ) EQ la_undef() or (data_ext(testqex).sigextQU) EQ la_undef(), ctcovqex)
if ctcovqex NE 0 then message, 'Covariance data is missing in the extinction pertaining to the StokesQ parameter. Data, errors and covariances do not match. Aborting...'
;Matching both filter and spectrum points
match2,uext.wav,qext.wav,qexpts,uexpts
;removing -1 values from the indices (qpts & upts)
indnotqex = where(qexpts eq -1, cnotqex)
if cnotqex ne 0 then remove, indnotqex, qexpts
indnotuex = where(uexpts eq -1, cnotuex)
if cnotuex ne 0 then remove, indnotuex, uexpts
match2, ext.wav, qexsed.wav(qexpts), qextmp, extpts
indnotx = where(extpts eq -1, cnotx)
if cnotx ne 0 then remove, indnotx, extpts
instru_names = uext.instru_names(uexpts)
filt_names = uext.filt_names(uexpts)
wav = uext.wav(uexpts)
values = 0.5*atan(uext.values(uexpts),qext.values(qexpts))/degtorad
polar_variance_iqu2ippsi,ext.values(sedpts),qexsed.values(qexpts),uext.values(uexpts),(data_ext(testx).sigextII)(extpts),(data_ext(testqex).sigextQQ)(qexpts),(data_ext(testuex).sigextUU)(uexpts),(data_ext(testx).sigextIQ)(extpts),(data_ext(testx).sigextIU)(extpts),(data_ext(testx).sigextQU)(extpts),variance_smallp,variance_psi,variance_largep
sigma = sqrt(variance_psi)
psi_ext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma}
goto, pexzone ; the variances for polextPSI_EXT were already computed.
pexzone:
polext = psi_ext
polext.values = sqrt(uext.values(uexpts)^2+qexsed.values(qexpts)^2)
polext.sigma = sqrt(variance_largep)
goto, end_ext
end_ext:
ENDIF
;stop
RETURN, 0.
END