PRO dustem_init,dir=dir,wraptest=wraptest,plot_it=plot_it,model=model,help=help,noerase=noerase,st_model=st_model,pol=pol,verbose=verbose,kwords=kwords

;+
; NAME:
;   dustem_init
; CALLING SEQUENCE:
;   dustem_init
; PURPOSE:
;   Defines system variables to run the DUSTEM code
; INPUTS:

;    model    = Selects one of the dust mixture used by dustem

;              'COMPIEGNE_ETAL2010' from Compiegne et al 2010 (default)
;              'DBP90' from Desert et al 1990
;              'DL01' from Draine & Li 2001
;              'DL07' from Draine & Li 2007
;   dir      = overrides default directory where to read data
;   wraptest = test running dustem f90 through the wrapper
;   plot_it  = plots run result (only if wraptest)
;   help     = If set print this help
; OPTIONAL INPUT:
;   None
; OUTPUTS:
;   None
; PROCEDURE AND SUBROUTINE USED
;    1/ The folowing system variables are initialized to null values
;	or empty pointers:
;       !dustem_fit
;       !dustem_data
;       !dustem_inputs
;       !dustem_params
;       !run_ionfrac
;       !dustem_idl_continuum
;       !dustem_filters
;       !dustem_verbose
;       !dustem_show_plot
;    2/ optionally runs the fortran code (if wraptest)
;       This step calls the Fortran code to run a standard model
;    3/ !dustem_filters is initialized calling dustem_filter_init.pro.
; RESTRICTIONS:
;    The dustem fortran code must be installed
;    The dustem idl wrapper must be installed
; EXAMPLES
;    dustem_init,/wrap,/plot_it
; MODIFICATION HISTORY:
;    See cvs logs
;    Written by JPB,NF,DP Jan-2007
;    see evolution details on the dustem cvs maintained at CESR
;-

;IC 
;I think this procedure should be throughly revised. (and also the rest of the code corrected)
;Pointers are used when the need does not present itself. 
;ie: defsysv, '!dustem_fit', ptr_new(dustem_fit_st), defsysv,'!dustem_inputs',ptr_new() or defsysv,'!dustem_params',ptr_new() which don't require initialization for instance

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_init'
  goto,the_end
ENDIF

dustem_define_la_common 

defsysv,'!dustem_verbose',0             ;1=verbose
IF keyword_set(verbose) THEN !dustem_verbose=1

defsysv,'!run_pol',0             ;0=no polar
IF keyword_set(pol) THEN !run_pol=1

;=== color correction system variables
defsysv,'!dustem_previous_cc',ptr_new()   ;store previous color correction values
defsysv,'!dustem_do_cc',1   ;indicates to dustem if color corrections are needed. Will be changed when fitting
defsysv,'!dustem_never_do_cc',0   ;set to 1 to never do color corrections (not recommanded, but faster).

;This is for polarization
defsysv, '!run_circ', 0.  
defsysv, '!run_anis', 0.  
defsysv, '!run_rrf', 0.    
defsysv, '!run_univ', 0.    

defsysv, '!run_lin', 0.

defsysv, '!dustem_redshift', 0.    ;in case the SED needs to be shifted for a non zero redshift

;This is just to make sure that plugin system variables are created even if no plugins are used
pd_bidon=['NONE']
dustem_init_plugins,pd_bidon

defsysv,'!dustem_nocatch',0   ;set to catch errors in sed fitting.

;stop

;==== define a stellar modyfier to the ISRF

;toto=dustem_define_isrf_stellar_modifyer_variable()

file=!dustem_wrap_soft_dir+'instrument_description.xcat'
st=read_xcat(file,/silent)
defsysv,'!dustem_instrument_description',st

defsysv,'!indef', exists=exists_indef
IF NOT exists_indef THEN defsysv,'!indef',-32768.

