dustem_compute_polsed.pro
4.08 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
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
dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values)
stp=dustem_run(p_dim)
ENDIF
If keyword_set(result) THEN result=stp
fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/(stp.polsed).wav)*1.e20/1.e7
P = stp.polsed.em_tot*fact;1.00e17/(4.*!pi) ;This is Polarized intensity P
stI = stp.sed.em_tot*fact;1.00e17/(4.*!pi) ;This is Total intensity I
out=0.
Nwaves=(size(P))[1]
frac=P/stI
tes=where(finite(frac) eq 0)
frac(tes)=0.
;IF ptr_valid(!dustem_plugin) THEN BEGIN
scopes=tag_names((*!dustem_scope))
IF scopes[0] NE 'NONE' THEN BEGIN
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)
ENDIF
IF isa(!dustem_data.polsed) THEN dustem_polsed = (*!dustem_data.polsed).values * 0. ELSE dustem_polsed = (*!dustem_data.polfrac).values * 0.
;stop
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
IF isa(!dustem_data.polsed) THEN ind_polsed=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_polsed) ELSE ind_polsed=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_polsed)
IF count_polsed NE 0 THEN BEGIN
IF isa(!dustem_data.polsed) THEN filter_names=((*!dustem_data.polsed).filt_names)(ind_polsed) ELSE filter_names=((*!dustem_data.polfrac).filt_names)(ind_polsed)
spolsed=dustem_cc(stp.polsed.wav,P,filter_names,cc=cc)
dustem_polsed[ind_polsed]=spolsed
ENDIF
;stop
;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).
IF isa(!dustem_data.polsed) THEN ind_spec=where((*!dustem_data.polsed).filt_names EQ 'SPECTRUM',count_spec) ELSE ind_spec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',count_spec)
IF count_spec NE 0 THEN BEGIN
IF isa(!dustem_data.polsed) THEN dustem_polsed(ind_spec)=interpol(P,stp.polsed.wav,(((*!dustem_data.polsed).wav)(ind_spec))) ELSE dustem_polsed(ind_spec)=interpol(P,stp.polsed.wav,(((*!dustem_data.polfrac).wav)(ind_spec)))
ENDIF
out_st=stp
;clean pointers
heap_gc
the_end:
RETURN,dustem_polsed
END