dustem_compute_polext.pro
4.6 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
FUNCTION dustem_compute_polext ,p_dim,sti,POLEXT_spec,out_st=out_st,dustem_qext,dustem_uext,Q_ext,U_ext,SPEXT_spec,dustem_fpolext,dustem_ext=dustem_ext,_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=dustem_run(p_dim)
ENDIF
POLEXT_spec = sti.polext.ext_tot * (*!dustem_HCD)/1.0e21 ;
out_st = sti
EXT_spec = st.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
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]
SPEXT_spec=(*(*!dustem_plugin).(i).spec)[*,3]
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]
SPEXT_spec+=(*(*!dustem_plugin).(i).spec)[*,3]
ENDIF
ENDFOR
ENDIF
POLEXT_spec=sqrt(QEXT_spec^2+UEXT_spec^2)
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,st.polext.wav,(((*!dustem_data.qext).wav)(ind_spec)))
out_st=st
;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,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