;define !dustem_fit
dustem_fit_st={data:ptr_new(),wavelength:ptr_new(), $
               param_descs:ptr_new(),param_init_values:ptr_new(),param_func:ptr_new(), $
               fixed_param_descs:ptr_new(),fixed_param_init_values:ptr_new(), $
               chi2:0.D0,rchi2:0.D0, $
               current_param_values:ptr_new(), current_param_errors:ptr_new() $
               }
defsysv, '!dustem_fit', ptr_new(dustem_fit_st)    ;Data to fit

if keyword_set(pol) then begin
	defsysv, '!dustem_data', {    $ ;Data to fit

                                sed:    ptr_new(),      $ ; 
                                ext:    ptr_new(),      $ ;

                                polext: ptr_new(),      $ ; FREE TO USE
                                polsed: ptr_new(),      $ ; FREE TO USE
                                polfrac: ptr_new(),     $ ; FREE TO USE  
                                psi_em: ptr_new(),      $ ; FREE TO USE
                                psi_ext: ptr_new(),     $ ; FREE TO USE

                                qsed: ptr_new(),        $ ; 
                                used: ptr_new(),        $ ; 
                                qext: ptr_new(),        $ ; 
                                uext: ptr_new()         $ ; 
                         }
    
         ;rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.,qsed:0.,used:0.,qext:0.,uext:0.} ;polsed and polfrac newly deprecated
         rchi2_weight = {sed: 0.,ext: 0.,polext:0.,polsed:0.,polfrac:0.,psi_em:0.,psi_ext:0.,qsed:0.,used:0.,qext:0.,uext:0.}
         
    defsysv,'!dustem_show', {    $ ;Data to fit

                                sed:    ptr_new(),      $

                                ext:    ptr_new(),      $

                                polext: ptr_new(),      $
                                polsed: ptr_new(),      $ 
                                polfrac: ptr_new(),     $
                                psi_em: ptr_new(),      $
                                psi_ext: ptr_new(),     $ 

                                qsed: ptr_new(),        $
                                used: ptr_new(),        $
                                qext: ptr_new(),        $
                                uext: ptr_new()         $

                         }     
         

endif else begin

	defsysv, '!dustem_data', {    $ ;Data to fit
                                sed:    ptr_new(),      $
                                ext:    ptr_new()       $
                         }
	rchi2_weight = {sed: 0.,ext: 0.}

	
	defsysv, '!dustem_show', {    $ ;Data to fit
                                sed:    ptr_new(),      $
                                ext:    ptr_new()       $
                         }

endelse

if !run_pol then begin
    tagnames=tag_names(!dustem_data)
    testpol = tagnames EQ 'POLSED' or tagnames EQ 'POLEXT' or tagnames EQ 'POLFRAC' or tagnames EQ 'POLSED' or tagnames EQ 'PSI_EM' or tagnames EQ 'PSI_EXT'
    ind_data=where(~testpol,ctestpol)
    if ctestpol ne 0 then tagnames =(tag_names(!dustem_data))[ind_data] 

endif else tagnames=tag_names(!dustem_data)

instr="defsysv, '!fit_rchi2_weight', {"
FOR i = 0, n_elements(tagnames)-1 DO BEGIN
  instr+=tagnames(i)+': rchi2_weight.'+tagnames(i)
  IF i NE n_elements(tagnames)-1 THEN instr+=','
ENDFOR
instr+='}'
toto=execute(instr)

;#########system variables necessary for cgwindow plotting##########     
defsysv, '!dustemcgwin_id', {    $ ;IDs of windows to plot                  
                            sed:    la_undef(),      $                          
                            ext:    la_undef()       $                 
                     }                                                      
;I am not sure about the layout of this system variable/structure   
                                                                                                
