PRO dustem_make_polarization_ext_example,model=model $
                                         ,wave=wave $
                                         ,next=next $
                                         ,outfile=outfile $
                                         ,help=help $
                                         ,verbose=verbose


;+
; NAME:
;    dustem_make_polarization_ext_example
;
; PURPOSE:
;    This is an example of how to generate an extinction curve 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_ext_example[,model=][,filters=][,/help]
;
; INPUTS:
;    wave = wavelength range specified as [minwave,maxwave] in microns over which
;           to construct the extinction curve. The resulting curve
;           will be defined on the vector :
;           wave_vector=10^(alog10(wave[0])+alog10(wave[1])*findgen(Next)/float(Next))
;           default is [0.01,50]
;    Next = number of measurements in the extinction curve. Default is 100.
;
; OPTIONAL INPUT PARAMETERS:
;    None
;
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUT PARAMETERS:
;    outfile = path+filename for .xcat output of the EXT
;
; 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_ext_example
;    dustem_make_polarization_ext_example,model='G17_MODELC',outfile='my_polarized_ext.xcat'
;
; MODIFICATION HISTORY:
;    Written by AH Nov-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_ext_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(outfile) THEN BEGIN
  use_outfile=outfile
ENDIF ELSE BEGIN
  use_outfile='./polarization_ext.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
use_next=100 ; number of extinction measurements
use_wave=10^(alog10(0.01)+alog10(50)*findgen(use_next)/float(use_next)) ; wavelengths at which extinction is defined
if keyword_set(verbose) then use_verbose=1
if keyword_set(next) then use_next=next
if keyword_set(wave) then use_wave=10^(alog10(wave[0])+alog10(wave[1])*findgen(use_next)/float(use_next))

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

;=== initialize wavelengths for the EXTINCTION
ext=dustem_initialize_ext(use_next)
ext.instru='EXTINCTION'
ext.filter='SPECTRUM'
ext.wave=use_wave

;=== initialize IQU and associated errors to avoid problems when checking EXT in dustem_set_data.pro
ext[*].EXT_I=1.
ext.EXT_Q=ext.EXT_I/100.
ext.EXT_U=ext.EXT_I/100.
ext.SIGEXTII=ext.EXT_I/1000.
ext.SIGEXTQQ=ext.EXT_I/1000.
ext.SIGEXTUU=ext.EXT_I/1000.
ext.SIGEXTIQ=ext.EXT_I/1000.
ext.SIGEXTIU=ext.EXT_I/1000.
ext.SIGEXTQU=ext.EXT_I/1000.

;=== Set which model parameters to use for the EXT, this will depend
;=== on the selected dust model. For model parameters not specified
;=== here, the default Dustem fortran values will be used
pd = [ $
     '(*!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_polx_2' $                 ;Dust polarization angle in extinction
     ]                   
true_vals=[7.8000E-04,7.8000E-04,7.8000E-04, 52.]


Npar=n_elements(pd)  
ulimed=replicate(0,Npar)
llimed=replicate(1,Npar) & llimed[3]=0
llims=replicate(0,Npar)

;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT 
;== AND ACTIVATE ANY PLUGINS
;== FOR REASONS, WE NEED TO DO THIS DUSTEM_SET_DATA CALL
dustem_init_params,use_model,pd,true_vals,pol=use_polarization,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims


;=== initialize the !dustem_data structure that will be used internally
;== note that m_fit and m_show do not exist -- this is deliberate
;== since we are not treating emission data
dustem_set_data,m_fit,m_show,ext,ext

;=== compute the extinction curve (I then QU) using the model and specified
;=== true_vals of parameters
dustem_iext=dustem_compute_ext(true_vals,st=ssti) ;sst is not set on purpose, for dustem_compute_ext to compute the ext using fortran code
toto=dustem_compute_stokext(true_vals,st=sstqu) ;this procedure also allows for the extraction of the spectra
dustem_qext = toto[0]
dustem_uext = toto[1]

ext.EXT_I=dustem_Iext
ext.EXT_Q=dustem_Qext
ext.EXT_U=dustem_Uext

;=== set uncertainties to I,Q,U
ext.SIGEXTII=abs(ext.EXT_I*0.000001)
ext.SIGEXTQQ=abs(ext.EXT_Q*0.000001)
ext.SIGEXTUU=abs(ext.EXT_U*0.000001)

;=== set covariances to 0.
ext.SIGEXTIQ=0.
ext.SIGEXTIU=0.
ext.SIGEXTQU=0.

;==== fill in dependent columns of the EXT.
ext=dustem_fill_ext_dependent_columns(ext)

;======== save the EXT to a file
write_xcat,ext,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_EXT_###.xcat'
;str='cp '+use_outfile+' '+filename_final
;message,'Do '+str+' to make change permanent',/continue

the_end:

END