dustem_fit_sed_orion.pro
10.8 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
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
PRO dustem_fit_sed_orion,ps=ps,png=png
dustem_define_la_common
;=== initialise dustem
use_mode='G17_MODELC';'DBP90';'G17_MODELA';'THEMIS';;'USER_MODEL';'THEMIS';;'THEMIS';'G17_MODELA';'THEMIS';'NY_MODELA';'G17_MODELA';'DBP90';'USER_MODEL';'THEMIS';'NY_MODELA';'G17_MODELA';'DBP90';'COMPIEGNE_ETAL2010'
;,/pol
dustem_init,mode=use_mode,kwords=['logn-chrg-spin','plaw-pol','plaw-pol'],/pol;,kwords=['plaw-ed','logn','logn'];,/pol;kwords=['plaw-ed-chrg-spin'],/pol;,kwords=['plaw-ed'];kwords=['plaw-ed-chrg-spin'];,/pol;
;
;stop
!dustem_verbose=0
!dustem_show_plot=1
!dustem_nocatch=1
;specify level
lev=['300.000','400.000','500.000','600.000','800.000','1000.00','1200.00','1500.00','3000.00','5000.00','10000.0','15000.0']
indxlev=6 ; meaning level 400, 1 is for 600 etc...;THIS IS FOR THE PLOTTING OF ONE LEVEL ONLY
msk=['mask1','mask2','mask3'] ;this is used to be able to plot the new data that JP generated.
;1,2,3 correspond to the North-West, Center and South-East part of the OrionA integral shape filament
msk=msk(1) ;mask2 center of ORION_A
plckprd=['PR3','PR4','SROLL2']
plckprd=plckprd(1) ;PR4
polangmatch=['ang5','ang10','ang15']
polangmatch=polangmatch(1) ;best match. Lowest difference.
co_corr = ['noCO0','noCO1','noCO2']
co_corr = co_corr(2)
;opening corresponding sed file
smooth_str='_smoothedLFI3'
; dir_seds='/Users/ilyeschoubani/Downloads/SC110012_V2p2_PRO_5000001211_2/Orion/seds/' ;new SEDS have a problem
; ;dir_seds='/Users/ilyeschoubani/Downloads/SC110012_V2p2_PRO_5000001211/Orion/seds/'
dir_seds = '/Users/ilyeschoubani/Downloads/SC110012_V2p2_PRO_5000001210/Orion/seds/'
;with the mask - new batch still has the error
;file_sed=dir_seds+'extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'_'+msk+'_'+co_corr+'.xcat'
;without the mask - new batch does not have the error
file_sed=dir_seds+'extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'_'+co_corr+'.xcat'
;;THIS WAS PRIOR;file_sed=dir_seds+'extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'.xcat'
data_struct=read_xcat(file_sed,/silent) ;This is the new general structure
;stop
;
;===reading and applying the HCD:
openr,unit,file_sed,/get_lun
str=''
WHILE not eof(unit) DO BEGIN
readf,unit,str
if strtrim(strmid(str,13,3),2) EQ '857' then *!dustem_HCD=float(strtrim(strmid(str,21),2))/5.25*1.00E+26
ENDWHILE
close,unit
free_lun,unit
;===================
;stop
;generating fake errors for SED and Q/USED - here 10%
;stop
data_struct.sigmaII = (double((data_struct.StokesI*0.1)^2)) > 1d-20 ;this is false. Squaring it does not make sense. The square root is the one that should be taken
;NOT THE REAL REASON NOW;because of a bigger error in the variance computation in the xcat file. CHIPASS has an error value for StokesQ and StokesU but not for their associated flux.
;CHIPASS filter isn't well defined
;I looked it up and found a bandpass of 64 MHz (in 1024 channels).
data_struct(-1).sigmauu = la_undef()
data_struct(-1).sigmaqq = la_undef()
data_struct(-1).sigmaii = la_undef()
data_struct(-1).stokesu = la_undef()
data_struct(-1).stokesq = la_undef()
data_struct(-1).stokesi = la_undef()
data_struct(-1).psi = la_undef()
data_struct(-1).sigma_psi = la_undef()
; data_struct(-1).stokesi = la_undef()
; tt=where(data_struct.sigmauu EQ -39.8637)
; ;stop
;
; data_struct[tt].sigmauu = la_undef()
; data_struct[tt].sigmaqq = la_undef()
;data_struct.sigma_smallp = (double((data_struct.smallp*0.1)^2)) > 1d-20
idsigq=where(data_struct.sigmaQQ NE la_undef())
idsigu=where(data_struct.sigmaUU NE la_undef())
;stop
data_struct[idsigq].sigmaQQ = (double((data_struct[idsigq].StokesQ*0.1)^2)) > 1d-20
data_struct[idsigu].sigmaUU = (double((data_struct[idsigu].StokesU*0.1)^2)) > 1d-20
;#######FOR THE FITTING OF Q AND U
;locating the negative values in stokesQ data so that filters are omitted from the fit
indq=where(data_struct.stokesQ gt 0 and data_struct.stokesQ ne la_undef())
rmvq=(data_struct(indq).filter)
;locating the negative values in StokesU data so that filters are omitted form the fit
indu=where(data_struct.stokesU lt 0 and data_struct.stokesU ne la_undef())
rmvu=(data_struct(indu).filter)
rmvpsi=[rmvq,rmvu]
;removing LFI2 manually
remove,1,rmvpsi
;stop
;BECAUSE SPASS band has no errors:
;NOW IT DOES
; (data_struct(-2).sigma_largeP)=(data_struct(-3).sigma_largeP)
; (data_struct(-2).sigma_smallp)=(data_struct(-3).sigma_smallp)
; (data_struct(-2).sigmaqq)=(data_struct(-3).sigmaqq)
; (data_struct(-2).sigmauu)=(data_struct(-3).sigmauu)
; (data_struct(-2).sigma_psi)=(data_struct(-3).sigma_psi)
;stop
fit_struct=data_struct
;stop
;uncomment this block afterwards - because you are now just fitting the sed.
; fit_struct=dustem_mask_data(fit_struct,omit=list(rmvq,rmvu,rmvpsi,rmvpsi,rmvpsi),data_set=['QSED','USED','PSI','POLSED','POLFRAC']) ;works
; rmvq=(data_struct(-2).filter)
;
;
; fit_struct=dustem_mask_data(fit_struct,omit=list(rmvq,rmvq,rmvq,rmvq,rmvq),data_set=['QSED','USED','PSI','POLSED','POLFRAC'])
;
;
; fit_struct=dustem_mask_data(fit_struct,omit=rmvq,data_set=['QSED'])
; fit_struct=dustem_mask_data(fit_struct,omit=rmvu,data_set=['USED'])
;stop
;since we're using q and u I will have to null polsed and polfrac. Otherwise set_data will return an error
;fit_struct.largeP=la_undef()
;fit_struct.smallp=la_undef()
;doing the same for the errors as dustem_check will notice some erros have undefined associated data points.
;fit_struct.sigma_largeP=la_undef()
;fit_struct.sigma_smallp=la_undef()
;because I am now fitting only the sed
; fit_struct.sigmaQQ = la_undef()
; fit_struct.sigmaUU = la_undef()
;stop
;show structure
show_st=data_struct
;stop
; show_st.largeP=la_undef()
; show_st.sigma_largeP=la_undef()
; show_st.stokesQ=la_undef()
; show_st.sigmaQQ=la_undef()
; show_st.stokesU=la_undef()
; show_st.sigmaUU=la_undef()
;#################################################
;params_TH=[5.44E-05,5.60E-04,1.22E-04,3.0,4.13E+06,2.25E+00];,0.30,53]
params_G17C=[8.85E-06,1.25E-05,5.30E-03,100.0,3.00E3,6.43,1.00E+07,2.73E+00,58]
;params_DBP=[1.06E-05,7.00E-04,2.32E-03,6.43,1.00E+07,2.57E+00]
;=== Set which parameters you want to fit
pd = [ $
;'(*!dustem_params).gas.G0', $ ;G0
'(*!dustem_params).grains(0).mdust_o_mh', $
'(*!dustem_params).grains(1).mdust_o_mh', $
'(*!dustem_params).grains(2).mdust_o_mh', $
;'(*!dustem_params).grains(3).mdust_o_mh', $
;'(*!dustem_params).grains(4).mdust_o_mh', $
'(*!dustem_params).gas.Tgas',$
'(*!dustem_params).gas.nh',$
'dustem_plugin_stellar_population_O6V3',$
;'dustem_plugin_mbbdy_isrf_1',$
;'dustem_plugin_mbbdy_isrf_2', $
;'dustem_plugin_mbbdy_1',$
;'dustem_plugin_mbbdy_2', $
;'dustem_plugin_mbbdy_3', $
;'dustem_plugin_mbbdy_isrf_3', $
;'dustem_plugin_mbbdy_isrf_4', $
;'dustem_plugin_synchrotron_1', $
;'dustem_plugin_synchrotron_2', $
;'dustem_plugin_synchrotron_3', $
;'dustem_plugin_synchrotron_4', $
'dustem_plugin_freefree_1', $ ;4--->6
'dustem_plugin_freefree_2', $ ;5
;'dustem_plugin_freefree_3', $
;'dustem_plugin_modify_dust_pol_1',$
'dustem_plugin_modify_dust_pol_2'$
;'dustem_plugin_mdfy_spinning_pol_1'$
;'dustem_plugin_continuum_2', $
;'dustem_plugin_continuum_1' $
;'dustem_plugin_mdfy_spinning_pol_1', $
;'dustem_plugin_mdfy_spinning_pol_2' $
]
iv=params_G17C
Npar=n_elements(pd)
ulimed=replicate(0,Npar);[0]
ulims=replicate(0.0,Npar)
ulimed[6]=1
;uncomment
;ulimed[6]=1
;ulimed[9]=1
ulims[6]=1.00E+07
;ulims[8]=5.0E-03
;uncomment
;ulims[6]=9.00E-01;becomes 9
;stop
llimed=replicate(1,Npar);[1]
;llims=[1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20];,1.00e-20,1.00e-20]
llims=replicate(5.3e-6,Npar);[10]
;llims[0]=1.0e-9
;llims[3]=1.0
llims[6]=1.00E+04
llims[7]=1.00E-04
dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims;,tied=tied
dustem_init_plugins,pd
; fpd = [ $
; '(*!dustem_params).grains(1).mdust_o_mh'$
; ;'(*!dustem_params).grains(1).mdust_o_mh'$
; ;'(*!dustem_params).grains(4).mdust_o_mh'$
; ]
; fiv = [1.00E-20]
;
; dustem_init_fixed_params, fpd, fiv
;=== SET THE OBSERVATION STRUCTURE
ind=where(data_struct.wave GE 0.1)
st_ie=dustem_set_data(st_sed_fit=fit_struct[ind],st_sed_show=show_st[ind],rchi2_weight=rchi2_weight,f_HI=f_HI) ; comment if second line is to be used
;stop
;=== RUN fit
tol=1.e-10
xtol=1.e-10
Nitermax=3;0; just an initialization
xr=[1.,5e5]
yr=[5e-8,1.00e6]
loadct,13
;!y.range=[1e-8,100] ;This is to adjust plot range from outside the routine
yyy="restrng=textoidl('[MJy/sr]_{"+strmid(smooth_str,9)+'|'+plckprd+"}')"
yyy=execute(yyy)
title=textoidl('OrionA I_{240\mum} > ')+strtrim(lev[indxlev],2)+restrng
sttr=' I_{\nu} (MJy/sr) for N_{H}='+strmid(string(*!dustem_HCD,format='(1E10.2)'),2)+'H/cm^2'
ytit=textoidl(sttr)
;xtit=''
t1=systime(0,/sec)
res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,xtit=xtit,ytit=ytit,title=title,_extra=_extra);
t2=systime(0,/sec)
;help,(*!dustem_fit),/STRUCTURE
;=== SAVE FIT RESULTS
;NB: THIS LINE IS WRONG BECAUSE I AM NOT USING THE ANGLE CORRELATION CRITERIA
sedid='extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'_masked_'+polangmatch+'_'+msk
file='/tmp/DUSTEM_lastfit_orionA_'+sedid+'.sav'
dustem_save_system_variables,file
message,'Saved '+file,/continue
;======================================
;====You can exit IDL here and re-enter - THIS PART APPEARS TO BE FAULTY. NEEDS INVESTIGATION
;======================================
dustem_restore_system_variables,file
;=== Plot best fit
win=1 ; who wrote this?
window,win & win=win+1
;=== recover best fit values
res=*(*!dustem_fit).current_param_values
dustemwrap_plot,res,stp,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=title
print, 'DONE!!! IT WORKED'
; chi2=(*!dustem_fit).chi2
; rchi2=(*!dustem_fit).rchi2
; errors=*(*!dustem_fit).current_param_errors
; errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
; ;=== Plot best fit
; tit='DUSTEM polarization dust Example (Saved)'
; ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
; xtit=textoidl('\lambda (\mum)')
;
; loadct,13
;
; ps_file=!dustem_wrap_soft_dir+'Docs/'+'Last_dustem_fit.ps'
; dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=title,yr=yr,xr=xr,/ysty,/xsty,res=res,chi2=chi2,rchi2=rchi2,/xlog,/ylog,/pol,ps=ps_file,png=png,fpol=fpol
print,'dustem_mpfit_sed executed in ',t2-t1,' sec'
the_end:
END