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

;+
; 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:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:

;    outfile = path+filename for .xcat output of the SED

; 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 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_MODELD'

; 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


if keyword_set(verbose) then use_verbose=1

dustem_define_la_common

;dustem_init,mode=use_model,grain_keywords=['logn-chrg-spin','?','plaw-pol'],polarization=1
dustem_init,mode=use_model,polarization=1

!dustem_verbose=1
!dustem_show_plot=1
!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

;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_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, 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_params,use_model,pd,pd_true,pol=1
;dustem_init_parinfo,pd,pd_true,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
;dustem_init_plugins,pd

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]

;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 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