defsysv, '!dustemcgwin_ncmds', {    $ ;Keeping track of number of commands run.
                            sed:     {pl:0, nrm:0},      $                              
                            ext:     {pl:0, nrm:0},       $                         
                            polext:  {pl:0, nrm:0},      $                          
                            polsed:  {pl:0, nrm:0},      $                                      
                            polfrac: {pl:0, nrm:0},     $                       
                            psi_em:  {pl:0, nrm:0},      $                      
                            psi_ext: {pl:0, nrm:0},     $                           
                            qsed:    {pl:0, nrm:0},        $                    
                            used:    {pl:0, nrm:0},        $                        
                            qext:    {pl:0, nrm:0},        $                    
                            uext:    {pl:0, nrm:0}         $                        
                     }                                                                      
                                                                                                
defsysv, '!dustem_mlog', 0  ;necessary for the plotting of negative values on the axis in dstmwrp_exp.pro                                                                 
;###################################################################


defsysv, '!dustem_psi', 0.

;IDL> help,*!dustem_data,/str
;** Structure <2173b20>, 4 tags, length=448, data length=448, refs=1:
;   INSTRU_NAMES    STRING    Array[14]
;   FILT_NAMES      STRING    Array[14]
;   VALUES          FLOAT     Array[14]
;   SIGMA           FLOAT     Array[14]
defsysv,'!dustem_inputs',ptr_new()    ;dustem input files

;=== This is not a system variable anymore
;defsysv, '!dustem_ext', ptr_new()    ;Contains the values of extinction
;=== This is not a system variable anymore
;defsysv, '!dustem_specem', ptr_new() ;Contains the values of emission

defsysv, '!dustem_params', ptr_new() ;Contains the values of all Desert Model parameters (except extinction !)
;IDL> help,*!dustem_params,/str
;** Structure <2926e10>, 7 tags, length=13312, data length=13282, refs=1:
;   GRAINS          STRUCT    -> <Anonymous> Array[3]
;   GEMISSIV        STRUCT    -> <Anonymous> Array[300]
;   ISRF            STRUCT    -> <Anonymous> Array[64]
;   QABS            STRUCT    -> <Anonymous> Array[64]
;   CALOR           STRUCT    -> <Anonymous> Array[5]
;   SPECEM          STRUCT    -> <Anonymous> Array[1]
;   IONFRAC         STRUCT    -> <Anonymous> Array[10]

;defsysv, '!dustem_fit_sed',ptr_new()          ;Best fit SED

;The following is not used in the WEB3p8 version of the code
;Only in the VERSTRAETE version of the code.
;IF getenv('DUSTEM_WHICH') EQ 'VERSTRAETE' THEN BEGIN
  defsysv, '!run_ionfrac', 0.     ; 0:no call to dustem_create_ionfrac UNLESS it is called as a param
                                  ; 1:call to dustem_create_ionfrac with the MIX_xxx.DAT files
                                  ; 2:call to dustem_create_ionfrac with the SIZE_xxx.DAT files
;ENDIF

defsysv,'!dustem_filters',ptr_new()     ;filter information (filled by dustem_filter_init.pro)
;; IDL> help,*!dustem_filters,/str
;; ** Structure <1b95608>, 17 tags, length=8232, data length=8118, refs=1:
;;    IRAC            STRUCT    -> <Anonymous> Array[1]
;;    MIPS            STRUCT    -> <Anonymous> Array[1]
;;    MSX             STRUCT    -> <Anonymous> Array[1]
;;    IRAS            STRUCT    -> <Anonymous> Array[1]
;;    DIRBE           STRUCT    -> <Anonymous> Array[1]
;;    SPM             STRUCT    -> <Anonymous> Array[1]
;;    ISOCAM          STRUCT    -> <Anonymous> Array[1]
;;    ISOPHOTP        STRUCT    -> <Anonymous> Array[1]
;;    ISOPHOTC        STRUCT    -> <Anonymous> Array[1]
;;    EFIRAS          STRUCT    -> <Anonymous> Array[1]
;;    ARCHEOPS        STRUCT    -> <Anonymous> Array[1]
;;    HFI             STRUCT    -> <Anonymous> Array[1]
;;    LFI             STRUCT    -> <Anonymous> Array[1]
;;    WMAP            STRUCT    -> <Anonymous> Array[1]
;;    SPIRE           STRUCT    -> <Anonymous> Array[1]
;;    PACS            STRUCT    -> <Anonymous> Array[1]
;;    AKARI           STRUCT    -> <Anonymous> Array[1]

