PRO dustem_write_sed_grid,model $
						 						 ,parameters_description $
						 						 ,iv_min $
						 						 ,iv_max $
						 						 ,iv_Nvalues $
						 						 ,plog $
						 						 ,parameter_values $
												 ,seds $
												 ,em_spectra=em_spectra $
												 ,ext_spectra=ext_spectra $
						 						 ,fpd=fpd $
						 						 ,fiv=fiv $
						 						 ,filename=filename $
						 						 ,use_polarisation=use_polarisation $
						 						 ,double=double $
						 						 ,help=help

;+
; NAME:
;    dustem_write_sed_grid
; PURPOSE:
;    Writes a model GRID fits file
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_write_sed_grid,model,parameters_description,iv_min,iv_max,iv_Nvalues,plog $
;						 										,parameter_values,seds,em_spectra=em_spectra,ext_spectra=ext_spectra $
;																,fpd=fpd,fiv=fiv,filename=filename,use_polarisation=use_polarisation
; INPUTS:
;    model                  : dustemwrap model name to be used
;    parameters_description : dustemwrap parameter description array
;    iv_min                 : minimum values for each parameter
;    iv_max                 : maximum values for each parameter
;    iv_Nvalues             : number of parameter values for each parameter
;    plog                   : array indicating if a given parameter is to be sampled in linear (0) or log scale (1)
;    parameter_values       : array of parameter values as computed with dustem_param_range2param_values.pro
;    seds                   : array of seds
; OPTIONAL INPUT PARAMETERS:
;    em_spectra             : Emission spectra
;    ext_spectra            : Extinction spectra
;    fpd                    : Fixed parameter descritions
;    fiv                    : fixed parameter values
;    filename               : output fits file name for the table (default ='/tmp/dustem_seds.fits')
;    use_polarisation       : flag indicating if the model has polarisation predictions.
;    double                 : if set, writes the table in double accuracy
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    A file is written
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
; PROCEDURE:
;    Filters in the table are ordered by increasing wavelength
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard (2023)
;    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_sed_grid'
  goto,the_end
ENDIF

polar=0
IF keyword_set(use_polarisation) THEN polar=1

IF keyword_set(em_spectra) AND keyword_set(ext_spectra) THEN BEGIN
  	include_spectra=1
ENDIF ELSE BEGIN
	include_spectra=0
	message,'either em_spectra or ext_spectra is unset',/continue
	message,'Spectra will not be saved in the Grid',/continue
	;stop
ENDELSE

;stop

Nc=n_elements(parameter_values)
Nparams=n_elements(*parameter_values[0])
use_filters=(*seds[0]).filter
Nfilters=n_elements(use_filters)

wavs=dustem_get_wavelengths()
Nwavs=n_elements(wavs)

format_template='0.E0'
double4header=0
IF keyword_set(double) THEN BEGIN
	format_template='0.D0'
	double4header=1
ENDIF

;=== do the parameter value structure
str='one_param_st={'
FOR i=0L,Nparams-1 DO BEGIN
  str=str+'P'+strtrim(i+1,2)+':'+format_template
  IF i NE Nparams-1 THEN str=str+','
ENDFOR
str=str+'}'

toto=execute(str)

;=== do the SED value structure
str='one_Iemsed_st={'
FOR i=0L,Nfilters-1 DO BEGIN
	istr=strtrim(use_filters[i],2)
	str=str+'I'+istr+':'+format_template
	str=str+','
ENDFOR
str=str+'total:0.d0}'
toto=execute(str)

;=== do the IR wavelengths structure
one_wave_st={wav:0.D0}

;=== do the Emission total spectrum structure
str='one_Iemspec_st={'
FOR i=0L,Nwavs-1 DO BEGIN
	istr=strtrim(i+1,2)
	str=str+'W'+istr+':'+format_template
	IF i NE Nwavs-1 THEN str=str+','
ENDFOR
str=str+'}'
toto=execute(str)

params_array=replicate(one_param_st,Nc)
Iem_seds_array=replicate(one_Iemsed_st,Nc)
wav_array=replicate(one_wave_st,Nwavs)
IF include_spectra THEN BEGIN
	;=== get the number of grains
	;Ngrain=n_tags(st.sed)-2
	Ngrains=(*!dustem_params).Ngrains
	Iem_spec_array=replicate(one_Iemspec_st,Nc)
	Iext_spec_array=replicate(one_Iemspec_st,Nc)
	Iem_spec_array_per_grain=replicate(one_Iemspec_st,[Nc,Ngrains])
	Iext_spec_array_per_grain=replicate(one_Iemspec_st,[Nc,Ngrains])
ENDIF

