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:

;    filename   = 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.

;    All dustem system variables are reinitialized

;

; 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

;stop

used_model=!dustem_model

emission_sed_present=ptr_valid((*!dustem_data).sed)
extinction_sed_present=ptr_valid((*!dustem_data).ext)
;dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st  
;fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7 ;this is correct despite the way it is presented (1.e4/1e.7*1e.20=1.e17 which is the factor to convert from ergs to Mega Janskys)
;SED_spec = st.sed.em_tot * fact

;stop

;===== run dustem to get predicted spectrum and sed

IF emission_sed_present THEN BEGIN
	dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,st=dustem_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]
	ENDIF

    str_predicted_SED=dustem_make_fits_predicted_sed(*!dustem_data,dustem_predicted_sed $
    												,str_input_sed=str_input_sed $
    												,dustem_predicted_Qsed=dustem_predicted_Qsed $
    												,dustem_predicted_Used=dustem_predicted_Used)
    str_predicted_shown_SED=dustem_make_fits_predicted_sed(*!dustem_show,dustem_predicted_sed $
    												,str_input_sed=str_shown_sed $
    												,dustem_predicted_Qsed=dustem_predicted_Qsed $
    												,dustem_predicted_Used=dustem_predicted_Used)

ENDIF

IF extinction_sed_present THEN BEGIN
	;stop
	dustem_predicted_ext=dustem_compute_ext(*(*!dustem_fit).current_param_values,st=dustem_st)
	IF !run_pol THEN BEGIN
	;This is needed to compute Q,U of both SED and spectra
	   dustem_predicted_polext=dustem_compute_stokext((*!dustem_fit).current_param_values,st=dustem_st)

	   dustem_predicted_Qext=dustem_predicted_polext[0]
	   dustem_predicted_Uext=dustem_predicted_polext[1]
	ENDIF

    str_predicted_EXT=dustem_make_fits_predicted_ext(*!dustem_data,dustem_predicted_ext $
    												,dustem_predicted_Qext=dustem_predicted_Qext $
    												,dustem_predicted_Uext=dustem_predicted_Uext $
    												,str_input_ext=str_input_ext)
    str_predicted_shown_EXT=dustem_make_fits_predicted_ext(*!dustem_show,dustem_predicted_ext $
    												,dustem_predicted_Qext=dustem_predicted_Qext $
    												,dustem_predicted_Uext=dustem_predicted_Uext $    	
    												,str_input_ext=str_shown_ext)

ENDIF

Ngrains=n_tags(dustem_st.sed)-2

;===== 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

IF ptr_valid((*!dustem_fit).fixed_param_descs) THEN BEGIN
	there_are_fixed_params=1
	fixed_parameters_desc=*(*!dustem_fit).fixed_param_descs
	fixed_parameters_values=*(*!dustem_fit).fixed_param_init_values
	Nfixedparams=n_elements(fixed_parameters_desc)
ENDIF ELSE BEGIN
	there_are_fixed_params=0
	Nfixedparams=0
ENDELSE

Nparams=n_elements(parameters_desc)

file='./dustemwrap_tmp.fits'

;==== Write the observed SED in extension 1

unit=1L

;Yes, toto does not exist. This is just to write a dummy input to allow for a general header
mwrfits,toto,file,/create
IF emission_sed_present THEN BEGIN
	;==== add the input emission SED
	mwrfits,str_input_SED,file
	unit_input_sed=unit
	unit=unit+1
	;==== add the predicted emission SED
	mwrfits,str_predicted_SED,file
	unit_predicted_sed=unit
	unit=unit+1

	;==== add the shown emission SED
	mwrfits,str_shown_SED,file
	unit_shown_sed=unit
	unit=unit+1

ENDIF ELSE BEGIN
	unit_input_sed=-1
	unit_predicted_sed=-1

	unit_shown_sed=-1

