PRO dustem_write_fits_table,filename=filename,help=help

;+
; NAME:
;    dustem_write_fits_table
;
; PURPOSE:
;    Write a FITS table with content of current DustEMWrap system variables
;
; CATEGORY:
;    DustEMWrap, Distributed, High-Level, User Convenience
;
; CALLING SEQUENCE:
;    dustem_write_fits_table,file[,/help]
;
; INPUTS:
;    file      = File name to be written (default='./dustemwrap_results.fits')
;
; OPTIONAL INPUT PARAMETERS:
;    None
;
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUT PARAMETERS:
;    None
;
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    Output FITS file is written to disk.
;
; 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_write_fits_table'
  goto,the_end
ENDIF

;===== run dustem to get predicted spectrum and sed
IF ptr_valid((*!dustem_data).sed) THEN $
	dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,out_st=dustem_spectra_st)

IF !run_pol THEN BEGIN
;This is needed to compute Q,U of both SED and spectra
   dustem_predicted_polsed=dustem_compute_stokes(*(*!dustem_fit).current_param_values, $
                                                 Q_spec=dustem_predicted_Q_spec,U_spec=dustem_predicted_U_spec,PSI_spec=dustem_predicted_PSI_spec, $
                                                 dustem_psi_em=dustem_psi_em)
   dustem_predicted_Qsed=dustem_predicted_polsed[0]
   dustem_predicted_Used=dustem_predicted_polsed[1]
;Same in extinction, but does not work in polarization extinction was not fitted.
;   dustem_predicted_polext=dustem_compute_stokext(*(*!dustem_fit).current_param_values,bidon, $
;											  dustem_predicted_Qext,dustem_predicted_Uext, $
;											  dustem_predicted_Qext_spec,dustem_predicted_Uext_spec,dustem_predicted_PSI_ext, $
;											  dustem_psi_ext,out_st=out_st)
ENDIF

used_model=!dustem_model

Ngrains=n_tags(dustem_spectra_st.sed)-2

;===== define structures containing results
one_str_input_SED={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5),varianceII:la_undef(5),varianceQQ:la_undef(5),varianceUU:la_undef(5)}
one_str_predicted_SED={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5)}
;one_str_predicted_spectrum_tot={Wavelength:0.,I:0.,Q:0.,U:0.}
;one_str_predicted_spectrum_grain1={Wavelength:0.,I:0.,Q:0.,U:0.}

Nsed=n_elements((*(*!dustem_data).sed).filt_names)
str_input_SED=replicate(one_str_input_SED,Nsed)
str_input_SED.FILTER=(*(*!dustem_data).sed).filt_names
str_input_SED.Wavelength=(*(*!dustem_data).sed).wav
str_input_SED.I=(*(*!dustem_data).sed).values
str_input_SED.varianceII=(*(*!dustem_data).sed).sigma   ;should it be squared ??

N_predicted_sed=n_elements(dustem_predicted_sed)
str_predicted_SED=replicate(one_str_predicted_SED,N_predicted_sed)
str_predicted_EXT=replicate(one_str_predicted_SED,N_predicted_sed)
fully_undefined_sed=str_predicted_SED
str_predicted_SED.FILTER=(*(*!dustem_data).sed).filt_names    ;taken from input SED
str_predicted_SED.Wavelength=(*(*!dustem_data).sed).wav       ;taken from input SED
str_predicted_SED.I=dustem_predicted_sed    
;stop
IF !run_pol THEN BEGIN
	str_predicted_SED.Q=dustem_predicted_Qsed
	str_predicted_SED.U=dustem_predicted_Used
	;str_predicted_EXT.Q=dustem_predicted_Qext
	;str_predicted_EXT.U=dustem_predicted_Uext
ENDIF ELSE BEGIN
	str_predicted_SED.Q=fully_undefined_sed.Q
	str_predicted_SED.U=fully_undefined_sed.U
	;str_predicted_EXT.Q=fully_undefined_sed.Q
	;str_predicted_EXT.U=fully_undefined_sed.U
ENDELSE

N_predicted_spectra=n_elements(dustem_spectra_st.sed)

