PRO dustem_make_polarization_sed_example,model=model,help=help

;+
; NAME:
;    dustem_make_polarization_sed_example
; PURPOSE:
;    This is an example of how to generate an SED using dustemwrap.
;    It is meant to be an example to follow when writing your own
;    programs using the dustem IDL/GDL wrapper.
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_make_polarization_sed_example[,model=][,/help]
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    model     = if set, name of the Dustem dust model to use (default='G17_MODELC')
;    help      = If set print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    dustem_make_polarization_sed_example
;    dustem_make_polarization_sed_example,model='G17_MODELD'
; MODIFICATION HISTORY:
;    Written by J.P. Bernard June 29th 2011
;    see evolution details on the dustem cvs maintained at IRAP
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_fit_sed_polsed_readme'
  goto,the_end
ENDIF

dustem_define_la_common

;=== initialise dustem
use_mode='G17_MODELC'
dustem_init,mode=use_mode,kwords=['logn-chrg-spin','plaw-pol','plaw-pol'],/pol

!dustem_verbose=1
!dustem_show_plot=1
!dustem_nocatch=1

;=== define filters for the SED
filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
Nfilt=n_elements(filters)
sed=dustem_initialize_sed(Nfilt)
sed.filter=filters
sed.wave=dustem_filter2wav(filters)
sed.instru=dustem_filter2instru(filters)

;=== initializing IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro
sed[*].StokesI=1.
sed.StokesQ=sed.StokesI/100.
sed.StokesU=sed.StokesI/100.
sed.SigmaII=sed.StokesI/100.
sed.SigmaQQ=sed.StokesI/100.
sed.SigmaUU=sed.StokesI/100.
sed.SigmaIQ=sed.StokesI/100.
sed.SigmaIU=sed.StokesI/100.
sed.SigmaQU=sed.StokesI/100.

;=== initializing the !dustem_data structure that will be internally used 
dustem_set_data,sed

;help,!dustem_data

;=== Set which model parameters to use for the SED
pd = [ $
             '(*!dustem_params).G0', $                   ;G0
             '(*!dustem_params).grains(0).mdust_o_mh',$  ;PAH0 mass fraction
             '(*!dustem_params).grains(1).mdust_o_mh',$  ;PAH0 mass fraction   
             '(*!dustem_params).grains(2).mdust_o_mh', $    ;PAH1 mass fraction
             'dustem_plugin_modify_dust_pol_2',  $     ;This will set the polarization angle
             'dustem_plugin_synchrotron_1', $   ;Synchrotron spectra index
             'dustem_plugin_synchrotron_2', $   ;Synchrotron amplitude at 10 mm
             'dustem_plugin_synchrotron_3', $    ;Synchrotron polarization fraction
             'dustem_plugin_synchrotron_4'  $    ;Synchrotron polarization angle
             ]                   
p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,10.,3., 1.e-2, 0.3, 45.]

Npar=n_elements(pd)  
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar)
llims=replicate(0,Npar)

dustem_init_parinfo,pd,p_truth,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
dustem_init_plugins,pd

dustem_Ised=dustem_compute_sed(p_truth,ssti)   ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
toto=dustem_compute_stokes(p_truth,sstqu,dustem_Qsed,dustem_Used) ;this procedure also allows for the extraction of the spectra

;print,dustem_Qsed,dustem_Used

sed.stokesI=dustem_Ised
sed.stokesQ=dustem_Qsed
sed.stokesU=dustem_Used

;=== set arbitrary uncertainties to I,Q,U
sed.sigmaII=abs(sed.StokesI*0.000001)
sed.sigmaQQ=abs(sed.StokesQ*0.000001)
sed.sigmaUU=abs(sed.StokesU*0.000001)

sed.sigmaIQ=0.
sed.sigmaIU=0.
sed.sigmaQU=0.

;print,sed.stokesI,sed.stokesQ,sed.stokesU
;print,sed.sigmaII,sed.sigmaQQ,sed.sigmaUU

;==== fill in dependent columns of the SED.
sed=dustem_fill_sed_dependent_columns(sed)

;======== save the SED
filename='/tmp/fake_polarization_sed.xcat'
write_xcat,sed,filename

message,'Wrote '+filename,/continue

filename_final=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/example_sed3.xcat'
str='cp '+filename+' '+filename_final
message,'Do '+str+' to make change permanent',/continue

the_end:

END