dustem_compute_polsed.pro
4.54 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_polsed,p_dim,stp,_extra=extra,out_st=out_st,dustem_qsed,dustem_used,Q_sed,U_sed
;+
; 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:
; None
; OUTPUTS:
; sed = computed SED for filters in !dustem_data
; OPTIONAL OUTPUT PARAMETERS:
; st = 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),res=res,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
spec = stp.polsed.em_tot * fact
;stop
;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 spec+=(*(*!dustem_plugin).(i))[*,0]
endfor
;-----------------------------------------
if not isarray(spec) THEN stop
;stop
;COMPUTE THE MODEL SED
dustem_polsed = (*!dustem_data.polsed).values * 0.
;MODIFYING VIA PLUGIN-------------- THIS BLOCK NEEDS TO BE UNDER THE DUSTEM_POLSED LINE
for i=0L,n_tags(*!dustem_scope)-1 do begin
if total(strsplit((*(*!dustem_scope).(i)),'+',/extract) eq 'MODIFY_POLSED') then begin
;making sure it's the stokes plugin that gets used here
iidx=where(tag_names(*!dustem_scope) eq 'SYNCHROTRON',countiidx)
IF countiidx ne 0 THEN BEGIN
Q_dust=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/4-1)
U_dust=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/4:n_elements((*(*!dustem_plugin).(i)))/2-1)
Q_synch=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:3*n_elements((*(*!dustem_plugin).(i)))/4-1)
U_synch=(*(*!dustem_plugin).(i))(3*n_elements((*(*!dustem_plugin).(i)))/4:*)
ENDIF ELSE BEGIN
Q_dust=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/2-1)
U_dust=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:*)
Q_synch=Q_dust*0
U_synch=U_dust*0
ENDELSE
;THE ABOVE VARIABLES SHOULD BE ACCESSED DIRECTLY BY THE MPFIT_RUN.PRO
Q_sed=Q_dust*fact+Q_synch
U_sed=U_dust*fact+U_synch
;data the dustem_wrap will use
dustem_qsed = dustem_polsed
dustem_used = dustem_polsed
endif
endfor
;-------------------------------------------
;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,spec,filter_names,cc=cc) & dustem_polsed(ind_polsed)=spolsed
;PLUGIN 'STOKES' W/ scope='MODIFY' PART------------------------------------
if isa(dustem_qsed) then begin
sqsed=dustem_cc(stp.sed.wav,Q_sed,filter_names,cc=cc) & dustem_qsed(ind_polsed)=sqsed
endif
if isa(dustem_used) then begin
sused=dustem_cc(stp.sed.wav,U_sed,filter_names,cc=cc) & dustem_used(ind_polsed)=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_spec=where((*!dustem_data.polsed).filt_names EQ 'SPECTRUM',count_spec)
IF count_spec NE 0 THEN BEGIN
dustem_polsed(ind_spec)=interpol(spec,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
;STOKES PARAMS TREATMENT
if isa(dustem_qsed) then dustem_qsed(ind_spec)=interpol(Q_sed,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
if isa(dustem_used) then dustem_used(ind_spec)=interpol(U_sed,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
ENDIF
out_st=stp
;clean pointers
heap_gc
the_end:
RETURN,dustem_polsed
END