defsysv,'!dustem_verbose',0                                  ;1=verbose

defsysv,'!dustem_show_plot',1           ;0=show no plot

;apparently not needed anymore. Was removed.
;defsysv,'!ki',0.  ;apparently needed by mpfit, but what is it ?

;== This is for the default mode
dustem_inputs={GRAIN:'GRAIN.DAT',ISRF:'ISRF.DAT',ALIGN:'ALIGN.DAT'}
!dustem_inputs=ptr_new(dustem_inputs)

defsysv, '!dustem_HCD', ptr_new(1.00E20) ;I do not see the reason behind using a pointer in this case.

;===This is the Desert option

IF keyword_set(model) THEN BEGIN

defsysv, '!dustem_model', model

  CASE model OF

    'COMPIEGNE_ETAL2010':BEGIN
      (*!dustem_inputs).grain='GRAIN_MC10.DAT'
    END
    'DBP90':BEGIN
      (*!dustem_inputs).grain='GRAIN_DBP90.DAT'
    END
    'DL01':BEGIN
      (*!dustem_inputs).grain='GRAIN_DL01.DAT'
    END
    'DL07':BEGIN
      (*!dustem_inputs).grain='GRAIN_DL07.DAT'

    END
    'AJ13':BEGIN

      (*!dustem_inputs).grain='GRAIN_J13.DAT'

    END
    'G17_MODELA':BEGIN
      (*!dustem_inputs).grain='GRAIN_G17_MODELA.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELA.DAT'
    END
    'G17_MODELB':BEGIN
      (*!dustem_inputs).grain='GRAIN_G17_MODELB.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELB.DAT'
    END
    'G17_MODELC':BEGIN
      (*!dustem_inputs).grain='GRAIN_G17_MODELC.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELC.DAT'
    END
    'G17_MODELD':BEGIN
      (*!dustem_inputs).grain='GRAIN_G17_MODELD.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELD.DAT'
    END

    'THEMIS':BEGIN
      (*!dustem_inputs).grain='GRAIN_THEMIS.DAT'

      (*!dustem_inputs).align='ALIGN_THEMIS.DAT'

    END

    'NY_MODELA':BEGIN
      (*!dustem_inputs).grain='GRAIN_NY_MODELA.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELA.DAT' 

    END

    ;This is a user-defined model ;This will prolly have to get modified.    

    'USER_MODEL':BEGIN
      (*!dustem_inputs).grain='GRAIN_USER.DAT' 
    END

    ELSE:BEGIN

      message,'model '+model+' unknown',/continue

      message,'Known models are COMPIEGNE_ETAL2010,DBP90,DL01,DL07,AJ13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD,THEMIS',/continue

      stop

    END
  ENDCASE
ENDIF

;stop

;Initialize !dustem_params structure
IF not keyword_set(dir) THEN BEGIN
;  CASE getenv('DUSTEM_WHICH') OF
  CASE !dustem_which OF
    'DESERT':dir_in=!dustem_wrap_soft_dir+'/Data/DESERT_POST_ISO_MORE/'
    'COMPIEGNE': dir_in=!dustem_wrap_soft_dir+'/Data/MC_DAT/'
;    'VERSTRAETE': dir_in=!dustem_wrap_soft_dir+'/Data/d_3.5/'
    'VERSTRAETE': dir_in=!dustem_wrap_soft_dir+'/Data/LV_DAT/'
    'WEB3p8': dir_in=!dustem_soft_dir

    ELSE: dir_in=!dustem_wrap_soft_dir+'/Data/les_DAT/';!dustem_wrap_soft_dir+'/Data/les_DAT/'

  ENDCASE
