PRO dustem_run_prediction_example, model $
                                   ,filters $
                                   ,wave=wave $
                                   ,next=next $
                                   ,sed_outfile=sed_outfile $
                                   ,ext_outfile=ext_outfile $
                                   ,show=show $
                                   ,fits_save=fits_save $
                                   ,noobj=noobj $
                                   ,help=help

;+
; NAME:
;    dustem_run_prediction_example
;
; PURPOSE:
;    This is an example of how to run the DustEM fortran code using
;    DustEMWrap and use it to make prediction for fluxes in a set of filters
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Example
;
; CALLING SEQUENCE:
;    dustem_run_prediction_example,model,filters,show=show,postscript=postscript,help=help
;
; INPUTS:
;    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.
;    filters = specifies the filters for which flux predictions should
;              be calculated. The filters will be shown on the plot.
;  
; OPTIONAL INPUT PARAMETERS:
;    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.

; OUTPUTS:
;    
; OPTIONAL OUTPUT PARAMETERS:
;   sed_outfile = name of file to write predicted SED (in .xcat
;   format). Default is './dustemwrap_predicted_sed.xcat'
;   ext_outfile = name of file to write predicted extinction curve (in .xcat
;   format). Default is './dustemwrap_predicted_ext.xcat'
;
; ACCEPTED KEY-WORDS:
;    help      = if set, print this help
;    show      = if set, show a plot of the predicted SED and the
;                model spectra
;    fits_save = if set, save the fit results in a binary
;               FITS file. 
;    noobj     = if set, runs with no object-oriented graphics
;
; 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_run_prediction_example,'DBP90',['MIRI1','SPIRE3','NIKA21','IRAS1','IRAS2','IRAS3','SPIRE2','SPIRE1','PACS1','WISE3']
;
; MODIFICATION HISTORY:
;    Written JPB Apr-2011
;    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_run_prediction_example'
  goto,the_end
ENDIF

exists=dustem_test_model_exists(model)
if exists ne 1 then $
   message,'Unknown dust model'
polct=dustem_test_model_exists(model,/pol,/silent)

;== INITIALISE DUSTEM
dustem_init,model=model,pol=1
!dustem_verbose=1
!dustem_show_plot=1
!dustem_which='RELEASE'
dir_in=!dustem_soft_dir

;== READ DEFAULT DUST MODEL VALUES
st_model=dustem_read_all(dir_in)

IF keyword_set(sed_outfile) THEN BEGIN
  use_sed_outfile=sed_outfile
ENDIF ELSE BEGIN
  use_sed_outfile='./dustmodel_prediction_sed.xcat'
ENDELSE
IF keyword_set(ext_outfile) THEN BEGIN
  use_ext_outfile=ext_outfile
ENDIF ELSE BEGIN
  use_ext_outfile='./dustmodel_prediction_ext.xcat'
ENDELSE

;=== initialize filters for the SED
Nfilt=n_elements(filters)
sed=dustem_initialize_sed(Nfilt)
sed.filter=filters
sed.wave=dustem_filter2wav(filters)
sed.instru=dustem_filter2instru(filters)

;=== initialize IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro
sed[*].StokesI=1.e-10
sed.StokesQ=sed.StokesI/100.
sed.StokesU=sed.StokesI/100.
sed.SigmaII=sed.StokesI/100.
sed.SigmaQQ=sed.StokesI/100.
sed.SigmaUU=sed.StokesI/100.
sed.SigmaIQ=sed.StokesI/100.
sed.SigmaIU=sed.StokesI/100.
sed.SigmaQU=sed.StokesI/100.

;=== initialize wavelengths for the EXTINCTION
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(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))

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.e-10
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.

dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext

;=== initialize at least one parameter as 'free'
;=== this is only to keep dustemwrap happy in its management of variables
;=== no fitting is going to be done
pd = ['(*!dustem_params).G0'] ; G0
pval = [st_model.g0]

;=== declare some other parameters as fixed
;=== this shows how to add any non-dust plugins to the observed spectra/SED
;=== for the dust-model parameters, it is not necessary, but allows the
;=== user to (i) change default values and (2) see the adopted values in the graphical output windows

fpd=['(*!dustem_params).grains(0).mdust_o_mh',$          ;PAH0 mass fraction
     '(*!dustem_params).grains(1).mdust_o_mh',$          ;VSG mass fraction
     '(*!dustem_params).grains(2).mdust_o_mh'];,$
;     'dustem_plugin_continuum_1', $ ;
;     'dustem_plugin_continuum_2']                      ;                   
fpval = [st_model.grains(0).mdust_o_mh,$
         st_model.grains(1).mdust_o_mh,$
         st_model.grains(2).mdust_o_mh];,$
 ;        800,$
 ;        0.05]

;=== initialise all this information into dustemwrap's brain
dustem_init_params,model,pd,pval,fpd=fpd,fiv=fpval,pol=1

;=== compute the predicted SED and extinction curve for the dust-model
;=== and any plugins 
dustem_Ised=dustem_compute_sed(pval,st=dummy)
dustem_Qsed=dustem_Ised*0.
dustem_Used=dustem_Ised*0.
; uncomment following lines if running a polarisation model and you
;want predictions of QU
;toto=dustem_compute_stokes(pval,st=dummy) ;this procedure also allows for the extraction of the spectra
;dustem_qsed = toto[0]
;dustem_used = toto[1]

dustem_iext=dustem_compute_ext(pval,st=dummy)
dustem_Qext=dustem_Iext*0.
dustem_Uext=dustem_Iext*0.
; uncomment following lines if running a polarisation model and you
;want predictions of QU
;toto=dustem_compute_stokext(pval,st=dummy) ;this procedure also allows for the extraction of the spectra
;dustem_qext = toto[0]
;dustem_uext = toto[1]

;== Now we have a predicted SED (I), we fill up a structure that we
;== will write to the sed_outfile
;== if you 
sed.stokesI=dustem_Ised
sed.stokesQ=dustem_Qsed
sed.stokesU=dustem_Used
;=== set uncertainties for I,Q,U to 0.
sed.sigmaII=0.
sed.sigmaQQ=0.
sed.sigmaUU=0.
;=== set covariances to 0.
sed.sigmaIQ=0.
sed.sigmaIU=0.
sed.sigmaQU=0.

;==== fill in dependent columns of the SED.
sed=dustem_fill_sed_dependent_columns(sed)

;======== save the predicted SED to a file
write_xcat,sed,use_sed_outfile
message,'Wrote '+use_sed_outfile,/continue

;== Now we have a predicted extinction (IQU), we fill up a structure that we
;== will write to the ext_outfile
ext.EXT_I=dustem_Iext
ext.EXT_Q=dustem_Qext
ext.EXT_U=dustem_Uext
;=== set uncertainties to 0.
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_ext_outfile
message,'Wrote '+use_ext_outfile,/continue

;=== MAKE A PLOT TO VISUALISE THE OUTPUT
IF !dustem_noobj THEN BEGIN
  dustemwrap_plot_noobj,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit
 ENDIF ELSE BEGIN
  dustemwrap_plot,*(*!dustem_fit).CURRENT_PARAM_VALUES,st=dummy,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit
ENDELSE


IF keyword_set(fits_save) THEN BEGIN
   message,'Writing out results structure: '+fits_save,/info
   dustem_write_fits_table,filename=fits_save,help=help
ENDIF

the_end:
message,'Finished dustem_run_prediction_example',/info

END