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