PRO dustem_init,model=model $
                ,wraptest=wraptest $
                ,show_plots=show_plots $
                ,noerase=noerase $
                ,noobj=noobj $
                ,st_model=st_model $
                ,polarization=polarization $
                ,verbose=verbose $
                ,help=help $
                ,grain_keywords=grain_keywords

;+
; NAME:
;   dustem_init
;  
; PURPOSE:
;   Defines system variables to run the DUSTEM code
;  
; CATEGORY:
;    DustEMWrap, Distributed, HighLevel, Initialization
;  
; CALLING SEQUENCE:
;   dustem_init,model=model
;
; INPUTS:
;    None
;
; OPTIONAL INPUT PARAMETERS:
;    None
;
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUT PARAMETERS:
;    ;   st_model
;
; 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.
;   wraptest       = test running dustem f90 through the wrapper
;   help           = if set, print this help
;   verbose        = if set, output code diagnostic info (default OFF)
;   show_plots     = if set, a plot will be produced for each fitting step
;   polarization   = set this variables for runs taking polarization into account
;   noerase        : set this keyword to not erase the temporary DAT files used by DustEMWrap
;   noobj          : set this keyword to turn off object-oriented
;                    graphical output when using IDL. Object-oriented
;                    graphics are OFF by default for GDL and Fawlty.
;   grain_keywords : Fortran grain keywords. Must be an array of Ngrains strings. Use '?' to keep default values.
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    1/ The following 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 DustEMWrap idl code must be installed
;
; PROCEDURES AND SUBROUTINES USED:
;   
; EXAMPLES
;    dustem_init,model="MC10",/wrap
;
; MODIFICATION HISTORY:
;    Written by JPB,NF,DP Jan-2007
;    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_init'
  goto,the_end
ENDIF

dustem_define_la_common 
device,decomposed=0

;define which code is running, idl, gdl or fawlty
defsysv,'!fl',exist=fawlty
defsysv,'!gdl',exist=gdl
IF fawlty THEN defsysv,'!dustemwrap_which_language','fawlty'
IF gdl THEN defsysv,'!dustemwrap_which_language','gdl'
IF not fawlty and not gdl THEN defsysv,'!dustemwrap_which_language','idl'

;Control of object-oriented graphics
IF fawlty OR gdl THEN defsysv,'!dustem_noobj',1 ELSE defsysv,'!dustem_noobj',0
IF keyword_set(noobj) then defsysv,'!dustem_noobj',1

;see if the user has set their dustem_which
defsysv,'!dustem_which',exist=dwhich
if not dwhich then begin
   defsysv,'!dustem_which','RELEASE'
   message,'!dustem_which not known from idl_startup, using RELEASE version',/info
end

un_mouchard={iteration:0L, $
            chi2:0.d0,rchi2:0.d0, $
            qchi2:0.d0,qrchi2:0.d0, $
            uchi2:0.d0,urchi2:0.d0, $
            parameters_descriptions:ptr_new(),parameters_values:ptr_new()}
mouchard=ptr_new(replicate(un_mouchard,1))
defsysv,'!iteration_mouchard',mouchard

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

defsysv,'!dustem_show_plot',0   ;1=show all iterations of fitting
IF keyword_set(show_plots) THEN !dustem_show_plot=1

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

;version={version:'V1.2',version_number:1.2}
;version={version:'V2.0',version_number:2.0}
version={version:'V4.3',version_number:4.3} ; March 2023
defsysv, '!dustem_version',version       ;This is the current dustemwrap version

;=== 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 recommended, but faster).

;=== polarization system variables
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

; default column density assumed for translating fortran output to an emission spectrum
defsysv, '!dustem_HCD', ptr_new(1.00E20) ;I do not see the reason behind using a pointer in this case.

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

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

;=== read info about instrument
file=!dustem_wrap_soft_dir+'instrument_description.xcat'
st=read_xcat(file,/silent)
defsysv,'!dustem_instrument_description',st

defsysv, '!dustem_isrf_file', ptr_new() ;but there is only one system variable for this... Is this ok?