ENDELSE
IF extinction_sed_present THEN BEGIN
	;stop
	mwrfits,str_input_EXT,file
	unit_input_ext=unit
	unit=unit+1

	;==== add the predicted extinction SED

	mwrfits,str_predicted_EXT,file
	unit_predicted_ext=unit
	unit=unit+1

	;==== add the shown extinction SED
	mwrfits,str_shown_EXT,file
	unit_shown_ext=unit
	unit=unit+1

ENDIF ELSE BEGIN
	unit_input_ext=-1
	unit_predicted_ext=-1

	unit_shown_ext=-1

ENDELSE
;==== add the predicted per-grain emission spectrum
mwrfits,dustem_st.sed,file

unit_predicted_spectra=unit
unit=unit+1

;==== add the predicted per-grain extinction spectrum
mwrfits,dustem_st.ext,file

unit_predicted_extinctions=unit

unit=unit+1
IF !run_pol EQ 1 THEN BEGIN
	mwrfits,dustem_st.polsed,file
	unit_predicted_spectra_pol=unit
	unit=unit+1
	mwrfits,dustem_st.polext,file
	unit_predicted_extinctions_pol=unit
	unit=unit+1
ENDIF ELSE BEGIN
	unit_predicted_spectra_pol=-1
	unit_predicted_extinctions_pol=-1
ENDELSE

;==== add the pluggins spectra
plugin_tags=tag_names(*!dustem_plugin)
ind=where(plugin_tags NE 'NONE',Nplugins)
IF Nplugins NE 0 THEN BEGIN
	unit_plugins=lonarr(Nplugins)
	str_plugins=strarr(Nplugins)
	FOR i=0L,Nplugins-1 DO BEGIN
		mwrfits,*((*!dustem_plugin).(i).spec),file
		unit_plugins[i]=unit
		str_plugins[i]=plugin_tags[i]
		unit=unit+1
	ENDFOR

ENDIF ELSE BEGIN  ;This is the case where there are no plugins involved
	unit_plugins=[-1]
	str_plugins=['']

ENDELSE

extension_numbers=[unit_input_sed,unit_predicted_sed,unit_shown_sed, $
				   unit_input_ext,unit_predicted_ext,unit_shown_ext, $

				   unit_predicted_spectra,unit_predicted_extinctions, $
				   unit_predicted_spectra_pol,unit_predicted_extinctions_pol]

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']

extension_numbers=[extension_numbers,unit_plugins]
extensions_strs=[extensions_strs,'PLUGIN_'+str_plugins]

ind=where(extension_numbers NE -1,count)
extension_numbers=extension_numbers[ind]
extensions_strs=extensions_strs[ind]
order=sort(extension_numbers)
extension_numbers=extension_numbers[order]
extensions_strs=extensions_strs[order]
Nextensions=n_elements(extension_numbers)

;===== read back fits file to get proper extension headers

file='./dustemwrap_tmp.fits'

toto=mrdfits(file,0,hdr)      ;this would be to get the main header

IF unit_input_sed NE -1 THEN sst_input_sed=mrdfits(file,unit_input_sed,header_input_sed)
IF unit_predicted_sed NE -1 THEN sst_predicted_SED=mrdfits(file,unit_predicted_sed,header_predicted_SED)

IF unit_shown_sed NE -1 THEN sst_shown_SED=mrdfits(file,unit_shown_sed,header_shown_SED)

IF unit_input_ext NE -1 THEN sst_input_ext=mrdfits(file,unit_input_ext,header_input_ext)
IF unit_predicted_ext NE -1 THEN sst_predicted_ext=mrdfits(file,unit_predicted_ext,header_predicted_ext)

IF unit_shown_ext NE -1 THEN sst_shown_ext=mrdfits(file,unit_shown_ext,header_shown_ext)

dustem_st.sed=mrdfits(file,unit_predicted_spectra,header_predicted_spectra)
dustem_st.ext=mrdfits(file,unit_predicted_extinctions,header_predicted_extinctions)
IF !run_pol EQ 1 THEN BEGIN
	dustem_st.polsed=mrdfits(file,unit_predicted_spectra_pol,header_predicted_spectra)
	dustem_st.polext=mrdfits(file,unit_predicted_extinctions_pol,header_predicted_extinctions)
