dustem_compute_stokes.pro
5.73 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
FUNCTION dustem_compute_stokes,p_dim,$
st=st,$
Q_spec=Q_spec,$
U_spec=U_spec,$
PSI_spec=PSI_spec,$
dustem_psi_em=dustem_psi_em,$
help=help
;+
; NAME:
; dustem_compute_stokes
;
; PURPOSE:
; Computes Stokes Parmeters SEDs for a given Dustem spectrum
;
; CATEGORY:
; Dustem
;
; CALLING SEQUENCE:
; out=dustem_compute_stokes(p_dim[,st=][,Q_spec=][,U_spec=][,PSI_spec=][,dustem_psi_em=][,/help])
;
; INPUTS:
; p_dim = parameter values
;
; OPTIONAL INPUT PARAMETERS:
; st = Dustem structure
;
; OUTPUTS:
; out = List Containing the Stokes parameters emission for spectrum points in !dustem_data
; out[0] = Computed Stokes Q emission for spectrum points in !dustem_data
; out[1] = Computed Stokes U emission for spectrum points in !dustem_data
;
; OPTIONAL OUTPUT PARAMETERS:
; st = Dustem output structure
; Q_spec = Dustem Stokes Q emission
; U_spec = Dustem Stokes U emission
; PSI_spec = Dustem Polarization angle output in emission
; dustem_psi_em = Computed Polarization angle ouput for spectrum points in !dustem_data
;
; ACCEPTED KEY-WORDS:
; help = If set, print this help
;
; COMMON BLOCKS:
; None
;
; SIDE EFFECTS:
; None
;
; RESTRICTIONS:
; The DustEMWrap IDL code must be installed
;
; PROCEDURES:
; None
;
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by J.-Ph. Bernard 2021
; Evolution details on the DustEMWrap gitlab
; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_compute_stokes'
out=0.
goto,the_end
ENDIF
IF not keyword_set(st) THEN $
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st
fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.e20/1.e7
stI = st.sed.em_tot * fact ;This is Total intensity I
P = st.polsed.em_tot * fact ;This is Polarized intensity P
Nwaves=(size(stI))[1]
frac=P/stI
tes=where(finite(frac) eq 0)
frac(tes)=0.
scopes=tag_names((*!dustem_plugin))
IF scopes[0] NE 'NONE' THEN BEGIN
FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN
IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'REPLACE_POLSED') THEN BEGIN
Q_spec = (*(*!dustem_plugin).(i).spec)[*,1]
U_spec = (*(*!dustem_plugin).(i).spec)[*,2]
ENDIF
ENDFOR
ENDIF
IF ~isa(Q_spec) && ~isa(U_spec) THEN BEGIN
polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(!dustem_psi,Nwaves);
ENDIF
IF scopes[0] NE 'NONE' THEN BEGIN
FOR i=0L,n_tags(*!dustem_plugin)-1 DO BEGIN
IF total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) EQ 'ADD_POLSED') THEN BEGIN
Q_spec+=(*(*!dustem_plugin).(i).spec)[*,1]
U_spec+=(*(*!dustem_plugin).(i).spec)[*,2]
ENDIF
ENDFOR
ENDIF
;-----------------------------------------
ind=where(finite(U_spec) NE 1 OR finite(Q_spec) NE 1,count)
IF count NE 0 THEN BEGIN
message,'Q or U is NaN',/continue
stop
ENDIF
;Hard-coded block so that atan(U,Q) does not return wrong values when Q eq -0
testnzero = where(Q_spec EQ -0., ct_nzero)
if ct_nzero ne 0 then Q_spec[testnzero] = 0
testnzero = where(U_spec EQ -0., ct_nzero)
if ct_nzero ne 0 then U_spec[testnzero] = 0
PSI_spec = 0.5*atan(U_spec,Q_spec)/!dtor
;polar_iqu2ippsi,stI,Q_spec,U_spec,p_bidon,PSI_spec ;same but this line runs several commands for each iteration of this current procedure
;INITIALIZING THE STOKES SEDS
dustem_qsed = (*(*!dustem_data).qsed).values * 0.
dustem_used = (*(*!dustem_data).used).values * 0.
;using a list because the arrays can have different dimensions
;this doesn't change anything ...
out=list(fltarr(n_elements(dustem_qsed)),fltarr(n_elements(dustem_used)))
if not isarray(stI) THEN stop ; I don't get this test. It is unnecessary.
;Performing color corrections on the dustem spectra and generating the SEDs at the given filters
IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
message,'DOING color correction calculations',/info
ENDIF ELSE BEGIN
message,'SKIPPING color correction calculations',/info
ENDELSE
ind_qsed=where((*(*!dustem_data).qsed).filt_names NE 'SPECTRUM',count_qsed)
ind_used=where((*(*!dustem_data).used).filt_names NE 'SPECTRUM',count_used)
IF count_qsed NE 0 THEN BEGIN
filter_names=((*(*!dustem_data).qsed).filt_names)[ind_qsed]
if isa(dustem_qsed) then begin
sqsed=dustem_cc(st.polsed.wav,Q_spec,filter_names,cc=cc)
dustem_qsed[ind_qsed]=sqsed
endif
ENDIF
IF count_used NE 0 THEN BEGIN
filter_names=((*(*!dustem_data).used).filt_names)[ind_used]
if isa(dustem_used) then begin
sused=dustem_cc(st.polsed.wav,U_spec,filter_names,cc=cc)
dustem_used[ind_used]=sused
endif
ENDIF
;For spectrum data points, interpolate in log-logdustem_compute_stokes.pro
;Linear interpolation leads to wrong values, in particular where few
;wavelengths points exist in the model (long wavelengths).
ind_qspec=where((*(*!dustem_data).qsed).filt_names EQ 'SPECTRUM',count_qspec)
ind_uspec=where((*(*!dustem_data).used).filt_names EQ 'SPECTRUM',count_uspec)
IF count_qspec NE 0 THEN BEGIN
if isa(dustem_qsed) then dustem_qsed[ind_qspec]=interpol(Q_spec,st.polsed.wav,(((*(*!dustem_data).qsed).wav)[ind_qspec]))
ENDIF
IF count_uspec NE 0 THEN BEGIN
if isa(dustem_used) then dustem_used[ind_uspec]=interpol(U_spec,st.polsed.wav,(((*(*!dustem_data).used).wav)[ind_uspec]))
ENDIF
;generating interpolates for Psi_em
dustem_psi_em = 0.5*atan(dustem_used,dustem_qsed)/!dtor
out[0]=dustem_qsed
out[1]=dustem_used
;clean pointers
heap_gc
the_end:
RETURN,out
END