PRO dustem_init,dir=dir,wraptest=wraptest,plot_it=plot_it,mode=mode,help=help,noerase=noerase,st_model=st_model,pol=pol

;+
; NAME:
;   dustem_init
; CALLING SEQUENCE:
;   dustem_init
; PURPOSE:
;   Defines system variables to run the DUSTEM code
; INPUTS:
;    mode    = 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
;-

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

;=== 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).

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(),       $
                                polsed: ptr_new(),      $
                                polfrac: ptr_new()       $
                         }
        rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.}
    
endif else begin
	defsysv, '!dustem_data', {    $ ;Data to fit
                                sed:    ptr_new(),      $
                                ext:    ptr_new()       $
                         }
	rchi2_weight = {sed: 0.,ext: 0.}
endelse

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)

;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_idl_continuum', 0.

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)
;===This is the Desert option
IF keyword_set(mode) THEN BEGIN
  CASE mode 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


    ELSE:BEGIN
      message,'mode '+mode+' unknown'
    END
  ENDCASE
ENDIF

;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/'
  ENDCASE
ENDIF ELSE BEGIN
  dir_in=dir
ENDELSE
;stop
;dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/'+dir
;dir_in=!dustem_soft_dir+'/src/'+dir
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 getenv('DUSTEM_WHICH') OF
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
;    spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/les_DAT'
;    spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/les_QABS'
;    spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/les_CAPA'
;    spawn,'rm -rf '+getenv('DUSTEM_RES')+'/les_RES'
  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'
    ENDIF
;    spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/data'
;    spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/hcap'
;    spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/oprop'
;    spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/out'
  END
  ELSE: BEGIN
  END
ENDCASE

;=== create storage directories
;=== This is needed for iterative use of DUSTEM from IDL
;CASE getenv('DUSTEM_WHICH') OF
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
;    spawn,'mkdir '+getenv('DUSTEM_DAT')+'/les_DAT'
;    spawn,'mkdir '+getenv('DUSTEM_DAT')+'/les_QABS'
;    spawn,'mkdir '+getenv('DUSTEM_DAT')+'/les_CAPA'
;    spawn,'mkdir '+getenv('DUSTEM_RES')+'/les_RES'
  END
  'WEB3p8':BEGIN
    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'
    ENDIF
;    spawn,'mkdir '+getenv('DUSTEM_DAT')+'/data'
;    spawn,'mkdir '+getenv('DUSTEM_DAT')+'/hcap'
;    spawn,'mkdir '+getenv('DUSTEM_DAT')+'/oprop'
;    spawn,'mkdir '+getenv('DUSTEM_DAT')+'/out'
  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