ENDIF

IF Nplugins NE 0 THEN BEGIN
	header_plugins=ptrarr(Nplugins)
	FOR i=0L,Nplugins-1 DO BEGIN
		a=mrdfits(file,unit_plugins[i],header_plugin)
		header_plugins[i]=ptr_new(header_plugin)
	ENDFOR
ENDIF

;===== write content of extensions in main header
FOR i=0L,Nextensions-1 DO BEGIN
	sxaddpar,hdr,'EXT'+strtrim(i+1,2),extensions_strs[i],'fits file extension'
ENDFOR

;====== write dustem wrap parameter best fit values
sxaddpar,hdr,'COMMENT','======= Dustemwrap dust model used ======',after='COMMENT'
sxaddpar,hdr,'MODEL',used_model,'Dust model used'
sxaddpar,hdr,'WRAP_V',!dustem_version.version,'Dustwrap version used'
sxaddpar,hdr,'COMMENT','===================================================',after='MODEL'
sxaddpar,hdr,'CHI2',best_fit_chi2,'Dustemwrap best fit chi2 [-]'
sxaddpar,hdr,'RCHI2',best_fit_rchi2,'Dustemwrap best fit reduced chi2 [-]'
sxaddpar,hdr,'NGRAINS',Ngrains,'Number of grain species'

sxaddpar,hdr,'COMMENT','======= Dustemwrap plugins names ======',after='NGRAINS'
sxaddpar,hdr,'NPLUG',Nplugins,'Number of plugins involved'
FOR i=0L,Nplugins-1 DO BEGIN
	sxaddpar,hdr,'PLUG'+strtrim(i+1,2),str_plugins[i],'PLUGIN name'
ENDFOR

;====== write dustem wrap parameter names

sxaddpar,hdr,'COMMENT','======= Dustemwrap parameter names ======',after='PLUG'+strtrim(i+1,2)

sxaddpar,hdr,'NPARAMS',Nparams,' Number of DustemWrap free parameters'

FOR i=0L,Nparams-1 DO BEGIN

	sxaddpar,hdr,'PARN'+strtrim(i+1,2),parameters_desc[i],' DustemWrap parameter name'

ENDFOR

sxaddpar,hdr,'COMMENT','===================================================',after='PARN'+strtrim(i,2)
sxaddpar,hdr,'COMMENT','======= Dustemwrap dust model used in polarisation ? ======',after='COMMENT'
sxaddpar,hdr,'POL',!run_pol,'1 if pol was used, 0 otherwise'
sxaddpar,hdr,'COMMENT','===================================================',after='POL'

;====== write dustem wrap parameter best fit values

sxaddpar,hdr,'COMMENT','======= Dustemwrap parameter best fit values ======',after='COMMENT'

FOR i=0L,Nparams-1 DO BEGIN

	sxaddpar,hdr,'PARV'+strtrim(i+1,2),best_parameters[i],'Param. '+parameters_desc[i]

ENDFOR

sxaddpar,hdr,'COMMENT','===================================================',after='PARV'+strtrim(i,2)

;====== write dustem wrap parameter best fit values uncertainties

sxaddpar,hdr,'COMMENT','======= Dustemwrap best fit parameter uncertainties ======',after='COMMENT'

FOR i=0L,Nparams-1 DO BEGIN

	sxaddpar,hdr,'PARU'+strtrim(i+1,2),best_parameters_uncertainties[i],'Param. '+parameters_desc[i]

ENDFOR

sxaddpar,hdr,'COMMENT','===================================================',after='PARU'+strtrim(i,2)

;====== write dustem wrap parameter initial fit values

sxaddpar,hdr,'COMMENT','======= Dustemwrap parameter initial fit values ======',after='COMMENT'