defsysv, '!dustem_iter', {prv:1,act:1}

;define !dustem_fit
dustem_fit_st={data:ptr_new(), $ ;,wavelength:ptr_new(), ;because wavelength arrays in the different datasets (data) can be different
               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, $
               qchi2:0.D0,qrchi2:0.D0, $
               uchi2:0.D0,urchi2:0.D0, $
               current_param_values:ptr_new(), current_param_errors:ptr_new() $
               }
defsysv, '!dustem_fit', ptr_new(dustem_fit_st)    ;Data to fit

dustem_data_st = dustem_define_dustem_data(polarization=polarization,rchi2_weight =rchi2_weight)
defsysv, '!dustem_data', ptr_new(dustem_data_st) 
defsysv,'!dustem_show', ptr_new(dustem_data_st)   

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 'FPOLEXT' 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', {"
instr="fit_rchi2_weight_str={"
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)
defsysv, '!fit_rchi2_weight',ptr_new(fit_rchi2_weight_str)

  ;#########system variables necessary for cgwindow plotting##########     
defsysv, '!dustemcgwin_id', {    $ ;IDs of windows to plot                  
                            sed:    la_undef(),      $                          
                            ext:    la_undef(),      $
                            prms:   la_undef(),      $
                            plgns:  la_undef()       $
                     }  
                     
                     
defsysv, '!dustem_dim', 0                                                                        

;after the new modifications I think that some of these pointers are not needed
defsysv, '!dustemcgwin_ncmds', {    $ ;Data to fit
                            sed:     {pl:0, nrm:0},      $ ; 
                            ext:     {pl:0, nrm:0},      $ ;
                            qsed:    {pl:0, nrm:0},      $ ; 
                            used:    {pl:0, nrm:0},      $ ; 
                            polsed:  {pl:0, nrm:0},      $ ; 
                            polfrac: {pl:0},             $ ;   
                            psi_em:  {pl:0},             $ ; 
                            qext:    {pl:0, nrm:0},      $ ; 
                            uext:    {pl:0, nrm:0},      $ ;
                            polext:  {pl:0, nrm:0},      $ ;
                            fpolext: {pl:0, nrm:0},      $ ;  
                            psi_ext: {pl:0},             $ ;
                            chi2 :   {pl:0, txtwdth:0.}, $ ; ;beacause chi2 needs to be updated (needing the command index)
                            rchi2 :  {pl:0, txtwdth:0.}, $ ;
                            pltit :  {pl:0, txtwdth_m:0.,txtwdth_x:0.}, $ ; plot title
                            prms :   {pl:0},             $ ;
                            plgns :  {pl:0},             $ ;
                            runs :   {pl:0, txtwdth:0.}  $ ;
                     }


;NB: #This is definitely not a robust way to handle the _extra tags but it's a primariy workaround for the plotting tags of the data sets. 
defsysv, '!dustem_plot_range', {    $

                            sed:     {xr:[1.,5.0e5], yr:[5e-8,1.00e6]},      $ ; 
                            qsed:    {xr:[1.,5.0e5], yr:[5e-8,1.00e6]},      $ ; 
                            used:    {xr:[1.,5.0e5], yr:[5e-8,1.00e6]},      $ ; 
                            polsed:  {xr:[1.,5.0e5], yr:[5e-8,1.00e6]},      $ ; 
                            polfrac: {xr:[1.,5.0e5], yr:[1.00E-04,50.00]},      $ ;   
                            psi_em:  {xr:[1.,5.0e5], yr:[-90.,90.]},   $ ;
                            ext:     {xr:[0.01,30], yr:[5.00E-8,10]},      $ ; 
                            qext:    {xr:[0.01,30], yr:[5.00E-8,10]},      $ ; 
                            uext:    {xr:[0.01,30], yr:[5.00E-8,10]},      $ ;
                            polext:  {xr:[0.01,30], yr:[5.00E-8,10]},      $ ;
                            fpolext: {xr:[0.01,30], yr:[1.00E-04,10.00]},      $ ; 
                            psi_ext: {xr:[0.01,30], yr:[-90.,90.]},   $ ;
                            title_m : 'DUSTEMWRAP '+'('+!dustem_version.version+') '+'Spectral Energy Distribution', $ ; default title if no title is specified
                            title_x : 'DUSTEMWRAP '+'('+!dustem_version.version+') '+'Dust Optical Depth' $
                     }

