dustem_compute_stokext.pro
4.79 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
FUNCTION dustem_compute_stokext,p_dim,st,dustem_qext,dustem_uext,QEXT_spec,UEXT_spec,PSIEXT_spec,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=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.d21 ;This is Total extinction (for a default value of 10^20 NH/cm^2)
POLEXT_spec = st.polext.ext_tot * (*!dustem_HCD)/1.d21 ;This is Polarized extinction
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]
PSIEXT_spec = (*(*!dustem_plugin).(i).spec)[*,4]
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
PSIEXT_spec = 0.5*atan(UEXT_spec,QEXT_spec)/degtorad
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]
PSIEXT_spec+= (*(*!dustem_plugin).(i).spec)[*,4]
ENDIF
ENDFOR
ENDIF
;-----------------------------------------
;INITIALIZING THE STOKES SEDS
dustem_qext = (*!dustem_data.qext).values * 0.
dustem_uext = (*!dustem_data.uext).values * 0.
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
;clean pointers
heap_gc
the_end:
RETURN,out
END