PRO dustem_read_fits_table,filename=filename,dustem_st=dustem_st,help=help

;+
; 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_st = emission and extinction dustem output structure (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

toto=mrdfits(file,0,main_header)

;get extention numbers from main header
current_ext=1
ext_names=['']
exts=[-1]
finished=0
WHILE not finished DO BEGIN
	ext_name=sxpar(main_header,'EXT'+strtrim(current_ext,2),count=count)
	IF count NE 0 THEN BEGIN
		ext_names=[ext_names,ext_name]
		exts=[exts,current_ext]
		current_ext=current_ext+1
	ENDIF ELSE BEGIN
		finished=1
	ENDELSE
ENDWHILE
ext_names=ext_names[1:*]
exts=exts[1:*]
Nextensions=n_elements(exts)

message,'================================================',/info
message,'Reading '+file,/info
message,'File extensions detected are: ',/info
FOR i=0L,Nextensions-1 DO BEGIN
	message,'extension '+strtrim(exts[i],2)+' : '+strtrim(ext_names[i],2),/info
ENDFOR
message,'================================================',/info

;extensions_strs=['INPUT_EMISSION_SED','DUSTEM_PREDICTED_SED', $
;				 'INPUT_EXTINCTION_SED','DUSTEM_PREDICTED_EXTINCTION_SED', $
;				 'DUSTEM_PREDICTED_EMISSION_SPECTRA','DUSTEM_PREDICTED_EXTINCTION_SPECTRA', $
;				 'DUSTEM_PREDICTED_EMISSION_SPECTRA_POL','DUSTEM_PREDICTED_EXTINCTION_SPECTRA_POL']

ind=where(ext_names EQ 'INPUT_EMISSION_SED',count)
IF count NE 0 THEN BEGIN
	unit_input_emission_sed=exts[ind[0]]
	;print,unit_input_emission_sed
ENDIF ELSE BEGIN
    unit_input_emission_sed=-1
ENDELSE
ind=where(ext_names EQ 'DUSTEM_PREDICTED_SED',count)
IF count NE 0 THEN BEGIN
	unit_dustem_predicted_sed=exts[ind[0]]
	;print,unit_dustem_predicted_sed
ENDIF ELSE BEGIN
    unit_dustem_predicted_sed=-1
ENDELSE

ind=where(ext_names EQ 'INPUT_EXTINCTION_SED',count)
IF count NE 0 THEN BEGIN
	unit_input_extinction_sed=exts[ind[0]]
	;print,unit_input_extinction_sed
ENDIF ELSE BEGIN
    unit_input_extinction_sed=-1
ENDELSE
ind=where(ext_names EQ 'DUSTEM_PREDICTED_EXTINCTION_SED',count)
IF count NE 0 THEN BEGIN
	unit_dustem_predicted_extinction_sed=exts[ind[0]]
	;print,unit_dustem_predicted_extinction_sed
ENDIF ELSE BEGIN
	unit_dustem_predicted_extinction_sed=-1   ;That would be abnormal
ENDELSE

ind=where(ext_names EQ 'DUSTEM_PREDICTED_EMISSION_SPECTRA',count)
IF count NE 0 THEN BEGIN
	unit_dustem_predicted_emission_spectra=exts[ind[0]]
	;print,unit_dustem_predicted_emission_spectra
ENDIF ELSE BEGIN
	unit_dustem_predicted_emission_spectra=-1   ;That would be abnormal
ENDELSE
ind=where(ext_names EQ 'DUSTEM_PREDICTED_EXTINCTION_SPECTRA',count)
IF count NE 0 THEN BEGIN
	unit_dustem_predicted_extinction_spectra=exts[ind[0]]
	;print,unit_dustem_predicted_extinction_spectra
ENDIF ELSE BEGIN
	unit_dustem_predicted_extinction_spectra=-1   ;That would be abnormal
ENDELSE

ind=where(ext_names EQ 'DUSTEM_PREDICTED_EMISSION_SPECTRA_POL',count)
IF count NE 0 THEN BEGIN
	unit_dustem_predicted_emission_spectra_pol=exts[ind[0]]
	;print,unit_dustem_predicted_emission_spectra_pol
ENDIF ELSE BEGIN
	unit_dustem_predicted_emission_spectra_pol=-1
ENDELSE
ind=where(ext_names EQ 'DUSTEM_PREDICTED_EXTINCTION_SPECTRA_POL',count)
IF count NE 0 THEN BEGIN
	unit_dustem_predicted_extinction_spectra_pol=exts[ind[0]]
	;print,unit_dustem_predicted_extinction_spectra_pol
ENDIF ELSE BEGIN
	unit_dustem_predicted_extinction_spectra_pol=-1
ENDELSE

;stop

; unit_input_sed=sxpar(main_header,'INPUT_EMISSION_SED')
; unit_predicted_sed=sxpar(main_header,'DUSTEM_PREDICTED_SED')
; unit_predicted_spectra=sxpar(main_header,'DUSTEM_PREDICTED_EMISSION')
; unit_predicted_extinctions=sxpar(main_header,'DUSTEM_PREDICTED_EXTINCTION')
; unit_predicted_spectra_pol=sxpar(main_header,'DUSTEM_PREDICTED_EMISSION_POL')
; unit_predicted_extinctions_pol=sxpar(main_header,'DUSTEM_PREDICTED_EXTINCTION_POL')

;==== read the observed SED (if any)
IF unit_input_emission_sed NE -1 THEN BEGIN
	str_input_SED=mrdfits(file,unit_input_emission_sed,header_input_sed)
ENDIF ELSE BEGIN
ENDELSE
;==== read the predicted SED (if any)
IF unit_dustem_predicted_sed NE -1 THEN BEGIN
	str_predicted_SED=mrdfits(file,unit_dustem_predicted_sed,header_predicted_SED)
ENDIF ELSE BEGIN
ENDELSE
IF unit_input_extinction_sed NE -1 THEN BEGIN
	str_input_EXT=mrdfits(file,unit_input_extinction_sed,header_input_ext)
ENDIF ELSE BEGIN
ENDELSE
IF unit_dustem_predicted_extinction_sed NE -1 THEN BEGIN
	str_predicted_EXT=mrdfits(file,unit_input_extinction_sed,header_input_ext)
ENDIF ELSE BEGIN
ENDELSE

used_model=strtrim(sxpar(main_header,'MODEL'),2)
used_pol=long(strtrim(sxpar(main_header,'POL'),2))
used_version=strtrim(sxpar(main_header,'WRAP_V'),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
IF !dustem_version.version NE used_version THEN BEGIN
	message,'Caution: fits file was not created with the same Dustemwrap version as your current version',/continue
	message,'Fits dustemwrap version is    :'+used_version,/continue
	message,'current dustemwrap version is :'+!dustem_version.version,/continue
	stop
ENDIF

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

;==== get parameter values
Nparams=sxpar(main_header,'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(main_header,'PARN'+strtrim(i+1,2))
	best_parameters[i]=sxpar(main_header,'PARV'+strtrim(i+1,2))
	best_parameters_uncertainties[i]=sxpar(main_header,'PARU'+strtrim(i+1,2))
	best_parameters_init_values[i]=sxpar(main_header,'PARI'+strtrim(i+1,2))
	parameters_func_values[i]=sxpar(main_header,'PARF'+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)

;==== This is to create (*!dustem_data).sed in case it does not exist already
; it doesn't exist at all times because dustem_init has been run
IF unit_input_emission_sed NE -1 THEN BEGIN
	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=strtrim(str_input_SED.filter,2)
	(*(*!dustem_data).sed).wav=str_input_SED.wavelength
	(*(*!dustem_data).sed).values=str_input_SED.STOKESI
	(*(*!dustem_data).sed).sigma=la_power(str_input_SED.varianceII,0.5)

	IF used_pol 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.STOKESQ
		(*(*!dustem_data).qsed).sigma=la_power(str_input_SED.varianceQQ,0.5)
		(*(*!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.STOKESU
		(*(*!dustem_data).used).sigma=la_power(str_input_SED.varianceUU,0.5)
	ENDIF
ENDIF

IF unit_input_extinction_sed NE -1 THEN BEGIN
	Next=n_elements(str_input_EXT)
	ext=dustem_initialize_internal_sed(Next)
	(*!dustem_data).ext=ptr_new(ext)

	(*(*!dustem_data).ext).instru_names=dustem_filter2instru(str_input_EXT.filter)
	(*(*!dustem_data).ext).filt_names=str_input_EXT.filter
	(*(*!dustem_data).ext).wav=str_input_EXT.wavelength
	(*(*!dustem_data).ext).values=str_input_EXT.I
	(*(*!dustem_data).ext).sigma=la_power(str_input_EXT.varianceII,0.5)

	IF used_pol THEN BEGIN
		;stop
		(*!dustem_data).qext=ptr_new(ext)
		(*!dustem_data).uext=ptr_new(ext)
		(*(*!dustem_data).qext).instru_names=dustem_filter2instru(str_input_EXT.filter)
		(*(*!dustem_data).qext).filt_names=str_input_EXT.filter
		(*(*!dustem_data).qext).wav=str_input_EXT.wavelength
		(*(*!dustem_data).qext).values=str_input_EXT.Q
		(*(*!dustem_data).qext).sigma=la_power(str_input_EXT.varianceQQ,0.5)
		(*(*!dustem_data).uext).instru_names=dustem_filter2instru(str_input_EXT.filter)
		(*(*!dustem_data).uext).filt_names=str_input_EXT.filter
		(*(*!dustem_data).uext).wav=str_input_EXT.wavelength
		(*(*!dustem_data).uext).values=str_input_EXT.U
		(*(*!dustem_data).uext).sigma=la_power(str_input_EXT.varianceUU,0.5)
	ENDIF
ENDIF

;=== 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 ....
IF unit_input_emission_sed NE -1 THEN BEGIN
	recomputed_dustem_predicted_emission_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,st=dustem_st)
ENDIF ELSE BEGIN
	recomputed_dustem_predicted_ext_sed=dustem_compute_ext(*(*!dustem_fit).current_param_values,st=dustem_st)
ENDELSE

dustem_predicted_em=mrdfits(file,unit_dustem_predicted_emission_spectra,header_predicted_emission)
dustem_predicted_ext= mrdfits(file,unit_dustem_predicted_extinction_spectra,header_predicted_extinction)

dustem_st.sed=dustem_predicted_em   
dustem_st.ext=dustem_predicted_ext

IF used_pol THEN BEGIN
	toto=mrdfits(file,unit_dustem_predicted_emission_spectra_pol,header_predicted_spectra_pol)
	;help,dustem_st.polsed,toto,/str
	dustem_st.polsed=toto
	toto=mrdfits(file,unit_dustem_predicted_extinction_spectra_pol,header_predicted_extinctions_pol)
    dustem_st.polext=toto
ENDIF

the_end:

END