PRO dustem_init,model=model $
                ,wraptest=wraptest $
                ,plot_it=plot_it $
                ,noerase=noerase $
                ,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, High-Level, Initialization
; CALLING SEQUENCE:
;   dustem_init,model=model,
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    ;   st_model *** COMMENT AH -- what is this? END COMMENT AH ***
; ACCEPTED KEY-WORDS:
;    model = specifies the interstellar dust mixture used by DustEM
;           'MC10' model from Compiegne et al 2010 (default)
;           'DBP90' model from Desert et al 1990
;           'DL01' model from Draine & Li 2001
;           'WD01_RV5P5B' model from Weingartner & Draine 2002 with Rv=5.5
;           'DL07' model from Draine & Li 2007
;           'J13' model from Jones et al 2013, as updated in
;                 Koehler et al 2014
;           'G17_ModelA' model A from Guillet et al (2018). Includes
;                 polarisation. See Tables 2 and 3 of that paper for details.
;           'G17_ModelB' model B from Guillet et al (2018)
;           'G17_ModelC' model C from Guillet et al (2018)
;           'G17_ModelD' model A from Guillet et al (2018)
;   wraptest       = test running dustem f90 through the wrapper
;   plot_it        = plots result of wraptest (ignored if wraptest not set)
;   help           = if set, print this help
;   verbose        = if set, output code diagnostic info (default OFF)
;   polarization   = set this variables for runs taking polarization into account
;   noerase        : set this keyword to not erase the temporary memory used by dustemwrap
;   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 dustem idl wrapper must be installed
; PROCEDURES AND SUBROUTINES USED:
;    *** COMMENT AH --> is this really NONE? ****
; EXAMPLES
;    dustem_init,model="MC10",/wrap,/plot_it
; 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 

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

defsysv,'!run_pol',0             ;0=no polar
IF keyword_set(polarization) 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).

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

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


;=== read info about instrument
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.



;NOT REALLY SURE ABOUT THIS - getting the iteration number
defsysv, '!dustem_iter', {prv:1,act:1}
; defsysv, '!dustem_fit_perror', 0.

;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, $
               current_param_values:ptr_new(), current_param_errors:ptr_new() $
               }
defsysv, '!dustem_fit', ptr_new(dustem_fit_st)    ;Data to fit

if keyword_set(polarization) then begin
    
    dustem_data_st ={    $ ;Structure contaning the data to fit
                                sed:    ptr_new(),      $ ; 
                                ext:    ptr_new(),      $ ;
                                qsed: ptr_new(),        $ ; 
                                used: ptr_new(),        $ ; 
                                polsed: ptr_new(),      $ ; 
                                polfrac: ptr_new(),     $ ;   
                                psi_em: ptr_new(),      $ ; 
                                qext: ptr_new(),        $ ; 
                                uext: ptr_new(),        $ ;
                                polext: ptr_new(),      $ ;
                                fpolext: ptr_new(),     $ ;  
                                psi_ext: ptr_new()      $ ; 
                                 
                         }

	defsysv, '!dustem_data', ptr_new(dustem_data_st) 
    
    defsysv,'!dustem_show', ptr_new(dustem_data_st)   
    
    rchi2_weight = {sed: 0.,ext: 0.,qsed:0.,used:0.,polsed:0.,polfrac:0.,psi_em:0.,qext:0.,uext:0.,polext:0.,psi_ext:0.}
         
endif else begin

    dustem_data_st = {    $ 
                                sed:    ptr_new(),      $
                                ext:    ptr_new()       $
                         }
	
	defsysv, '!dustem_data', ptr_new(dustem_data_st)
	
	
	defsysv, '!dustem_show', ptr_new(dustem_data_st)
    
    rchi2_weight = {sed: 0.,ext: 0.}

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 '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', {"
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(),      $
                            prms:   la_undef(),      $
                            plgns:  la_undef()       $
                     }                                                      

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:0.}, $ ; plot title
                            prms :   {pl:0},             $ ;
                            plgns :  {pl:0},             $ ;
                            runs :   {pl:0, txtwdth:0.}  $ ;
                            
                     }