;===== 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
parameters_func_values=*(*!dustem_fit).PARAM_FUNC

Nparams=n_elements(parameters_desc)

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

;stop
;This is what will be saved
;help,str_input_SED,/str
;help,str_predicted_SED,/str
;help,dustem_spectra_st.sed,/str
;help,dustem_spectra_st.ext,/str

file='./dustemwrap_tmp.fits'
;==== Write the observed SED in extension 1
unit=1L
mwrfits,str_input_SED,file,/create
unit_input_sed=unit
unit=unit+1
;==== add the predicted SED in extension 2
mwrfits,str_predicted_SED,file
unit_predicted_sed=unit
unit=unit+1
;==== add the predicted per-grain spectrum in extension 3
mwrfits,dustem_spectra_st.sed,file
unit_predicted_spectra=unit
unit=unit+1
;==== add the predicted per-grain extinction in extension 4
mwrfits,dustem_spectra_st.ext,file
unit_predicted_extinctions=unit

;print,unit_input_sed,unit_predicted_sed,unit_predicted_spectra,unit_predicted_extinctions
;stop

;===== read back extensions to get proper extension headers
file='./dustemwrap_tmp.fits'
sst_input_sed=mrdfits(file,unit_input_sed,header_input_sed)
sst_predicted_SED=mrdfits(file,unit_predicted_sed,header_predicted_SED)
;sst_predicted_spectrum_tot=mrdfits(file,3,header_predicted_spectrum_tot)
dustem_spectra_st.sed=mrdfits(file,unit_predicted_spectra,header_predicted_spectra)
dustem_spectra_st.ext=mrdfits(file,unit_predicted_extinctions,header_predicted_extinctions)

;===== complement extension header for input SED
sxaddpar,header_input_sed,'TITLE','OBSERVED INPUT SED TO DUSTEMWRAP',before='COMMENT'
sxaddpar,header_input_sed,'TTYPE1','FILTER','Dustemwrap Filter name'
sxaddpar,header_input_sed,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
sxaddpar,header_input_sed,'TTYPE3','I','Observed Total Stokes I Intensity [??]'
sxaddpar,header_input_sed,'TTYPE4','Q','Observed Stokes Q Intensity [??]'
sxaddpar,header_input_sed,'TTYPE5','U','Observed Stokes U Intensity [??]'
sxaddpar,header_input_sed,'TTYPE6','VARIANCEII','Observed Total Stokes I Variance [??]'
sxaddpar,header_input_sed,'TTYPE7','VARIANCEQQ','Observed Stokes Q Variance [??]'
sxaddpar,header_input_sed,'TTYPE8','VARIANCEUU','Observed Stokes U Variance [??]'

;===== complement extension header for predicted SED
sxaddpar,header_predicted_SED,'TITLE','OUTPUT BEST FIT SED BY DUSTEMWRAP'
sxaddpar,header_predicted_SED,'TTYPE1','FILTER','Dustemwrap Filter name'
sxaddpar,header_predicted_SED,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
sxaddpar,header_predicted_SED,'TTYPE3','I','Observed Total Stokes I Intensity [??]'
sxaddpar,header_predicted_SED,'TTYPE4','Q','Observed Stokes Q Intensity [??]'
sxaddpar,header_predicted_SED,'TTYPE5','U','Observed Stokes U Intensity [??]'
sxaddpar,header_predicted_SED,'CHI2',best_fit_chi2,'Dustemwrap best fit chi2 [-]'
sxaddpar,header_predicted_SED,'RCHI2',best_fit_rchi2,'Dustemwrap best fit reduced chi2 [-]'
sxaddpar,header_predicted_SED,'NGRAINS',Ngrains,'Number of grain species'
;====== write dustem wrap parameter names
sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter names ======',after='RCHI2'
sxaddpar,header_predicted_SED,'NPARAMS',Nparams,' Number of DustemWrap free parameters'
FOR i=0L,Nparams-1 DO BEGIN
	sxaddpar,header_predicted_SED,'PARN'+strtrim(i+1,2),parameters_desc[i],' DustemWrap parameter name'
ENDFOR
sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARN'+strtrim(i,2)