FOR i=0L,Nparams-1 DO BEGIN

	sxaddpar,hdr,'PARI'+strtrim(i+1,2),parameters_init_values[i],'Param. '+parameters_desc[i]

ENDFOR

sxaddpar,hdr,'COMMENT','===================================================',after='PARI'+strtrim(i,2)

;====== write dustem wrap parameter func values (not sure exactly what this is)

sxaddpar,hdr,'COMMENT','======= Dustemwrap parameter func values ======',after='COMMENT'

FOR i=0L,Nparams-1 DO BEGIN

	sxaddpar,hdr,'PARF'+strtrim(i+1,2),parameters_func_values[i],'Param. '+parameters_desc[i]

ENDFOR

sxaddpar,hdr,'COMMENT','===================================================',after='PARF'+strtrim(i,2)

sxaddpar,hdr,'NFPARAMS',Nfixedparams,' Number of DustemWrap fixed parameters'
FOR i=0L,Nfixedparams-1 DO BEGIN
	sxaddpar,hdr,'FPARN'+strtrim(i+1,2),fixed_parameters_desc[i],' Fixed DustemWrap parameter name'
ENDFOR
sxaddpar,hdr,'COMMENT','======= Dustemwrap fixed parameter values ======',after='FPARN'+strtrim(i,2)
FOR i=0L,Nfixedparams-1 DO BEGIN
	sxaddpar,hdr,'FPARV'+strtrim(i+1,2),fixed_parameters_values[i],'Param. '+fixed_parameters_desc[i]
ENDFOR
sxaddpar,hdr,'COMMENT','===================================================',after='FPARV'+strtrim(i,2)

IF unit_input_sed NE -1 THEN BEGIN

    ;===== 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','STOKESI','Observed Total Stokes I Intensity [??]'
    sxaddpar,header_input_sed,'TTYPE4','STOKESQ','Observed Stokes Q Intensity [??]'
    sxaddpar,header_input_sed,'TTYPE5','STOKESU','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 [??]'

    sxaddpar,header_input_sed,'NH_SED',*!dustem_hcd,'Column density [H/cm2]'
;    sxaddpar,header_input_sed,'NH',*!dustem_hcd,'Column density [H/cm2]'

    
    ;===== 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','STOKESI','Observed Total Stokes I Intensity [??]'
    sxaddpar,header_predicted_SED,'TTYPE4','STOKESQ','Observed Stokes Q Intensity [??]'
    sxaddpar,header_predicted_SED,'TTYPE5','STOKESU','Observed Stokes U Intensity [??]'

;ENDIF ELSE BEGIN;IC: hesitated between this and 'IF unit_input_ext NE -1 THEN BEGIN'. This is under the assumtion that we always fit something.

 END
IF unit_input_ext NE -1 THEN BEGIN

    ;===== complement extension header for input extinction SED
    sxaddpar,header_input_ext,'TITLE','OBSERVED INPUT EXTINCTION SED TO DUSTEMWRAP',before='COMMENT'
    sxaddpar,header_input_ext,'TTYPE1','FILTER','Dustemwrap Filter name'
    sxaddpar,header_input_ext,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
    sxaddpar,header_input_ext,'TTYPE3','STOKESI_EXT','Observed Total Stokes I optical depth []'
    sxaddpar,header_input_ext,'TTYPE4','STOKESQ_EXT','Observed Stokes Q optical depth []'
    sxaddpar,header_input_ext,'TTYPE5','STOKESU_EXT','Observed Stokes U optical depth []'
    sxaddpar,header_input_ext,'TTYPE6','VARIANCEII_EXT','Observed Total Stokes I optical depth Variance [??]'
    sxaddpar,header_input_ext,'TTYPE7','VARIANCEQQ_EXT','Observed Stokes Q optical depth Variance [??]'
    sxaddpar,header_input_ext,'TTYPE8','VARIANCEUU_EXT','Observed Stokes U optical depth Variance [??]'

    sxaddpar,header_input_ext,'NH_EXT',*!dustem_hcd,'Column density [H/cm2]'
