dustem_extract_sed_example.pro
1.49 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
PRO test_dustem_sed_extractor
define_la_common
dustem_init
Nmaps=2 ;This is total intensity and variances
filters=['IRAS3','IRAS4','PACS1','PACS2','PACS3','PILOT1','SPIRE1','SPIRE2','SPIRE3','HFI3','HFI4','HFI5','HFI6']
Ndata_sets=n_elements(filters)
maps=fltarr(100,100,Nmaps,Ndata_sets)
maps[*]=la_undef()
mask=maps[*,*,0,0]
mask[5:35,5:35]=1
mask[10:30,10:30]=mask[10:30,10:30]+1
index_on=where(mask EQ 2,count_on)
index_off=where(mask EQ 1,count_off)
indI=0
indII=1
maps_order=[indI,indII]
temp_on=40. ;K
temp_off=20.
FOR i=0L,Ndata_sets-1 DO BEGIN
;fill in intensity data
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
;fill in intensity uncertainties
maps[*,*,indII,i]=la_power(la_mul(maps[*,*,indI,i],0.3),2.)
ENDFOR
sed=dustem_sed_extractor(maps,index_on,filters,maps_order=maps_order,/total_intensity_only)
nsigma=10
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)
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'
stop
the_end:
END