defsysv, '!dustem_mlog', 0      ;necessary for the plotting of negative values on the axis in dstmwrp_exp.pro

;###################################################################


;These two sysvars allow for the polarization angle (whether in emission or extinction) to
;be updated in procedures (non-plugins) that create polarization data such as 
;dustem_compute_stokes and dustem_compute_stokext
defsysv, '!dustem_psi', 0.
defsysv, '!dustem_psi_ext', 0.


;useful pointers (current dustem output structure and current mpfit params values)
defsysv,'!dustem_current',ptr_new()
; defsysv,'!dustem_current_params',ptr_new() ;Is replaced by (*!dustem_fit).CURRENT_PARAM_VALUES 
; defsysv, '!dustem_current_errors', ptr_new() ;Is replaced by (*!dustem_fit).current_param_errors=ptr_new(p_er_dim)


defsysv, '!dustem_end', 0 ;Easiest and fastest way I could find to signal the end of the fit. If there is a better way please remove this.                     
;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 RELEASE version of the code, only in VERSTRAETE
;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]
;; ... etc


;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)
use_model='USER_MODEL'

IF keyword_set(model) THEN BEGIN
   use_model = strupcase(model)
   exists=dustem_test_model_exists(use_model,/silent)
   polexists=dustem_test_model_exists(use_model,/pol,/silent)
   if exists ne 1 then $
      message,'Unknown dust model '+use_model
   if keyword_set(polarization) and polexists eq 0 then begin
      message,'You have requested polarization mode, but the '+use_model+' model does not predict polarization',/info
      stop
   end
END

CASE use_model OF
      'MC10':BEGIN
         (*!dustem_inputs).grain='GRAIN_MC10.DAT'
      END
      'C10':BEGIN
         (*!dustem_inputs).grain='GRAIN_MC10.DAT' ; for convenience
      END
      'COMPIEGNE_ETAL10':BEGIN  ; for backwards compatibility
         (*!dustem_inputs).grain='GRAIN_MC10.DAT'
      END
      'DBP90':BEGIN
         (*!dustem_inputs).grain='GRAIN_DBP90.DAT'
      END
      'DL01':BEGIN              ; for backwards compatibility
         (*!dustem_inputs).grain='GRAIN_DL01.DAT'
      END
      'LD01':BEGIN              ; for backwards compatibility
         (*!dustem_inputs).grain='GRAIN_DL01.DAT'
      END
      'DL01_RV3P1_BC6':BEGIN
         (*!dustem_inputs).grain='GRAIN_DL01.DAT'
      END
;      'DL01_RV5P5B_BC3':BEGIN  ; removed from DUSTEM fortran in Feb2023 release
;         (*!dustem_inputs).grain='GRAIN_DL01_Rv5.5B_bc3.DAT'
;      END
;      'DL01_RV5P5B_BC0':BEGIN  ; removed from DUSTEM fortran in Feb2023 release
;         (*!dustem_inputs).grain='GRAIN_DL01_Rv5.5B_bc0.DAT'
;      END
      'WD01_RV5P5B':BEGIN       ; for backwards compatibility
;         (*!dustem_inputs).grain='GRAIN_DL01_Rv5.5B_bc3.DAT'
          (*!dustem_inputs).grain='GRAIN_WD01_Rv5.5B.DAT'
      END
      'WD01':BEGIN       ; for backwards compatibility
