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