PRO dustem_make_polarization_sed_example,model=model $
                                         ,filters=filters $
                                         ,outfile=outfile $
                                         ,help=help $
                                         ,verbose=verbose 

;+
; 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 DustEMWrap
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:
;    dustem_make_polarization_sed_example[,model=][,filters=][,/help]
;
; INPUTS:
;    None
;
; OPTIONAL INPUT PARAMETERS:
;    filters = instrument bands for which SED is created. Default is
;    ['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
;
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUT PARAMETERS:
;    outfile = path+filename for .xcat output of the SED. Default is './polarization_sed.xcat'
;
; ACCEPTED KEY-WORDS:
;    model     = if set, name of the DustEM dust model to use (default='G17_MODELD')
;    help      = If set print this help
;    verbose      = If set, subroutines will be more chatty
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
;
; PROCEDURE:
;    None
;
; EXAMPLES
;    dustem_make_polarization_sed_example
;    dustem_make_polarization_sed_example,model='G17_MODELA',outfile='mysed.xcat'
;
; MODIFICATION HISTORY:
;    Written by JPB June-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_make_polarization_sed_example'
  goto,the_end
ENDIF

dustem_define_la_common

IF keyword_set(model) THEN BEGIN
  use_model=strupcase(model)
ENDIF ELSE BEGIN
  use_model='G17_MODELD'    ;Default is one of the Guillet et al (2017) models since they treat dust polarisation properties
ENDELSE

IF keyword_set(filters) THEN BEGIN
  use_filters=strupcase(filters)
ENDIF ELSE BEGIN
  use_filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
ENDELSE

IF keyword_set(outfile) THEN BEGIN
  use_outfile=outfile
ENDIF ELSE BEGIN
  use_outfile='./polarization_sed.xcat'
ENDELSE

known_mdls=['MC10','DBP90','DL01','WD01_RV5P5B','DL07','J13','G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD'] 
pol_mdls=['G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD'] 
test_model = where(known_mdls eq use_model,ct)
if ct eq 0 then begin
   message,'ISM dust model '+use_model+' unknown',/continue
   message,'Known models are MC10,DBP90,DL01,WD01_RV5P5B,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
   stop
end
test_model = where(pol_mdls eq use_model,ct)
if ct eq 0 then begin
   message,'The only models with polarisation are G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
   stop
end


use_verbose=0
use_polarization=1
if keyword_set(verbose) then use_verbose=1

dustem_init,mode=use_model,polarization=use_polarization
!dustem_verbose=use_verbose
!dustem_nocatch=1

;=== initialize filters for the SED
Nfilt=n_elements(use_filters)
sed=dustem_initialize_sed(Nfilt)
sed.filter=use_filters
sed.wave=dustem_filter2wav(use_filters)
sed.instru=dustem_filter2instru(use_filters)

;=== initialize 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.

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

;=== Set which model parameters to use for the SED
pd = [ $
             '(*!dustem_params).G0', $                   ;G0
             '(*!dustem_params).grains(0).mdust_o_mh',$  ;PAH0_MC10 mass fraction
             '(*!dustem_params).grains(1).mdust_o_mh',$  ;amCBE_0.3333x mass fraction   
             '(*!dustem_params).grains(2).mdust_o_mh', $    ;aSil2001BE6pctG_0.4x mass fraction
             'dustem_plugin_modify_dust_pol_2',  $     ;Dust 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
             ]                   
pd_true=[1.,7.8000E-04,7.8000E-04,7.8000E-04, 48.,3., 1.e-2, 0.3, 45.]

;== SET INITIAL VALUES AND LIMITS 
;== AND ACTIVATE ANY PLUGINS
dustem_init_params,use_model,pd,pd_true,pol=use_polarization

;=== compute the SED (I then QU) using the model and specified
;=== true_vals of parameters
dustem_Ised=dustem_compute_sed(pd_true,st=ssti)   ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
toto=dustem_compute_stokes(pd_true,st=sstqu) ;this procedure also allows for the extraction of the spectra
dustem_qsed = toto[0]
dustem_used = toto[1]


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

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

;=== set covariances to 0.
sed.sigmaIQ=0.
sed.sigmaIU=0.
sed.sigmaQU=0.

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

;======== save the SED to a file
write_xcat,sed,use_outfile
message,'Wrote '+use_outfile,/continue

;======== move the file to the Data/EXAMPLE_OBSDATA/ subdirectory 
;filename_final=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/my_example_SED_###.xcat'
;str='cp '+use_outfile+' '+filename_final
;message,'Do '+str+' to make change permanent',/continue

the_end:

END