dustem_compute_stokext.pro
5.53 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
FUNCTION dustem_compute_stokext,p_dim,st,$
QEXT_spec=QEXT_spec,$
UEXT_spec=UEXT_spec,$
PSIEXT_spec=PSIEXT_spec,$
dustem_psi_ext=dustem_psi_ext,$
out_st=out_st,$
_extra=extra
;+
; NAME:
; dustem_compute_stokext
; PURPOSE:
; Provides Stokes Q and Stokes U "extinction" optical depths for a given dustem spectrum and wrapper SED
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; out=dustem_compute_stokext(p_dim[,st=][,_extra=][,/help])
; INPUTS:
; p_dim = parameter values
; OPTIONAL INPUT PARAMETERS:
; stp = Dustem structure
; OUTPUTS:
;
;
;
;
; OPTIONAL OUTPUT PARAMETERS:
; stp = 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_stokext'
out=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
out=0. ; no specific output is required.
If keyword_set(result) THEN result=st
degtorad = !pi/180
EXT_spec = st.ext.ext_tot * (*!dustem_HCD)/1.e21 ;This is Total extinction (for a default value of 10^20 NH/cm^2)
POLEXT_spec = st.polext.ext_tot * (*!dustem_HCD)/1.e21 ;This is Polarized extinction
;Hard-coded test to set to zero potential negative values in the polarized extinction arrays:
;so far this has been seen in absorption arrays
;This test needs to include all grain species that polarize
IF !run_pol THEN BEGIN
;THIS FALSE CORRECT IT USING POLEXT_spec
;POLEXT_spec.abs_grain[2,*]
ind_ngtv_plxt = where (POLEXT_spec LT 0, ct_ngtv_plxt)
IF ct_ngtv_plxt NE 0 THEN (POLEXT_spec)[ind_ngtv_plxt] = 0.
ENDIF
Nwaves=(size(EXT_SPEC))[1]
frac=POLEXT_spec/EXT_spec
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_POLEXT') THEN BEGIN
QEXT_spec = (*(*!dustem_plugin).(i).spec)[*,1]
UEXT_spec = (*(*!dustem_plugin).(i).spec)[*,2]
ENDIF
ENDFOR
ENDIF
IF ~isa(QEXT_spec) && ~isa(UEXT_spec) THEN BEGIN
polar_ippsi2iqu,EXT_spec,QEXT_spec,UEXT_spec,frac,replicate(!dustem_psi_ext,Nwaves);PSI_spec
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_POLEXT') THEN BEGIN
QEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,1]
UEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,2]
ENDIF
ENDFOR
ENDIF
;-----------------------------------------
;Hard-coded block to avoid cases where Q=-0.0 and U=0.
testnzero = where(QEXT_spec EQ -0., ct_nzero)
if ct_nzero ne 0 then QEXT_spec[testnzero] = 0
PSIEXT_spec = 0.5*atan(UEXT_spec,QEXT_spec)/degtorad
;INITIALIZING THE STOKES SEDS
dustem_qext = (*(*!dustem_data).qext).values * 0.
dustem_uext = (*(*!dustem_data).uext).values * 0.
;Using a list because dustem_qext and dustem_uext can have different dimensions
out=list(fltarr(n_elements(dustem_qext)),fltarr(n_elements(dustem_uext)))
if not isarray(EXT_spec) THEN stop
;COLOR CORRECTION NOT PERFORMED YET. ONLY SPECTRUM DATA POINTS ARE CONSIDERED FOR NOW
;
; ;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-log
;Linear interpolation leads to wrong values, in particular where few
;wavelengths points exist in the model (long wavelengths).
ind_qspec=where((*(*!dustem_data).qext).filt_names EQ 'SPECTRUM',count_qspec)
ind_uspec=where((*(*!dustem_data).uext).filt_names EQ 'SPECTRUM',count_uspec)
IF count_qspec NE 0 THEN BEGIN
if isa(dustem_qext) then dustem_qext(ind_qspec)=interpol(QEXT_spec,st.polext.wav,(((*(*!dustem_data).qext).wav)(ind_qspec)))
ENDIF
IF count_uspec NE 0 THEN BEGIN
if isa(dustem_uext) then dustem_uext(ind_uspec)=interpol(UEXT_spec,st.polext.wav,(((*(*!dustem_data).uext).wav)(ind_uspec)))
ENDIF
out_st=st
;generating interpolates for Psi_ext
dustem_psi_ext = 0.5*atan(dustem_uext,dustem_qext)/degtorad
out[0] = dustem_qext
out[1] = dustem_uext
;clean pointers
heap_gc
the_end:
RETURN,out
END