PRO dustem_write_fits_table,filename=filename,help=help

;+
; NAME:
;    dustem_write_fits_table
; PURPOSE:
;    Write a fits table with content of current Dustem system variables
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_write_fits_table,file[,/help]
; INPUTS:
;    file      = File name
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    File is stored.
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at IRAP
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

; !DUSTEM_DAT = '/tmp/ramdisk8MB/'       
; !DUSTEM_DATA = -> <Anonymous> Array[1]
; !DUSTEM_DO_CC =        1
; !DUSTEM_F90_EXEC = '$HOME/Soft_Librairies/dustem4.3_wk/src/dustem'
; !DUSTEM_FILTERS = <PtrHeapVar570>
; !DUSTEM_FIT = <PtrHeapVar1>
; !DUSTEM_F_HI =       1.00000   
; !DUSTEM_IDL_CONTINUUM =       0.00000
; !DUSTEM_IDL_FREEFREE =       0.00000
; !DUSTEM_IDL_SYNCHROTRON =       0.00000
; !DUSTEM_INPUTS = <PtrHeapVar2>
; !DUSTEM_INSTRUMENT_DESCRIPTION = -> <Anonymous> Array[120]
; !DUSTEM_NEVER_DO_CC =        0 
; !DUSTEM_PARAMS = <PtrHeapVar15>
; !DUSTEM_PARINFO = <PtrHeapVar943>
; !DUSTEM_PREVIOUS_CC = <PtrHeapVar3346>
; !DUSTEM_RES = '/tmp/ramdisk8MB/'
; !DUSTEM_SHOW_PLOT =        1 
; !DUSTEM_SOFT_DIR = '$HOME/Soft_Librairies/dustem4.3_wk/'       
; !DUSTEM_VERBOSE =        1
; !DUSTEM_WHICH = 'WEB3p8'
; !DUSTEM_WRAP_SOFT_DIR = '$HOME/Soft_Librairies/dustem-wrapper_idl/'
; !FIT_RCHI2_WEIGHT = -> <Anonymous> Array[1]
; !RUN_ANIS =       0.00000
; !RUN_CIRC =       0.00000      
; !RUN_IONFRAC =       0.00000
; !RUN_LIN =       1.00000
; !RUN_POL =        1
; !RUN_RRF =       0.00000
; !RUN_TLS =       0.00000
; !RUN_UNIV =       0.00000

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 not keyword_set(function_name) THEN BEGIN 
  dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,out_st=dustem_spectra_st)
ENDIF ELSE BEGIN
ENDELSE

;===== define structures containing results
one_str_input_SED={FILTER:'',Wavelength:0.,I:0.,Q:0.,U:0.,varianceII:0.,varianceQQ:0.,varianceUU:0.}
one_str_predicted_SED={FILTER:'',Wavelength:0.,I:0.,Q:0.,U:0.}
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_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    ;but where is Q,U ??
;str_predicted_SED.Q=
;str_predicted_SED.U=


N_predicted_spectra=n_elements(dustem_spectra_st.sed)
str_predicted_spectrum_tot=replicate(one_str_predicted_spectrum_tot,N_predicted_spectra)
str_predicted_spectrum_tot.Wavelength=dustem_spectra_st.SED.wav
str_predicted_spectrum_tot.I=dustem_spectra_st.SED.em_tot
;str_predicted_spectrum_tot.Q=dustem_spectra_st.polsed         ;not exactly, should have Q,U
;str_predicted_spectrum_tot.U=dustem_spectra_st.polsed         ;not exactly, should have Q,U

;===== 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=n_elements(parameters_desc)

;stop
file='/tmp/dustem_wrap_ouput.fits'
;==== Write the observed SED in extension 1
mwrfits,str_input_SED,file,/create
;==== add the predicted SED in extension 2
mwrfits,str_predicted_SED,file
;==== add the predicted total spectrum in extension 3
mwrfits,str_predicted_spectrum_tot,file

; mwrfits, xinput, file, header,              $
;         ascii=ascii,                            $
;        separator=separator,                    $
;        terminator=terminator,                  $
;        create=create,                          $
;        null=null,                              $
;        group=group,                            $
;        pscale=pscale, pzero=pzero,             $
;        alias=alias,                            $
;        use_colnum = use_colnum,                $
;        lscale=lscale, iscale=iscale,              $
;        no_copy = no_copy,                      $
;        bscale=bscale,                          $
;        no_types=no_types,                      $
;        silent=silent,                          $
;        no_comment=no_comment,                  $
;        logical_cols=logical_cols,              $
;        bit_cols=bit_cols,                      $
;        nbit_cols=nbit_cols,                    $
;        status = status,                        $
;        version=version

; mrdfits, file, extension, header,      $
;         structyp = structyp,                    $
;         use_colnum = use_colnum,                $
;         range = range,                          $
;         dscale = dscale, fscale=fscale,         $
;         silent = silent,                        $
;         columns = columns,                      $
;         no_tdim = no_tdim,                      $
;         error_action = error_action,            $
;  	compress=compress,                      $
; 	alias=alias,                            $
;         rows = rows,                        $
; 	unsigned=unsigned,                      $
; 	version=version,                        $
; 	pointer_var=pointer_var,                $
; 	fixed_var=fixed_var,                    $
;         status=status, extnum = extnum

;===== read back extensions to get proper extension headers
file='/tmp/dustem_wrap_ouput.fits'
sst_input_sed=mrdfits(file,1,header_input_sed)
sst_predicted_SED=mrdfits(file,2,header_predicted_SED)
sst_predicted_spectrum_tot=mrdfits(file,3,header_predicted_spectrum_tot)

;===== 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 [-]'
;====== write dustem wrap parameter names
sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter names ======',after='RCHI2'
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 valies
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)

;===== complement extension header for predicted spectrum total
sxaddpar,header_predicted_spectrum_tot,'TITLE','OUTPUT BEST FIT TOTAL SPECTRUM BY DUSTEMWRAP'

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

;==== Write the observed SED in extension 1
mwrfits,str_input_SED,file,header_input_sed,/create
;==== 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

message,'Wrote '+file,/info

stop
hprint,header_input_sed
hprint,header_predicted_SED


the_end:

END