dustem_compute_polsed.pro
6.18 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
FUNCTION dustem_compute_polsed,p_dim,st,P_spec,SP_spec,dustem_polfrac,out_st=out_st,dustem_sed=dustem_sed,_extra=extra
;+
;==========JPB
;THIS HELP IS NOT UP TO DATE
;===========
; NAME:
; dustem_compute_polsed
; PURPOSE:
; Computes an SED from a given Dustem spectrum
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; sed=dustem_compute_polsed(p_dim[,st=][,_extra=][,/help])
; INPUTS:
; p_dim = parameter values
; OPTIONAL INPUT PARAMETERS:
; st = Dustem structure
; OUTPUTS:
; sed : computed SED for filters in !dustem_data
; dustem_qsed : SED in Q
; dustem_used : SED in U
; OPTIONAL OUTPUT PARAMETERS:
; st = Dustem output structure
; ACCEPTED KEY-WORDS:
; help = If set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The dustem idl wrapper must be installed
; PROCEDURE:
; None
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by J.-Ph. Bernard
; see evolution details on the dustem cvs maintained at CESR
; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_compute_polsed'
dustem_polsed=0.
goto,the_end
ENDIF
IF not keyword_set(st) THEN BEGIN
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st
;st=dustem_run(p_dim)
ENDIF
If keyword_set(result) THEN result=st
fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.e20/1.e7
P_spec = st.polsed.em_tot*fact ;This is Polarized intensity P
stI = st.sed.em_tot*fact ;This is Total intensity I
out=0.
Nwaves=(size(P_spec))[1]
frac=P_spec/stI
tes=where(finite(frac) eq 0)
frac(tes)=0.
;This block was added to generate the polarization fraction spectrum with the plugins contributions.
;smallp (polfrac) is naturally not an additive quantity.
;NB: if plugins that replace the total emission given by the fortran executable exist,
;the need to add a block with a 'REPLPACE_SED' scope as for the 'REPLACE_POLSED' below, may arise.
scopes=tag_names((*!dustem_plugin))
IF scopes[0] NE 'NONE' THEN BEGIN
;IF ptr_valid(!dustem_plugin) THEN BEGIN
for i=0L,n_tags(*!dustem_plugin)-1 do begin
if total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) eq 'ADD_SED') then stI+=(*(*!dustem_plugin).(i).spec)[*,0]
endfor
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 'REPLACE_POLSED') THEN BEGIN
Q_spec=(*(*!dustem_plugin).(i).spec)[*,1]
U_spec=(*(*!dustem_plugin).(i).spec)[*,2]
;make sure sti hasn't been replaced. But this is not happening in the near future
ENDIF
ENDFOR
ENDIF
IF ~isa(Q_spec) && ~isa(U_spec) THEN BEGIN
polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(!dustem_psi,Nwaves) ;the user needs to set this pointer somewhere
SP_spec = frac
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
P_spec=sqrt(Q_spec^2+U_spec^2)
SP_spec=P_spec/Sti
dustem_polsed = (*(*!dustem_data).qsed).values * 0.
if not isarray(P_spec) THEN stop ;I don't understand this test. The only thing it can indicate is a problem with the dust parameters or the fortran executable. But in my opinion these tests would be
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_polsed=where((*(*!dustem_data).qsed).filt_names NE 'SPECTRUM',count_polsed)
IF count_polsed NE 0 THEN BEGIN
filter_names=((*(*!dustem_data).qsed).filt_names)(ind_polsed)
spolsed=dustem_cc(st.polsed.wav,P_spec,filter_names,cc=cc)
dustem_polsed[ind_polsed]=spolsed
ENDIF
;For spectrum data points, interpolate in log-log.
;Linear interpolation leads to wrong values, in particular where few
;wavelengths points exist in the model (long wavelengths).
;I believe this to be automatic because the provided dustem spectra are already sampled on a log-log grid.
ind_spec=where((*(*!dustem_data).qsed).filt_names EQ 'SPECTRUM',count_spec)
IF count_spec NE 0 THEN dustem_polsed(ind_spec)=interpol(P_spec,st.polsed.wav,(((*(*!dustem_data).qsed).wav)(ind_spec)))
out_st=st
;GENERATING THE INTERPOLATES FOR POLFRAC (dustem_polfrac)
if !run_lin then begin
if not keyword_set(dustem_sed) then dustem_sed = dustem_compute_sed(p_dim,st,SED_spec)
If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin
if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin
dustem_polsed_x = dustem_polsed
dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified
nwaves = n_elements(dustem_polsed)
for i=0L,nwaves-1 do begin
j=where((*(*!dustem_data).sed).wav EQ (*(*!dustem_data).polfrac).wav(i),testwav)
if testwav ne 0 then dustem_sed_x(i) = dustem_sed(j(0))
endfor
endif ELSE begin
dustem_polsed_x = dustem_sed
dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified
nwaves = n_elements(dustem_sed)
for i=0L,nwaves-1 do begin
j=where((*(*!dustem_data).sed).wav EQ (*(*!dustem_data).polfrac).wav(i),testwav)
if testwav ne 0 then dustem_polsed_x(i) = dustem_polsed(j(0))
endfor
ENDELSE
endif ELSE BEGIN
dustem_polsed_x = dustem_polsed
dustem_sed_x = dustem_sed
ENDELSE
dustem_polfrac = dustem_polsed_x/dustem_sed_x
ENDIF
;clean pointers
heap_gc
the_end:
RETURN,dustem_polsed
END