dustem_compute_polext.pro
6.52 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
FUNCTION dustem_compute_polext,p_dim,$
sti=sti,$
POLEXT_spec=POLEXT_spec,$
SPEXT_spec=SPEXT_spec,$
dustem_fpolext=dustem_fpolext,$
dustem_ext=dustem_ext,$
out_st=out_st,$
_extra=extra, $
help=help
;+
; NAME:
; dustem_compute_polext
; PURPOSE:
; Computes Polarized Extinction in units of optical depth for a given Dustem spectrum
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; polext=dustem_compute_polext(p_dim[,sti=][,polext_spec=][,spext_spec=][,dustem_fpolext=][,dustem_ext=][,out_st=][,_extra=][,/help])
; INPUTS:
; p_dim = parameter values
; OPTIONAL INPUT PARAMETERS:
; st = Dustem output structure
; dustem_ext = Computed polarized extinction for spectrum points in !dustem_data
; OUTPUTS:
; polext = computed Polarized Extinction for spectrum points in !dustem_data
; OPTIONAL OUTPUT PARAMETERS:
; out_st = Dustem output structure
; polext_spec = Dustem polarized extiction output
; spext_spec = Dustem polarized exinction fraction
; dustem_fpolext = Computed polarized extinction fraction for spectrum points in !dustem_data
; ACCEPTED KEY-WORDS:
; help = If set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The dustem idl wrapper must be installed
;-
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_compute_polext'
dustem_polext=0.
goto,the_end
ENDIF
IF not keyword_set(sti) THEN BEGIN
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=sti
ENDIF
POLEXT_spec = sti.polext.ext_tot * (*!dustem_HCD)/1.0e21 ;
;Hard-coded test to set to potential zero 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
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
EXT_spec = sti.ext.ext_tot * (*!dustem_HCD)/1.0e21 ;This is Total intensity I
out=0.
Nwaves=(size(POLEXT_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
;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_EXT') then EXT_spec+=(*(*!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_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)
SPEXT_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_POLEXT') THEN BEGIN
QEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,1]
UEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,2]
ENDIF
ENDFOR
ENDIF
POLEXT_spec=sqrt(QEXT_spec^2+UEXT_spec^2)
SPEXT_spec = POLEXT_SPEC/EXT_spec
dustem_polext = (*(*!dustem_data).qext).values * 0.
if not isarray(POLEXT_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
;NO COLOR CORRECTION FOR EXTINCTION. FOR NOW ...
; 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).qext).filt_names EQ 'SPECTRUM',count_spec)
IF count_spec NE 0 THEN dustem_polext(ind_spec)=interpol(POLEXT_spec,sti.polext.wav,(((*(*!dustem_data).qext).wav)(ind_spec)))
out_st=sti
;GENERATING THE INTERPOLATES FOR POLFRAC (dustem_polfrac)
if !run_lin then begin
if not keyword_set(dustem_ext) then dustem_ext = dustem_compute_ext(p_dim,st=sti,EXT_spec=EXT_spec)
If n_elements(dustem_ext) ne n_elements(dustem_polext) then begin
if n_elements(dustem_ext) gt n_elements(dustem_polext) then begin
dustem_polext_x = dustem_polext
dustem_ext_x = dustem_polext ;meaning dustem_sed needs to be modified
nwaves = n_elements(dustem_polext)
for i=0L,nwaves-1 do begin
j=where((*(*!dustem_data).ext).wav EQ (*(*!dustem_data).fpolext).wav(i),testwav)
if testwav ne 0 then dustem_ext_x(i) = dustem_ext(j(0))
endfor
endif ELSE begin
dustem_polext_x = dustem_ext
dustem_ext_x = dustem_ext ;meaning dustem_polsed needs to be modified
nwaves = n_elements(dustem_ext)
for i=0L,nwaves-1 do begin
j=where((*(*!dustem_data).ext).wav EQ (*(*!dustem_data).fpolext).wav(i),testwav)
if testwav ne 0 then dustem_polext_x(i) = dustem_polext(j(0))
endfor
ENDELSE
endif ELSE BEGIN
dustem_polext_x = dustem_polext
dustem_ext_x = dustem_ext
ENDELSE
dustem_fpolext = dustem_polext_x/dustem_ext_x
ENDIF
the_end:
return, dustem_polext
END