ENDIF ELSE BEGIN
  dir_in=dir
ENDELSE
;stop
;dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/'+dir
;dir_in=!dustem_soft_dir+'/src/'+dir

if keyword_set(kwords) then defsysv, '!dustem_kwords', ptr_new(kwords) else defsysv, '!dustem_kwords', ptr_new()

st_model=dustem_read_all(dir_in,/silent)
!dustem_params=ptr_new(st_model)

;=== create dynamical storage if needed
IF not file_test(!dustem_dat,/dir) THEN BEGIN

  spawn,'mkdir '+!dustem_dat

ENDIF
IF not file_test(!dustem_res,/dir) THEN BEGIN
  spawn,'mkdir '+!dustem_dat
ENDIF

;=== remove dynamical storage directories
;=== This is needed for iterative use of DUSTEM from IDL
if not keyword_set(noerase) then begin

CASE !dustem_which OF
  'VERSTRAETE':BEGIN
    IF file_test(!dustem_dat,/dir) THEN BEGIN
      spawn,'rm -rf '+!dustem_dat+'/les_DAT'
      spawn,'rm -rf '+!dustem_dat+'/les_QABS'
      spawn,'rm -rf '+!dustem_dat+'/les_CAPA'
    ENDIF
    IF file_test(!dustem_res,/dir) THEN BEGIN
      spawn,'rm -rf '+!dustem_res+'/les_RES'
    ENDIF

  END

  'WEB3p8':BEGIN 

    IF file_test(!dustem_dat,/dir) THEN BEGIN
      spawn,'rm -rf '+!dustem_dat+'/data'
      spawn,'rm -rf '+!dustem_dat+'/hcap'
      spawn,'rm -rf '+!dustem_dat+'/oprop'
      spawn,'rm -rf '+!dustem_dat+'/out'

      ;spawn,'rm -rf '+!dustem_dat+'/plots' TO BE TESTED. NOT SURE ABOUT WHAT I AM DOING HERE. ASK JP ABOUT THIS

    ENDIF

  END
  ELSE: BEGIN
  END
ENDCASE

;=== create storage directories
;=== This is needed for iterative use of DUSTEM from IDL

CASE !dustem_which OF
  'VERSTRAETE':BEGIN
    IF file_test(!dustem_dat,/dir) THEN BEGIN
      spawn,'mkdir '+!dustem_dat+'/les_DAT'
      spawn,'mkdir '+!dustem_dat+'/les_QABS'
      spawn,'mkdir '+!dustem_dat+'/les_CAPA'
    ENDIF
    IF file_test(!dustem_res,/dir) THEN BEGIN
      spawn,'mkdir '+!dustem_res+'/les_RES'
    ENDIF

  END

  'WEB3p8':BEGIN  ; IN CHANGED THIS AS WELL

    IF file_test(!dustem_dat,/dir) THEN BEGIN
      spawn,'mkdir '+!dustem_dat+'/data'
      spawn,'mkdir '+!dustem_dat+'/hcap'
      spawn,'mkdir '+!dustem_dat+'/oprop'
      spawn,'mkdir '+!dustem_dat+'/out'

      ;spawn,'mkdir '+!dustem_dat+'/plots' TO BE TESTED

    ENDIF

  END
  ELSE: BEGIN
  END
ENDCASE
endif ; /noerase

;=== If needed, test runing the model through the wrapper
IF keyword_set(wraptest) THEN BEGIN
  dir_out=!dustem_dat
  dustem_write_all,st_model,dir_out

 ; stop

  st=dustem_run()
  IF keyword_set(plot_it) THEN BEGIN
    loadct,13
    dustem_plot_nuinu_em,st,yr=[1e-14,1.e-6],/ysty,xr=[1,5e3],/xsty
  ENDIF
ENDIF

;stop

;Initialize filters
dustem_filter_init,dir=dir

the_end:

END