version={version:'V2.0',version_number:2.0}
defsysv, '!dustem_version',version       ;This is the current dustemwrap version


;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_extra', {    $

                            sed:     {pl:{xr:[0,0], yr:[0,0],xlog:1, ylog:1},nrm:{xr:[0,0], yr:[0,0],xlog:1, ylog:0}},      $ ; 
                            qsed:    {pl:{xr:[0,0], yr:[0,0],xlog:1, ylog:1},nrm:{xr:[0,0], yr:[0,0],xlog:1, ylog:0}},      $ ; 
                            used:    {pl:{xr:[0,0], yr:[0,0],xlog:1, ylog:1},nrm:{xr:[0,0], yr:[0,0],xlog:1, ylog:0}},      $ ; 
                            polsed:  {pl:{xr:[0,0], yr:[0,0],xlog:1, ylog:1},nrm:{xr:[0,0], yr:[0,0],xlog:1, ylog:0}},      $ ; 
                            polfrac: {pl:{xr:[0,0], yr:[0,0],xlog:1, ylog:1},nrm:{xr:[0,0], yr:[0,0],xlog:1, ylog:0}},      $ ;   
                            psi_em:  {pl:{xr:[0,0], yr:[-90,90],xlog:1, ylog:0}}                                     ,      $ ;
                            ext:     {pl:{xr:[0,0], yr:[0,0],xlog:1, ylog:1},nrm:{xr:[0,0], yr:[0,0],xlog:1, ylog:0}},      $ ; 
                            qext:    {pl:{xr:[0,0], yr:[0,0],xlog:1, ylog:1},nrm:{xr:[0,0], yr:[0,0],xlog:1, ylog:0}},      $ ; 
                            uext:    {pl:{xr:[0,0], yr:[0,0],xlog:1, ylog:1},nrm:{xr:[0,0], yr:[0,0],xlog:1, ylog:0}},      $ ;
                            polext:  {pl:{xr:[0,0], yr:[0,0],xlog:1, ylog:1},nrm:{xr:[0,0], yr:[0,0],xlog:1, ylog:0}},      $ ;
                            fpolext: {pl:{xr:[0,0], yr:[0,0],xlog:1, ylog:1},nrm:{xr:[0,0], yr:[0,0],xlog:1, ylog:0}},      $ ; 
                            psi_ext: {pl:{xr:[0,0], yr:[-90,90],xlog:1, ylog:0}}                                     ,      $ ;
                            title : 'DUSTEMWRAP '+'('+!dustem_version.version+') '+'Spectral Energy Distribution'           $ ; default title if no title is specified
                            
                     }


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]
;;    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
    'MC10':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
    'WD01_RV5P5B':BEGIN
      (*!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
    '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 ;TB removed?
      (*!dustem_inputs).grain='GRAIN_NY_MODELA.DAT'
      (*!dustem_inputs).align='ALIGN_G17_MODELA.DAT' 
    END
    ;This is a user-defined model ;This will probably need to be modified.    
    'USER_MODEL':BEGIN ;TB removed? 
      (*!dustem_inputs).grain='GRAIN_USER.DAT' 
    END
    ELSE:BEGIN
      message,'model '+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
  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/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

;stop

;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',/continue
    message,'Please provide '+strtrim(Ngrains,2)+' keywords, using ? for those left to their default value',/continue
    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
;reading fortran inputs again, but now with replaced keywords
defsysv,'!dustem_kwords', ptr_new(needed_keywords)
st_model=dustem_read_all(dir_in,/silent)
!dustem_params=ptr_new(st_model)
;stop

;=== 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
    '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'
        ;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
    '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'
        ;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
  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

;Initialize filters
dustem_filter_init

the_end:

END