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 DustEMWrap system variables
;
; CATEGORY:
;    DustEMWrap, High-Level, Distributed, User convenience
;
; CALLING SEQUENCE:
;    dustem_read_fits_table[,filename=][,/help]
;
; INPUTS:
;
; OPTIONAL INPUT PARAMETERS:
;    filename      = File name to be read (default='./dustemwrap_results.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:
;    DustEMWrap system variables are set
;
; RESTRICTIONS:
;    The DustEMWrap IDL code must be installed
;
; PROCEDURE:
;    None
;
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard July 2022
;    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_read_fits_table'
  goto,the_end
ENDIF

;write final file
IF keyword_set(filename) THEN BEGIN
	file=filename
ENDIF ELSE BEGIN
	file='./dustemwrap_results.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)

used_model=strtrim(sxpar(header_predicted_SED,'MODEL'),2)
used_pol=long(strtrim(sxpar(header_predicted_SED,'POL'),2))
IF used_model NE 'DEFAULT' THEN BEGIN
	dustem_init,model=used_model,pol=used_pol
ENDIF ELSE BEGIN
  	dustem_init,model=0
ENDELSE
!dustem_model=used_model

;stop
;==== This is to create (*!dustem_data).sed in case it does not exist already
Nsed=n_elements(str_input_SED)
sed=dustem_initialize_internal_sed(Nsed)
(*!dustem_data).sed=ptr_new(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
	;stop
	(*!dustem_data).qsed=ptr_new(sed)
	(*!dustem_data).used=ptr_new(sed)
	(*(*!dustem_data).qsed).instru_names=dustem_filter2instru(str_input_SED.filter)
	(*(*!dustem_data).qsed).filt_names=str_input_SED.filter
	(*(*!dustem_data).qsed).wav=str_input_SED.wavelength
	(*(*!dustem_data).qsed).values=str_input_SED.Q
	(*(*!dustem_data).qsed).sigma=str_input_SED.varianceQQ
	(*(*!dustem_data).used).instru_names=dustem_filter2instru(str_input_SED.filter)
	(*(*!dustem_data).used).filt_names=str_input_SED.filter
	(*(*!dustem_data).used).wav=str_input_SED.wavelength
	(*(*!dustem_data).used).values=str_input_SED.U
	(*(*!dustem_data).used).sigma=str_input_SED.varianceUU
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)
parameters_func_values=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))
	parameters_func_values[i]=sxpar(header_predicted_SED,'PARI'+strtrim(i+1,2))
ENDFOR

(*!dustem_fit).PARAM_DESCS=ptr_new(parameters_desc)
(*!dustem_fit).PARAM_INIT_VALUES=ptr_new(best_parameters_init_values)
(*!dustem_fit).CURRENT_PARAM_VALUES=ptr_new(best_parameters)
(*!dustem_fit).CURRENT_PARAM_ERRORS=ptr_new(best_parameters_uncertainties)
(*!dustem_fit).CHI2=best_fit_chi2
(*!dustem_fit).RCHI2=best_fit_rchi2
(*!dustem_fit).PARAM_FUNC=ptr_new(parameters_func_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 ....
;stop
dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,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