PRO dustem_read_fits_table,filename=filename,dustem_st=dustem_st,help=help,Q_spec=Q_spec,U_spec=U_spec

;+
; 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','DUSTEM_SHOWN_SED', $
;				 'INPUT_EXTINCTION_SED','DUSTEM_PREDICTED_EXTINCTION_SED','DUSTEM_SHOWN_EXTINCTION_SED', $
;				 'DUSTEM_PREDICTED_EMISSION_SPECTRA','DUSTEM_PREDICTED_EXTINCTION_SPECTRA', $
;				 'DUSTEM_PREDICTED_EMISSION_SPECTRA_POL','DUSTEM_PREDICTED_EXTINCTION_SPECTRA_POL']
ind_plugin=where(strmid(ext_names,0,strlen('PLUGIN_')) EQ 'PLUGIN_',Nplugin)
IF Nplugin NE 0 THEN BEGIN
	ext_names_plugins=strarr(Nplugin)
	;names_plugins=strarr(Nplugin)
        FOR i=0L,Nplugin-1 DO BEGIN
		ext_names_plugins[i]=ext_names[ind_plugin[i]]
		;names_plugins[i]=strmid(ext_names[ind_plugin[i]],strlen('PLUGIN_'),100)
	ENDFOR
ENDIF

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 'DUSTEM_SHOWN_SED',count)
IF count NE 0 THEN BEGIN
	unit_dustem_shown_sed=exts[ind[0]]
	;print,unit_dustem_predicted_sed
ENDIF ELSE BEGIN
    unit_dustem_shown_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_SHOWN_EXTINCTION_SED',count)
IF count NE 0 THEN BEGIN
	unit_dustem_shown_extinction_sed=exts[ind[0]]
	;print,unit_dustem_predicted_extinction_sed
ENDIF ELSE BEGIN
	unit_dustem_shown_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
;=== find out logical units for plugins
IF Nplugin NE 0 THEN BEGIN
	unit_plugins=lonarr(Nplugin)
	FOR i=0L,Nplugin-1 DO BEGIN
		ind=where(ext_names EQ ext_names_plugins[i],count)
		IF count NE 0 THEN BEGIN
			unit_plugins[i]=exts[ind[0]]
		ENDIF ELSE BEGIN
			unit_plugins[i]=-1
		ENDELSE
	ENDFOR
ENDIF

;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
;==== read the shown SED (if any)
;stop
IF unit_dustem_shown_sed NE -1 THEN BEGIN
	str_shown_SED=mrdfits(file,unit_dustem_shown_sed,header_shown_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
IF unit_dustem_shown_extinction_sed NE -1 THEN BEGIN
	str_shown_EXT=mrdfits(file,unit_dustem_shown_extinction_sed,header_shown_ext)
ENDIF ELSE BEGIN
ENDELSE
;=== read plugin informations and plugins headers
IF Nplugin NE 0 THEN BEGIN
	str_plugin=ptrarr(Nplugin)
	headers_plugin=ptrarr(Nplugin)
	FOR i=0L,Nplugin-1 DO BEGIN
		IF unit_plugins[i] NE -1 THEN BEGIN
			a=mrdfits(file,unit_plugins[i],header)
			str_plugin[i]=ptr_new(a)
			headers_plugin[i]=ptr_new(header)
		ENDIF
	ENDFOR
ENDIF

;stop

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

IF unit_input_emission_sed NE -1 THEN defsysv,'!dustem_hcd',ptr_new(float(sxpar(header_input_sed,'NH_SED')))
IF unit_input_extinction_sed NE -1 THEN defsysv,'!dustem_hcd',ptr_new(float(sxpar(header_input_ext,'NH_EXT')))
IF unit_input_emission_sed NE -1 AND unit_input_extinction_sed NE -1 THEN BEGIN
   defsysv,'!dustem_hcd',ptr_new(float(sxpar(header_input_sed,'NH_SED')))
   if float(sxpar(header_input_sed,'NH_SED')) ne float(sxpar(header_input_ext,'NH_EXT')) then begin
      message,'Assumed column density for emission in input FITS is [H/cm2]: '+string(sxpar(header_input_sed,'NH_SED')),/info
      message,'Assumed column density for extinction in input FITS is [H/cm2]: '+string(sxpar(header_input_sed,'NH_EXT')),/info
      message,'Restoring using NH_SED value for !DUSTEM_HCD',/info
   end
ENDIF

   
!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)
;==== get fixed parameter values, if any
Nfixedparams=sxpar(main_header,'NFPARAMS')
IF Nfixedparams NE 0 THEN BEGIN
	fixed_parameters_desc=strarr(Nfixedparams)
	fixed_parameters_values=dblarr(Nfixedparams)
	FOR i=0L,Nfixedparams-1 DO BEGIN
		fixed_parameters_desc[i]=sxpar(main_header,'FPARN'+strtrim(i+1,2))
		fixed_parameters_values[i]=sxpar(main_header,'FPARV'+strtrim(i+1,2))
	ENDFOR
	(*!dustem_fit).fixed_param_descs=ptr_new(fixed_parameters_desc)
	(*!dustem_fit).fixed_param_init_values=ptr_new(fixed_parameters_values)
