dustem_make_polarization_sed_example.pro
7.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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
PRO dustem_make_polarization_sed_example,model=model,help=help
;+
; NAME:
; dustem_make_polarization_sed_example
; PURPOSE:
; This is an example of how to generate an SED using dustemwrap.
; It is meant to be an example to follow when writing your own
; programs using the dustem IDL/GDL wrapper.
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; dustem_make_polarization_sed_example[,model=][,/help]
; INPUTS:
; None
; OPTIONAL INPUT PARAMETERS:
; None
; OUTPUTS:
; None
; OPTIONAL OUTPUT PARAMETERS:
; None
; ACCEPTED KEY-WORDS:
; model = if set, name of the Dustem dust model to use (default='G17_MODELC')
; help = If set print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; None
; RESTRICTIONS:
; The dustem idl wrapper must be installed
; PROCEDURE:
; None
; EXAMPLES
; dustem_make_polarization_sed_example
; dustem_make_polarization_sed_example,model='G17_MODELD'
; MODIFICATION HISTORY:
; Written by J.P. Bernard June 29th 2011
; 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_fit_sed_polsed_readme'
goto,the_end
ENDIF
dustem_define_la_common
;=== initialise dustem
;dustem_init,mode='COMPIEGNE_ETAL2010',/pol
;use_mode='G17_MODELA'
use_mode='G17_MODELC'
dustem_init,mode=use_mode,kwords=['logn-chrg-spin','plaw-pol','plaw-pol'],/pol
!dustem_verbose=1
!dustem_show_plot=1
!dustem_nocatch=1
filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
Nfilt=n_elements(filters)
sed=dustem_initialize_sed(Nfilt)
sed.filter=filters
sed.wave=dustem_filter2wav(filters)
sed.instru=dustem_filter2instru(filters)
;INITIALIZING IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro
sed[*].StokesI=1.
sed.StokesQ=sed.StokesI/100.
sed.StokesU=sed.StokesI/100.
sed.SigmaII=sed.StokesI/100.
sed.SigmaQQ=sed.StokesI/100.
sed.SigmaUU=sed.StokesI/100.
sed.SigmaIQ=sed.StokesI/100.
sed.SigmaIU=sed.StokesI/100.
sed.SigmaQU=sed.StokesI/100.
;INITIALIZING THE !dustem_data tags that will be internally used
dustem_set_data,sed
help,!dustem_data
;=== Set which model parameters to use for the SED
;=== This does not work, for some reason ....
pd = [ $
'(*!dustem_params).G0', $ ;G0
'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(2).mdust_o_mh', $ ;PAH1 mass fraction
'dustem_plugin_modify_dust_pol_1', $ ;This will set the polarization fraction
'dustem_plugin_modify_dust_pol_2' $ ;This will set the polarization angle
]
p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,0.2,45.]
;=== Set which model parameters to use for the SED
pd = [ $
'(*!dustem_params).G0', $ ;G0
'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(2).mdust_o_mh' $ ;PAH1 mass fraction
]
p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04]
pd = [ $
'(*!dustem_params).G0', $ ;G0
'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(2).mdust_o_mh', $ ;PAH1 mass fraction
'dustem_plugin_synchrotron_2', $ ;Synchrotron amplitude at 10 mm
'dustem_plugin_synchrotron_3', $ ;Synchrotron polarization angle
'dustem_plugin_synchrotron_4' $ ;Synchrotron polarization fraction
]
p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04, 1., 45., 0.3]
;===== This example adds synchrotron with a polarization angle of 45 degree
pd = [ $
'(*!dustem_params).G0', $ ;G0
'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(2).mdust_o_mh', $ ;PAH1 mass fraction
'dustem_plugin_synchrotron_1', $ ;Synchrotron spectra index
'dustem_plugin_synchrotron_2', $ ;Synchrotron amplitude at 10 mm
'dustem_plugin_synchrotron_3', $ ;Synchrotron polarization fraction
'dustem_plugin_synchrotron_4' $ ;Synchrotron polarization angle
]
p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,3., 1.e-2, 0.3, 45.]
pd = [ $
'(*!dustem_params).G0', $ ;G0
'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH0 mass fraction
'(*!dustem_params).grains(2).mdust_o_mh', $ ;PAH1 mass fraction
'dustem_plugin_modify_dust_pol_2', $ ;This will set the polarization angle
'dustem_plugin_synchrotron_1', $ ;Synchrotron spectra index
'dustem_plugin_synchrotron_2', $ ;Synchrotron amplitude at 10 mm
'dustem_plugin_synchrotron_3', $ ;Synchrotron polarization fraction
'dustem_plugin_synchrotron_4' $ ;Synchrotron polarization angle
]
p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,10.,3., 1.e-2, 0.3, 45.]
;iv=p_truth+[0.,8.e-4,8.e-4,8.e-4] ;shifted from solution
iv=p_truth
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
llims=replicate(0,Npar)
dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
dustem_init_plugins,pd
;Here fortran fails when pluggins are activated and/or type-keywods exist
dustem_Ised=dustem_compute_sed(p_truth,ssti) ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
toto=dustem_compute_stokes(p_truth,sstqu,dustem_Qsed,dustem_Used) ;this procedure also allows for the extraction of the spectra
;toto=dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em,out_st=out_st,_extra=extra)
print,dustem_Qsed,dustem_Used
sed.stokesI=dustem_Ised
sed.stokesQ=dustem_Qsed
sed.stokesU=dustem_Used
;goto,no_change
;=modifies the Fortran outputs Q,U SED values to reflect a fixed polarization fraction and angle
;normaly, this could also be doen using the plugin:
;toto=dustem_plugin_modify_dust_pol(ssti, key=key, val=val, scope=scope, paramtag=paramtag, help=help)
;frac_used=0.01 ;1%
;psi_used=10. ;10 degrees
;smallp=dustem_Ised*0.+frac_used
;psi=dustem_Ised*0.+psi_used
;polar_ippsi2iqu,dustem_Ised,Q,U,smallp,psi
;sed.stokesI=dustem_Ised
;sed.stokesQ=Q
;sed.stokesU=U
;no_change:
;=== set arbitrary uncertainties
sed.sigmaII=abs(sed.StokesI*0.000001)
sed.sigmaQQ=abs(sed.StokesQ*0.000001)
sed.sigmaUU=abs(sed.StokesU*0.000001)
sed.sigmaIQ=0.
sed.sigmaIU=0.
sed.sigmaQU=0.
print,sed.stokesI,sed.stokesQ,sed.stokesU
print,sed.sigmaII,sed.sigmaQQ,sed.sigmaUU
;==== fill in dependent columns of the SED.
sed=dustem_fill_sed_dependent_columns(sed)
;cgplot,sed.wave,sed.stokesI,psym=-3,/xlog,/ylog
;stop
;======== save the SED
filename='/tmp/fake_polarization_sed.xcat'
write_xcat,sed,filename
message,'Wrote '+filename,/continue
filename_final=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/example_sed3.xcat'
str='cp '+filename+' '+filename_final
message,'Do '+str+' to make change permanent',/continue
;spawn,str
;message,'Wrote '+filename_final,/continue
the_end:
END