;====== write dustem wrap parameter best fit values
sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap dust model used ======',after='COMMENT'
sxaddpar,header_predicted_SED,'MODEL',used_model,'Dust model used'
sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='MODEL'

sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap dust model used in polarisation ? ======',after='COMMENT'
sxaddpar,header_predicted_SED,'POL',!run_pol,'1 if pol was used, 0 otherwise'
sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='POL'

;====== write dustem wrap parameter best fit values
sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter best fit values ======',after='COMMENT'
FOR i=0L,Nparams-1 DO BEGIN
	sxaddpar,header_predicted_SED,'PARV'+strtrim(i+1,2),best_parameters[i],'Param. '+parameters_desc[i]
ENDFOR
sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARV'+strtrim(i,2)

;====== write dustem wrap parameter best fit values uncertainties
sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap best fit parameter uncertainties ======',after='COMMENT'
FOR i=0L,Nparams-1 DO BEGIN
	sxaddpar,header_predicted_SED,'PARU'+strtrim(i+1,2),best_parameters_uncertainties[i],'Param. '+parameters_desc[i]
ENDFOR
sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARU'+strtrim(i,2)

;====== write dustem wrap parameter initial fit values
sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter initial fit values ======',after='COMMENT'
FOR i=0L,Nparams-1 DO BEGIN
	sxaddpar,header_predicted_SED,'PARI'+strtrim(i+1,2),parameters_init_values[i],'Param. '+parameters_desc[i]
ENDFOR
sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARI'+strtrim(i,2)

;====== write dustem wrap parameter func values (not sure exactly what this is)
sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter func values ======',after='COMMENT'
FOR i=0L,Nparams-1 DO BEGIN
	sxaddpar,header_predicted_SED,'PARF'+strtrim(i+1,2),parameters_func_values[i],'Param. '+parameters_desc[i]
ENDFOR
sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARF'+strtrim(i,2)

;===== complement extension header for predicted spectrum total
;sxaddpar,header_predicted_spectrum_tot,'TITLE','OUTPUT BEST FIT TOTAL SPECTRUM BY DUSTEMWRAP'
sxaddpar,header_predicted_spectra,'TITLE','OUTPUT BEST FIT PER-GRAIN AND TOTAL EMISSION SPECTRA BY DUSTEMWRAP'
sxaddpar,header_predicted_extinctions,'TITLE','OUTPUT BEST FIT PER-GRAIN AND TOTAL EXTINCTION BY DUSTEMWRAP'

;FOR i=0L,Ngrains-1 DO BEGIN
;	hh=*headers_predicted_spectrum_grain[i]
;	sxaddpar,hh,'TITLE','OUTPUT BEST FIT SPECTRUM BY DUSTEMWRAP FOR GRAIN '+strtrim(i+1,2)
;	*headers_predicted_spectrum_grain[i]=hh
;ENDFOR

;stop
;hprint,header_input_sed
;hprint,header_predicted_SED

;write final file
IF keyword_set(filename) THEN BEGIN
	file=filename
ENDIF ELSE BEGIN
	file='./dustemwrap_results.fits'
ENDELSE

;==== Write the observed SED in extension 1
mwrfits,str_input_SED,file,header_input_sed,/create
unit=unit+1
;==== add the predicted SED in extension 2
mwrfits,str_predicted_SED,file,header_predicted_SED
;==== add the predicted total spectrum in extension 3
;mwrfits,str_predicted_spectrum_tot,file,header_predicted_spectrum_tot
;==== add the predicted spectrum for each grain type in extension 4 - 4+Ngrains
;FOR i=0L,Ngrains-1 DO BEGIN
;	mwrfits,*str_predicted_spectrum_grains[i],file,*headers_predicted_spectrum_grain[i]
;ENDFOR
;==== add the predicted per-grain spectrum in extension 3
;stop
mwrfits,dustem_spectra_st.sed,file,header_predicted_spectra
;==== add the predicted per-grain extinction in extension 4
mwrfits,dustem_spectra_st.ext,file,header_predicted_extinctions


message,'Wrote '+file,/info

;stop
;hprint,header_input_sed
;hprint,header_predicted_SED


the_end:

END