dustem_compute_polext.pro
5.09 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
FUNCTION dustem_compute_polext ,p_dim,sti,POLEXT_spec,SPEXT_spec,dustem_fpolext,dustem_ext=dustem_ext,out_st=out_st,_extra=extra
;CREATED AS A CONSEQUENCE OF THE INCLUSION OF THE PLUGINS
;written to homogonize DustEM
IF not keyword_set(sti) THEN BEGIN
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),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,sti,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
return, dustem_polext
END