PRO dustem_read_fits_table,filename=filename,help=help,dustem_spectra_st=dustem_spectra_st

;+
; NAME:
;    dustem_read_fits_table
; PURPOSE:
;    Reads a Dustem fits table into current Dustem system variables
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_read_fits_table[,filename=][,/help]
; INPUTS:
; OPTIONAL INPUT PARAMETERS:
;    filename      = File name (default='/tmp/dustem_wrap_ouput_final.fits')
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    dustem_spectra_st = emission and extinction spectra (total and for each grain type)
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    dustem system variables are set
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at IRAP
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

; !DUSTEM_DAT = '/tmp/ramdisk8MB/'       
; !DUSTEM_DATA = -> <Anonymous> Array[1]
; !DUSTEM_DO_CC =        1
; !DUSTEM_F90_EXEC = '$HOME/Soft_Librairies/dustem4.3_wk/src/dustem'
; !DUSTEM_FILTERS = <PtrHeapVar570>
; !DUSTEM_FIT = <PtrHeapVar1>
; !DUSTEM_F_HI =       1.00000   
; !DUSTEM_IDL_CONTINUUM =       0.00000
; !DUSTEM_IDL_FREEFREE =       0.00000
; !DUSTEM_IDL_SYNCHROTRON =       0.00000
; !DUSTEM_INPUTS = <PtrHeapVar2>
; !DUSTEM_INSTRUMENT_DESCRIPTION = -> <Anonymous> Array[120]
; !DUSTEM_NEVER_DO_CC =        0 
; !DUSTEM_PARAMS = <PtrHeapVar15>
; !DUSTEM_PARINFO = <PtrHeapVar943>
; !DUSTEM_PREVIOUS_CC = <PtrHeapVar3346>
; !DUSTEM_RES = '/tmp/ramdisk8MB/'
; !DUSTEM_SHOW_PLOT =        1 
; !DUSTEM_SOFT_DIR = '$HOME/Soft_Librairies/dustem4.3_wk/'       
; !DUSTEM_VERBOSE =        1
; !DUSTEM_WHICH = 'WEB3p8'
; !DUSTEM_WRAP_SOFT_DIR = '$HOME/Soft_Librairies/dustem-wrapper_idl/'
; !FIT_RCHI2_WEIGHT = -> <Anonymous> Array[1]
; !RUN_ANIS =       0.00000
; !RUN_CIRC =       0.00000      
; !RUN_IONFRAC =       0.00000
; !RUN_LIN =       1.00000
; !RUN_POL =        1
; !RUN_RRF =       0.00000
; !RUN_TLS =       0.00000
; !RUN_UNIV =       0.00000

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

;write final file
IF keyword_set(filename) THEN BEGIN
	file=filename
ENDIF ELSE BEGIN
	file='/tmp/dustem_wrap_ouput_final.fits'
ENDELSE
message,'Reading '+file,/info

unit_input_sed=1
unit_predicted_sed=2
unit_predicted_spectra=3
unit_predicted_extinctions=4

;==== Write the observed SED in extension 1
str_input_SED=mrdfits(file,unit_input_sed,header_input_sed)
;==== add the predicted SED in extension 2
str_predicted_SED=mrdfits(file,unit_predicted_sed,header_predicted_SED)
;==== add the predicted total spectrum in extension 3
;str_predicted_spectrum_tot=mrdfits(file,3,header_predicted_spectrum_tot)

Nsed=n_elements(str_input_SED)
(*!dustem_data.sed).instru_names=dustem_filter2instru(str_input_SED.filter)
(*!dustem_data.sed).filt_names=str_input_SED.filter
(*!dustem_data.sed).wav=str_input_SED.wavelength
(*!dustem_data.sed).values=str_input_SED.I
(*!dustem_data.sed).sigma=str_input_SED.varianceII

IF !run_pol NE 0 THEN BEGIN

ENDIF

;predicted_I_sed=str_predicted_spectrum_tot.I

best_fit_chi2=sxpar(header_predicted_SED,'CHI2')
best_fit_rchi2=sxpar(header_predicted_SED,'RCHI2')

;===== Fill in header with parameter values
parameters_desc=*(*!dustem_fit).PARAM_DESCS
parameters_init_values=*(*!dustem_fit).PARAM_INIT_VALUES
best_parameters=*(*!dustem_fit).CURRENT_PARAM_VALUES
best_parameters_uncertainties=*(*!dustem_fit).CURRENT_PARAM_ERRORS
best_fit_chi2=(*!dustem_fit).CHI2
best_fit_rchi2=(*!dustem_fit).RCHI2

Nparams=sxpar(header_predicted_SED,'NPARAMS')
parameters_desc=strarr(Nparams)
best_parameters=dblarr(Nparams)
best_parameters_init_values=dblarr(Nparams)
best_parameters_uncertainties=dblarr(Nparams)
FOR i=0L,Nparams-1 DO BEGIN
	parameters_desc[i]=sxpar(header_predicted_SED,'PARN'+strtrim(i+1,2))
	best_parameters[i]=sxpar(header_predicted_SED,'PARV'+strtrim(i+1,2))
	best_parameters_uncertainties[i]=sxpar(header_predicted_SED,'PARU'+strtrim(i+1,2))
	best_parameters_init_values[i]=sxpar(header_predicted_SED,'PARI'+strtrim(i+1,2))
ENDFOR

*(*!dustem_fit).PARAM_DESCS=parameters_desc
*(*!dustem_fit).CURRENT_PARAM_VALUES=best_parameters
*(*!dustem_fit).CURRENT_PARAM_ERRORS=best_parameters_uncertainties
*(*!dustem_fit).PARAM_INIT_VALUES=best_parameters_init_values

;=== Below is only to get the same output form as dustem_compute_sed
;Note: If the same model is used, the outputs should be the same as what is in the fits file ....
dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,out_st=dustem_spectra_st)

;==== get the predicted emission spectra from extension 4
toto=mrdfits(file,unit_predicted_spectra,header_predicted_spectra)
dustem_spectra_st.sed=toto
;==== get the predicted extinction spectra from extension 5
toto=mrdfits(file,unit_predicted_extinctions,header_predicted_extinctions)
dustem_spectra_st.ext=toto

;stop

the_end:

END