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 (Stokes
;    IQU) using one of the physical dust models in 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 = specifies the interstellar dust mixture used by
;            DustEM. See userguide or dustem_test_model_exists.pro
;            for more details about available models in current release.

;    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

exists=dustem_test_model_exists(use_model,/pol)
if exists ne 1 then $
   message,'Unknown/unpolarised dust model: '+use_model

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

;== 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_fit,m_show=m_show,x_fit=ext,x_show=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

;======== example of how to move the file to the Data/EXAMPLE_OBSDATA/ subdirectory 
;filename_final=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/my_example_EXT_G17MDL_D.xcat'

;str='cp '+use_outfile+' '+filename_final
;message,'Do '+str+' to make change permanent',/continue

the_end:

message,'Finished dustem_make_polarization_ext_example',/info

END