;=== Fill up parameters array
FOR i=0L,Nc-1 DO BEGIN
	FOR j=0L,Nparams-1 DO BEGIN
		params_array[i].(j)=(*parameter_values[i])[j]
	ENDFOR
ENDFOR

;=== Fill up Iem SED array
FOR i=0L,Nc-1 DO BEGIN
	FOR j=0L,Nfilters-1 DO BEGIN
		Iem_seds_array[i].(j)=((*seds[i]).stokesI)[j]
	ENDFOR
ENDFOR
;=== compute the sum of the Iem SED array
FOR i=0L,Nc-1 DO BEGIN
	 Iem_seds_array[i].total=total((*seds[i]).stokesI)
ENDFOR

;=== fill up the wavelengths array
;wav_array=wavs
;FOR j=0L,Nwavs-1 DO BEGIN
;		wav_array.wav=wavs[j]
;ENDFOR
wav_array.wav=wavs
fact_em=dustem_get_emission_conversion_factor()

;=== Fill up Iem spec array
;stop
IF include_spectra THEN BEGIN
	FOR i=0L,Nc-1 DO BEGIN
		FOR j=0L,Nwavs-1 DO BEGIN
			Iem_spec_array[i].(j)=((*em_spectra[i]).em_tot)[j]*fact_em[j]
		ENDFOR
	ENDFOR
	;=== Fill up Iext spec array
	FOR i=0L,Nc-1 DO BEGIN
		FOR j=0L,Nwavs-1 DO BEGIN
			Iext_spec_array[i].(j)=((*ext_spectra[i]).ext_tot)[j]
		ENDFOR
	ENDFOR
	;=== Fill up Iem spec array per grain size
	FOR k=0L,Ngrains-1 DO BEGIN
		FOR i=0L,Nc-1 DO BEGIN
			FOR j=0L,Nwavs-1 DO BEGIN
				Iem_spec_array_per_grain[i,k].(j)=((*em_spectra[i]).(k+1))[j]*fact_em[j]
			ENDFOR
		ENDFOR
	ENDFOR
	;stop
	;=== Fill up Iext spec array per grain size
	FOR k=0L,Ngrains-1 DO BEGIN
		FOR i=0L,Nc-1 DO BEGIN
			FOR j=0L,Nwavs-1 DO BEGIN
				Iext_spec_array_per_grain[i,k].(j)=((*ext_spectra[i]).abs_grain)[k]+((*ext_spectra[i]).sca_grain)[k]
			ENDFOR
		ENDFOR
	ENDFOR
ENDIF
;stop

;======== save SED tables

fits_file='/tmp/dustem_seds.fits'
IF keyword_set(filename) THEN fits_file=filename
bits=str_sep(fits_file,'/')
filename_only=bits(n_elements(bits)-1)

unit_content=['Nothing']
mwrfits,params_array,fits_file,h,/create ;ext1
unit_content=[unit_content,'Model parameter values']
unit_count=1L
mwrfits,Iem_seds_array,fits_file
unit_content=[unit_content,'Emission SED values [MJy/sr for NH=1e20 H/cm2]']
unit_count=unit_count+1    				;ext1
IF include_spectra THEN BEGIN
	mwrfits,wav_array,fits_file
	unit_count=unit_count+1            				;ext2
	unit_content=[unit_content,'Spectra Wavelengths [mic]']
	mwrfits,Iem_spec_array,fits_file
	unit_count=unit_count+1       				;ext3
	unit_content=[unit_content,'Total Emission Spectrum [MJy/sr for NH=1e20 H/cm2]']
	mwrfits,Iext_spec_array,fits_file
	unit_count=unit_count+1      				;ext4
	unit_content=[unit_content,'Total Extinction Spectrum [opt depth for NH=1e20 H/cm2]']
	FOR k=0L,Ngrains-1 DO BEGIN
		mwrfits,Iem_spec_array_per_grain[*,k],fits_file
		unit_count=unit_count+1
		unit_content=[unit_content,'Grain'+strtrim(k+1,2)+' Emission Spectrum [MJy/sr for NH=1e20 H/cm2]']
	ENDFOR
	FOR k=0L,Ngrains-1 DO BEGIN
		mwrfits,Iext_spec_array_per_grain[*,k],fits_file
		unit_count=unit_count+1
		unit_content=[unit_content,'Grain'+strtrim(k+1,2)+' Extinction Spectrum [opt depth for NH=1e20 H/cm2]']
	ENDFOR
ENDIF

;==== read back the file to get main header and various extensions
all_sts=ptrarr(unit_count)
all_headers=ptrarr(unit_count)

