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,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