;         (*!dustem_inputs).grain='GRAIN_DL01_Rv5.5B_bc3.DAT'
          (*!dustem_inputs).grain='GRAIN_WD01_Rv5.5B.DAT'
      END
      'DL07':BEGIN
         (*!dustem_inputs).grain='GRAIN_DL07.DAT'
      END
      'J13':BEGIN
         (*!dustem_inputs).grain='GRAIN_J13.DAT'
      END
      'AJ13':BEGIN              ; for backwards compatibility
         (*!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
      'Y24':BEGIN            
         (*!dustem_inputs).grain='GRAIN_Y24.DAT'
         (*!dustem_inputs).align='ALIGN_Y24.DAT'
      END
       'THEMIS':BEGIN              
         (*!dustem_inputs).grain='GRAIN_J13.DAT'
         (*!dustem_inputs).align='ALIGN_J13.DAT'
      END
       'THEMIS2':BEGIN       
         (*!dustem_inputs).grain='GRAIN_Y24.DAT'
         (*!dustem_inputs).align='ALIGN_Y24.DAT'
      END
       'HD22':BEGIN        
          (*!dustem_inputs).grain='GRAIN_HD22.DAT'
          (*!dustem_inputs).align='ALIGN_HD22.DAT' 
       END
       'ASTRODUST':BEGIN        
          (*!dustem_inputs).grain='GRAIN_HD22.DAT'
          (*!dustem_inputs).align='ALIGN_HD22.DAT' 
       END
; Next is a user-defined model. In this case,
; The user is responsible for modifying the fortran GRAIN.DAT and
; ALIGN.DAT files and providing ther content
      'USER_MODEL':BEGIN  
         (*!dustem_inputs).grain='GRAIN.DAT'
         (*!dustem_inputs).align='ALIGN.DAT'
      END
      ELSE:BEGIN
         message,'Code should never reach here',/continue
         stop
      END
   ENDCASE
;ENDIF
defsysv, '!dustem_model', use_model

;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/LV_DAT/'
    'RELEASE': 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

;reading fortran files to set up !dustem_params
defsysv, '!dustem_kwords', ptr_new()
st_model=dustem_read_all(dir_in,/silent)
!dustem_params=ptr_new(st_model)
Ngrains=(*!dustem_params).Ngrains
needed_keywords=strarr(Ngrains)
;replace dustem keywords requested by the user
IF keyword_set(grain_keywords) THEN BEGIN
  IF n_elements(grain_keywords) NE Ngrains THEN BEGIN
    message,'Grain_keywords provided do not match the number of grains ('+strtrim(Ngrains,2)+') in the model',/info
    message,'Please provide '+strtrim(Ngrains,2)+' keywords, using ? for the ones that should be left at their default value',/info
    stop
  ENDIF ELSE BEGIN
    FOR i=0L,Ngrains-1 DO BEGIN
      IF grain_keywords[i] NE '?' THEN BEGIN
        needed_keywords[i]=grain_keywords[i]
      ENDIF ELSE BEGIN
        needed_keywords[i]=(*!dustem_params).grains(i).TYPE_KEYWORDS
      ENDELSE
    ENDFOR
  ENDELSE
ENDIF ELSE BEGIN
  FOR i=0L,Ngrains-1 DO BEGIN
    needed_keywords[i]=(*!dustem_params).grains(i).TYPE_KEYWORDS
  ENDFOR
ENDELSE
;we reread the fortran inputs again
;the dustem_read_ subroutines take care of putting the user-defined keywords into
;the structure via the !dustem_kwords system variable
defsysv,'!dustem_kwords', ptr_new(needed_keywords)
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
  force_mkdir,!dustem_dat
ENDIF
IF not file_test(!dustem_res,/dir) THEN BEGIN
  ;spawn,'mkdir '+!dustem_dat
  force_mkdir,!dustem_dat
ENDIF

;=== remove dynamical storage directories from previous runs
;=== 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
    'RELEASE':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
    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
    'RELEASE':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'
      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
  st=dustem_run()
  yr=[1e-14,1.e-6]
  xr=[1,5e3]
  tit='dustem_init wrapper test '+!dustem_model
  xtit=textoidl('\lambda (\mum)')
  ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
  yr=[1e-11,6.e-7]
  xr=[1,5e3]
  dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit
ENDIF

;Initialize filters
dustem_filter_init

the_end: ; of dustem_init
END