st0=mrdfits(fits_file,0,h0)
FOR i=0L,unit_count-1 DO BEGIN
	st=mrdfits(fits_file,i+1,h)
	all_sts[i]=ptr_new(st)
	all_headers[i]=ptr_new(h)
ENDFOR

;stop
NFparams=n_elements(fpd)

;=== Modify the main header
sxaddpar,h0,'COMMENT','***** This file is a DustemWrap Grid table ****',after='EXTEND'
sxaddpar,h0,'FILE',filename_only,'Original File Name',before='END'
sxaddpar,h0,'DATE',systime(0),'Creation Date',before='END'
sxaddpar,h0,'AUTHOR',strtrim(getenv('USER'),2),'Who made that file',before='END'
sxaddpar,h0,'SOFT',!dustem_version.version,'Dustem Sofware Version',before='END'
sxaddpar,h0,'INCSP',include_spectra,'Does the grid include Emission/Extinction spectra',before='END'
sxaddpar,h0,'COMMENT','***** File unit description ****',after='SOFT'
;stop
FOR k=1L,unit_count DO BEGIN
	tag='Unit'+strtrim(k,2)
	sxaddpar,h0,tag,unit_content[k],before='END'
	last_tag=tag
ENDFOR
sxaddpar,h0,'COMMENT','***** Model description ****',after=last_tag
sxaddpar,h0,'MODEL',model,'dustemwrap model used',before='END'
sxaddpar,h0,'NGRAIN',Ngrains,'Number of grain types',before='END'
sxaddpar,h0,'NHREF',*!dustem_HCD,'dustemwrap reference NH [H/cm2]',after='MODEL'
sxaddpar,h0,'POLAR',polar,'dustemwrap polarization ?',after='NHREF'
sxaddpar,h0,'DOUBLE',double4header,'double accuracy ?',after='POLAR'
sxaddpar,h0,'COMMENT','**** Variable Dustemwrap PARAMETERS Used in the Grid ***',after='POLAR'
sxaddpar,h0,'NPARAM',Nparams,'Number of free parameters',before='END'
sxaddpar,h0,'NCOMB',Nc,'Number of parameter values combinations',before='END'
;stop
FOR j=0L,Nparams-1 DO BEGIN
	;value=sxpar(*all_headers[0],'TTYPE'+strtrim(j+1,2))
	;sxaddpar,h0,'TTYPE'+strtrim(j+1,2),value,'dustemwrap '+parameters_description[j],before='END'
	sxaddpar,h0,'PNAME'+strtrim(j+1,2),parameters_description[j],'dustemwrap param name',before='END'
	sxaddpar,h0,'PMIN'+strtrim(j+1,2),iv_min[j],'minimum param value',before='END'
	sxaddpar,h0,'PMAX'+strtrim(j+1,2),iv_max[j],'maximum param value',before='END'
	sxaddpar,h0,'PNVAL'+strtrim(j+1,2),iv_Nvalues[j],'Number of param values',before='END'
	sxaddpar,h0,'PLOG'+strtrim(j+1,2),plog[j],'1 if log scale sampling',before='END'
	last_tag='PLOG'+strtrim(j+1,2)
ENDFOR
;stop
IF keyword_set(fpd) THEN BEGIN
	sxaddpar,h0,'COMMENT','**** Fixed Dustemwrap PARAMETERS Used in the Grid ***',after='last_tag'
	sxaddpar,h0,'NFPARAM',NFparams,'Number of fixed parameters',before='END'
	FOR j=0L,NFparams-1 DO BEGIN
		sxaddpar,h0,'FPNAME'+strtrim(j+1,2),fpd[j],'fixed dustemwrap param name',before='END'
		sxaddpar,h0,'FPVALUE'+strtrim(j+1,2),fiv[j],'fixed dustemwrap param value',before='END'
		last_tag='FPVALUE'+strtrim(j+1,2)
	ENDFOR
ENDIF
sxaddpar,h0,'COMMENT','**** Dustemwrap FILTERS ***',after=last_tag
FOR j=0L,Nfilters-1 DO BEGIN
	;value=sxpar(*all_headers[2],'TTYPE'+strtrim(j+1,2))
	;sxaddpar,h0,'TTYPE'+strtrim(j+1,2),value,'Intensity in '+use_filters[j]
	sxaddpar,h0,'FNAME'+strtrim(j+1,2),use_filters[j],'dustemwrap filter name',before='END'
ENDFOR

;stop
mwrfits,bidon,fits_file,h0,/create
;FOR i=1L,unit_count-1 DO BEGIN
FOR i=0L,unit_count-1 DO BEGIN
	mwrfits,*all_sts[i],fits_file,*all_headers[i]
ENDFOR

message,'Wrote '+fits_file,/continue
;stop
;hprint,h0

the_end:

END