dustem_compute_polsed.pro
4.47 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
FUNCTION dustem_compute_polsed,p_dim,stp,out_st=out_st,_extra=extra
;+
; NAME:
; dustem_compute_polsed
; PURPOSE:
; Computes an SED from a given Dustem spectrum
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; sed=dustem_compute_polsed(p_dim[,st=][,_extra=][,/help])
; INPUTS:
; p_dim = parameter values
; OPTIONAL INPUT PARAMETERS:
; stp = Dustem structure
; OUTPUTS:
; sed : computed SED for filters in !dustem_data
; dustem_qsed : SED in Q
; dustem_used : SED in U
; 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_polsed'
dustem_polsed=0.
goto,the_end
ENDIF
IF not keyword_set(stp) THEN BEGIN
stp=dustem_run(p_dim)
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),stp
ENDIF
If keyword_set(result) THEN result=stp
fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(stp.polsed).wav)*1.e20/1.e7
P = stp.polsed.em_tot * fact ;This is Polarized intensity P
I = stp.sed.em_tot * fact ;This is Total intensity I
out=0.
Nwaves=(size(P))[1]
frac=P/I
tes=where(finite(frac) eq 0)
frac(tes)=0.
;
;
; ;ADDING PLUGIN TO SPECTRUM----------------
; FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
; IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_POLSED') THEN P+=(*(*!dustem_plugin).(i))[*,3] ;3 is for the polarized portion of the plugin because this is polsed
; ENDFOR
; ;-----------------------------------------
;
;
; ;MODIFYING PLUGIN----------------
; FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
; IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_POLSED') THEN BEGIN
; P=(*(*!dustem_plugin).(i))[*,3]
; P+=(*(*!dustem_plugin).(i))[*,3] ;3 is for the polarized portion of the plugin because this is polsed
; ENDIF
; ENDFOR
; ;-----------------------------------------
;if Q and U have already been computed elsewhere
; FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
;
; IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'MODIFY_DUST_POLAR') THEN BEGIN
;
; Q_spec=(*(*!dustem_plugin).(i))[*,1]
; U_spec=(*(*!dustem_plugin).(i))[*,2]
;
; ENDIF
; ENDFOR
FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_QSED') THEN BEGIN
Q_spec=(*(*!dustem_plugin).(i))[*,1]
ENDIF
IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_USED') THEN BEGIN
U_spec=(*(*!dustem_plugin).(i))[*,2]
ENDIF
ENDFOR
IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(0.,Nwaves)
FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_QSED') THEN BEGIN
Q_spec+=(*(*!dustem_plugin).(i))[*,1]
ENDIF
IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_USED') THEN BEGIN
U_spec+=(*(*!dustem_plugin).(i))[*,2]
ENDIF
ENDFOR
P=sqrt(Q_spec^2+U_spec^2)
dustem_polsed = (*!dustem_data.sed).values * 0.
if not isarray(P) THEN stop
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.polsed).filt_names NE 'SPECTRUM',count_polsed)
IF count_polsed NE 0 THEN BEGIN ; INCLUDE COLOR CORRECTIONS OF (Q,U) for dust & synch
filter_names=((*!dustem_data.polsed).filt_names)(ind_polsed)
spolsed=dustem_cc(stp.sed.wav,P,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).
ind_spec=where((*!dustem_data.polsed).filt_names EQ 'SPECTRUM',count_spec)
IF count_spec NE 0 THEN BEGIN
dustem_polsed(ind_spec)=interpol(P,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
ENDIF
out_st=stp
;clean pointers
heap_gc
the_end:
RETURN,dustem_polsed
END