dustem_extract_sed_example.pro
4.63 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
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 using the continuum plugin
;=== INFORMATION TO RUN THE FITS
tol=1.e-16 ;fit tolerence
use_Nitermax=20
;=== 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_plugin_continuum_1', $ ;Temperature of a BB
'dustem_plugin_continuum_2'] ;Peak amplitude of a BB
iv = [15., 1.e6] ; initial guess
Npar=n_elements(pd)
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
llims=replicate(1.e-15,Npar)
;; DustEM requires a physical dust model. Here we have used DBP90, but
;; fix all grain components to ~zero.
fpd = [ $
'(*!dustem_params).G0', $ ;G0
'(*!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.0, 1.e-12,1.e-12,1.e-12]
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