dustem_extract_sed_example.pro
4.62 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
PRO dustem_extract_sed_example
;+
; NAME:
; dustem_extract_sed_example
;
; PURPOSE:
; This routine contains a simple example of how to extract an SED to use with DustEMWrap.
;
; CATEGORY:
; DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:
; dustem_extract_sed_example
;
; INPUTS:
; None
;
; OPTIONAL INPUT PARAMETERS:
; None
;
; OUTPUTS:
; None
;
; OPTIONAL OUTPUT PARAMETERS:
;
; ACCEPTED KEY-WORDS:
; None
;
; COMMON BLOCKS:
; None
;
; SIDE EFFECTS:
; None
;
; RESTRICTIONS:
; The DustEM fortran code must be installed
; The DustEMWrap IDL code must be installed
;
; PROCEDURES AND SUBROUTINES USED:
;
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by JPB Jan 2022
; Evolution details on the DustEMWrap gitlab.
; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
;-
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_extract_sed_example'
goto,the_end
END
define_la_common
use_model='DBP90'
use_polarization=0
dustem_init,model=use_model,polarization=use_polarization
Nmaptypes=2 ; number of map types in the structure. Here 2 (Stokes I, variance in I)
filters=['IRAS3','IRAS4','PACS1','PACS2','PACS3','PILOT1','SPIRE1','SPIRE2','SPIRE3','HFI3','HFI4','HFI5','HFI6']
Ndata_sets=n_elements(filters)
maps=fltarr(100,100,Nmaptypes,Ndata_sets)
maps[*]=la_undef()
;; Define the region of interest (ON) and OFF region
mask=maps[*,*,0,0]
mask[5:35,5:35]=1 ; OFF will be a guard band of 5 pixel around the ON defined below
mask[10:30,10:30]=mask[10:30,10:30]+1 ; ON source
index_on=where(mask EQ 2,count_on)
index_off=where(mask EQ 1,count_off)
indI=0
indII=1
maps_order=[indI,indII]
;; Loop through the maps corresponding to different filters
;; Generate some fake data using a BB function
;; ON pixels will be set to blackbody emission at 40K
;; OFF pixels will be set to blackbody emission at 20K
;; Loop through the maps corresponding to different filters
temp_on=40. ;Kelvin
temp_off=20.
FOR i=0L,Ndata_sets-1 DO BEGIN
wave=dustem_filter2wav(filters[i])
this_map0=maps[*,*,indI,i]
this_map=this_map0
this_map[index_on]=dustem_planck_function(temp_on,wave)
this_map[index_off]=dustem_planck_function(temp_off,wave)
maps[*,*,indI,i]=this_map
;; assign a variance value to each valid pixel
maps[*,*,indII,i]=la_power(la_mul(maps[*,*,indI,i],0.3),2.)
ENDFOR
;; extract the default type of Stokes I only SED for the ON pixels (method = 'MEAN')
sed=dustem_sed_extractor(maps,index_on,filters,maps_order=maps_order,/total_intensity_only)
;; show SED with 3-sigma error bars
nsigma=3
cgplot,sed.wave,sed.stokesI,/xlog,/ylog,xtit='Wavelength [mic]',ytit='SED [MJy/sr]'
cgoplot,sed.wave,sed.stokesI,psym='Filled Circle',color='red'
cgerrplot,sed.wave,sed.stokesI-nsigma*sqrt(sed.sigmaII),sed.stokesI+nsigma*sqrt(sed.sigmaII)
;; extract a MEAN_ONOFF Stokes I only SED for the ON pixels (method = 'MEAN_ONOFF')
sed2=dustem_sed_extractor(maps,index_on,filters,maps_order=maps_order,/total_intensity_only,index_off=index_off,method='MEAN_ONOFF')
cgoplot,sed2.wave,sed2.stokesI
cgoplot,sed2.wave,sed2.stokesI,psym='Filled Circle',color='blue'
;; Try a simple BB fit of a dust model to SED
;=== INFORMATION TO RUN THE FITS
tol=1.e-16 ;fit tolerence
use_Nitermax=10
;=== INFORMATION TO MAKE THE PLOTS
yr=[1.00e1,8.00E8] ; y-axis limits
xr=[5.00E0,2.00e4] ; x-axis limits
tit='DUSTEM SED EXTRACT EXAMPLE' ; plot title
ytit=textoidl('I_\nu (MJy/sr)') ; y-axis title
xtit=textoidl('\lambda (\mum)') ; x-axis title
pd = [ $
'(*!dustem_params).G0', $ ;G0
'dustem_plugin_continuum_1', $ ;Temperature of a BB
'dustem_plugin_continuum_2'] ;Peak amplitude of a BB
iv = [1., 10., 1.e6] ; initial guess
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
llims=replicate(1.e-15,Npar)
fpd = [ $
'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH
'(*!dustem_params).grains(1).mdust_o_mh',$ ;VSG
'(*!dustem_params).grains(2).mdust_o_mh'] ;BG
fiv = [1.e-12,1.e-12,1.e-12]
ind=where(sed.sigmaII LT (0.2*sed.StokesI)^2,count)
IF count NE 0 THEN sed[ind].sigmaII=(0.2*sed[ind].StokesI)^2
dustem_set_data,m_fit=sed,m_show=sed
dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims
res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit $
,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot)
the_end:
END