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 CESR
;    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
;dustem_init,mode='COMPIEGNE_ETAL2010',/pol

;use_mode='G17_MODELA'

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

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 tags that will be internally used 
dustem_set_data,sed

help,!dustem_data

;=== Set which model parameters to use for the SED
;=== This does not work, for some reason ....
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_1', $    ;This will set the polarization fraction
             'dustem_plugin_modify_dust_pol_2'  $     ;This will set the polarization angle
              ]                   
p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04,0.2,45.]

;=== 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
             ]                   
p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04]

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_synchrotron_2', $   ;Synchrotron amplitude at 10 mm
             'dustem_plugin_synchrotron_3', $    ;Synchrotron polarization angle
             'dustem_plugin_synchrotron_4'  $    ;Synchrotron polarization fraction
             ]                   
p_truth=[1.,7.8000E-04,7.8000E-04,7.8000E-04, 1., 45., 0.3]

;===== This example adds synchrotron with a polarization angle of 45 degree
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_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,3., 1.e-2, 0.3, 45.]

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.]


;iv=p_truth+[0.,8.e-4,8.e-4,8.e-4]      ;shifted from solution
iv=p_truth

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

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

;Here fortran fails when pluggins are activated and/or type-keywods exist
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
;toto=dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em,out_st=out_st,_extra=extra)

print,dustem_Qsed,dustem_Used

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


;goto,no_change
;=modifies the Fortran outputs Q,U SED values to reflect a fixed polarization fraction and angle
;normaly, this could also be doen using the plugin:
;toto=dustem_plugin_modify_dust_pol(ssti, key=key, val=val, scope=scope, paramtag=paramtag, help=help)
;frac_used=0.01  ;1%
;psi_used=10.   ;10 degrees

;smallp=dustem_Ised*0.+frac_used
;psi=dustem_Ised*0.+psi_used
;polar_ippsi2iqu,dustem_Ised,Q,U,smallp,psi

;sed.stokesI=dustem_Ised
;sed.stokesQ=Q
;sed.stokesU=U

;no_change:

;=== set arbitrary uncertainties
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)

;cgplot,sed.wave,sed.stokesI,psym=-3,/xlog,/ylog

;stop
;======== 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
;spawn,str

;message,'Wrote '+filename_final,/continue

the_end:

END