ENDIF
;stop

;=== set !dustem_plugin according to fits file content
IF Nplugin NE 0 THEN BEGIN
  dustem_init_plugins,parameters_desc,fpd=fixed_parameters_desc
  FOR i=0L,Nplugin-1 DO BEGIN
  	(*!dustem_plugin).(i).spec=ptr_new(str_plugin[i])
  	header=*headers_plugin[i]
  	pluggin_name=strtrim(sxpar(header,'PLUGIN'),2)
  	pluggin_scope=strtrim(sxpar(header,'SCOPE'),2)
    paramtags=['']
    ii=1 & cc=1
  	WHILE cc NE 0 DO BEGIN
  		ss=strtrim(sxpar(header,'TAGN'+strtrim(ii,2),count=cc),2)
  		IF cc NE 0 THEN paramtags=[paramtags,ss]
  		ii=ii+1
  	ENDWHILE
  	paramtags=paramtags[1:*]
  	(*!dustem_plugin).(i).scope=ptr_new(pluggin_scope)
  	(*!dustem_plugin).(i).PARAMTAG=ptr_new(paramtags)
  ENDFOR
ENDIF

;==== 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_fillup_systvar_from_fits,*!dustem_data,str_input_SED,used_pol
ENDIF

;help,str_shown_SED

;==== This is to create (*!dustem_show).sed in case it does not exist already
IF unit_dustem_shown_sed NE -1 THEN BEGIN
	;Nsed=n_elements(str_shown_SED)
	;sed=dustem_initialize_internal_sed(Nsed)
	dustem_fillup_systvar_from_fits,*!dustem_show,str_shown_SED,used_pol
ENDIF

IF unit_input_extinction_sed NE -1 THEN BEGIN
	;Next=n_elements(str_input_EXT)
	;ext=dustem_initialize_internal_sed(Next)
	dustem_fillup_systvar_from_fits,*!dustem_data,str_input_EXT,used_pol
ENDIF

IF unit_dustem_shown_extinction_sed NE -1 THEN BEGIN
	dustem_fillup_systvar_from_fits,*!dustem_show,str_shown_EXT,used_pol
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)
	;This is to compute the spectrum of Q,U which is not stored in dustem_st (only polsed is)
	if used_pol then toto = dustem_compute_stokes(*(*!dustem_fit).current_param_values,st=dustem_st,Q_spec=Q_spec,U_spec=U_spec,PSI_spec=PSI_spec,dustem_psi_em=dustem_psi_em)
ENDIF ELSE BEGIN
	recomputed_dustem_predicted_ext_sed=dustem_compute_ext(*(*!dustem_fit).current_param_values,st=dustem_st)
ENDELSE
;IF !run_pol AND !run_lin AND tag_exist(*!dustem_show,'qsed') THEN BEGIN
;     dustem_polsed = dustem_compute_polsed(p_dim,st=st,P_spec=P_spec,SP_spec=SP_spec,dustem_polfrac=dustem_polfrac)
;     toto = dustem_compute_stokes(p_dim,st=st,Q_spec=Q_spec,U_spec=U_spec,PSI_spec=PSI_spec,dustem_psi_em=dustem_psi_em) 
;     dustem_qsed = toto[0]
;     dustem_used = toto[1]   
;ENDIF         

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