;    sxaddpar,header_input_ext,'NH',*!dustem_hcd,'Column density [H/cm2]'

    
    ;===== complement extension header for predicted extinction SED
    sxaddpar,header_predicted_EXT,'TITLE','OUTPUT BEST FIT EXTINCTION SED BY DUSTEMWRAP'
    sxaddpar,header_predicted_EXT,'TTYPE1','FILTER','Dustemwrap Filter name'
    sxaddpar,header_predicted_EXT,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
    sxaddpar,header_predicted_EXT,'TTYPE3','STOKESI_EXT','Observed Total Stokes I optical depth []'
    sxaddpar,header_predicted_EXT,'TTYPE4','STOKESQ_EXT','Observed Stokes Q optical depth []'
    sxaddpar,header_predicted_EXT,'TTYPE5','STOKESU_EXT','Observed Stokes U optical depth []'

END

;===== complement extension header for predicted spectrum total

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'

;===== complement plugin headers
FOR i=0L,Nplugins-1 DO BEGIN
	header=*header_plugins[i]
	sxaddpar,header,'PLUGIN',str_plugins[i],'Plugin name'
	sxaddpar,header,'SCOPE',*((*!dustem_plugin).(i).scope),'Plugin scope'
	tags=*((*!dustem_plugin).(i).paramtag)
	FOR j=0L,n_elements(tags)-1 DO BEGIN
		sxaddpar,header,'TAGN'+strtrim(j+1,2),(*((*!dustem_plugin).(i).paramtag))[j],'Plugin parameter tag'
	ENDFOR
	header_plugins[i]=ptr_new(header)
ENDFOR

;hprint,*header_plugins[0]
;stop

;============== write final fits file

IF keyword_set(filename) THEN BEGIN
	file=filename
ENDIF ELSE BEGIN

	file='./dustemwrap_results.fits'

ENDELSE

;==== Write the observed SED (if any)
mwrfits,0,file,hdr,create=1
IF unit_input_sed NE -1 THEN BEGIN

	mwrfits,str_input_SED,file,header_input_SED

ENDIF
;==== add the predicted SED (if any)

IF unit_predicted_sed NE -1 THEN BEGIN

	mwrfits,str_predicted_SED,file,header_predicted_SED
ENDIF

;==== add the shown SED (if any)
IF unit_shown_sed NE -1 THEN BEGIN
	mwrfits,str_shown_SED,file,header_shown_SED
ENDIF
;==== Write the observed extinction SED (if any)

IF unit_input_ext NE -1 THEN BEGIN

	mwrfits,str_input_EXT,file,header_input_EXT

ENDIF
;==== add the predicted SED (if any)
IF unit_predicted_ext NE -1 THEN BEGIN

	mwrfits,str_predicted_EXT,file,header_predicted_EXT
ENDIF
;==== add the show extinction SED (if any)
IF unit_shown_ext NE -1 THEN BEGIN
	mwrfits,str_shown_EXT,file,header_shown_EXT

ENDIF
;==== add the predicted per-grain emission spectra
mwrfits,dustem_st.sed,file,header_predicted_spectra
;==== add the predicted per-grain extinction
mwrfits,dustem_st.ext,file,header_predicted_extinctions
IF !run_pol EQ 1 THEN BEGIN
	mwrfits,dustem_st.polsed,file
	mwrfits,dustem_st.polext,file
ENDIF

FOR i=0L,Nplugins-1 DO BEGIN
  mwrfits,*((*!dustem_plugin).(i).spec),file,*header_plugins[i]
ENDFOR

message,'================================================',/info

message,'Wrote '+file,/info

message,'File extensions are: ',/info

FOR i=0L,Nextensions-1 DO BEGIN

	message,'extension '+strtrim(extension_numbers[i],2)+' : '+strtrim(extensions_strs[i],2),/info

ENDFOR
message,'================================================',/info

;stop

the_end:

END