PRO dustem_make_sed_table,model, $
											    parameters_description, $
											    iv_min, $
											    iv_max, $
											    iv_Nvalues, $
						  						fpd=fpd, $
						  						fiv=fiv, $
						  						filename=filename, $
						  						filters=filters, $
						  						plog=plog, $
						  						show_seds=show_seds, $
						  						print_params=print_params, $
						  						help=help

;+
; NAME:
;    dustem_make_sed_table
; PURPOSE:
;    makes an SED lookup table using dustemwrap for a given dust model and free parameters
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_make_sed_table,model,parameters_description,parvalues_min,parvalues_max,parvalues_Nvalues, $
;                     	   [,fpd=][,fiv=][,filename=][,filters=][,log=][,/show_seds]
; INPUTS:
;    model                  : dustemwrap model name to be used
;    parameters_description : dustemwrap parameter description array
;    parvalues_min          : minimum values for each parameter
;    parvalues_max          : maximum values for each parameter
;    parvalues_Nvalues      : number of parameter values for each parameter
; OPTIONAL INPUT PARAMETERS:
;    filters                : name of dustemwrap filters to be included in the SED calculations (default = ALL known filters)
;    fpd                    : fixed parameter description
;    fiv                    : fixed parameter values
;    filename               : output fits file name for the table (default ='/tmp/dustem_seds.fits')
;    plog                   : array indicating if a given parameter is to be sampled in linear (0) or log scale (1)
;    show_seds              : if set, plots SEDs as they are computed.
; 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
;    dustem_make_sed_table,'DBP90',['(*!dustem_params).G0','(*!dustem_params).grains(0).mdust_o_mh'], $
;																[0.1,1.e-3],[100.,1.e-1],[10,3],filename='/tmp/IRAS_GO.fits',plog=[1,0],/show_seds, $
;                              	filters=[dustem_instru2filters('IRAC'),dustem_instru2filters('SPIRE')]
; 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_make_sed_table'
  goto,the_end
ENDIF

Nparams=n_elements(parameters_description)

;== INITIALISE DUSTEM
use_polarization=0
dustem_init,model=model,pol=use_polarisation
!quiet=1
;!dustem_verbose=1
;!dustem_show_plot=1
!dustem_which='RELEASE'

;=== define default filters
use_filters=dustem_get_all_filter_names(all_instrus=all_instrus)
IF keyword_set(filters) THEN use_filters=filters
ind=where(use_filters NE 'SPECTRUM',count)
IF count NE 0 THEN BEGIN
	use_filters=use_filters[ind]
ENDIF ELSE BEGIN
	message,'No usable filter found',/info
	stop
ENDELSE

wavs=dustem_filter2wav(use_filters)
order=sort(wavs)
use_filters=use_filters[order]
Nfilters=n_elements(use_filters)

;==== compute all combinations of model parameter values
parameter_values=dustem_param_range2param_values(iv_min,iv_max,iv_Nvalues,Nc=Nc,log=plog)

;=== initialize filters for the SED
Nfilt=n_elements(use_filters)
;stop
sed=dustem_initialize_sed(Nfilt)

sed.filter=use_filters
sed.wave=dustem_filter2wav(use_filters)
sed.instru=dustem_filter2instru(use_filters)

;=== initialize IQU and associated errors to avoid problems when checking SED in dustem_set_data.pro
sed[*].StokesI=1.e-10
sed.StokesQ=sed.StokesI/100.
sed.StokesU=sed.StokesI/100.
sed.SigmaII=sed.StokesI/100.
sed.SigmaQQ=sed.StokesI/100.
sed.SigmaUU=sed.StokesI/100.
sed.SigmaIQ=sed.StokesI/100.
sed.SigmaIU=sed.StokesI/100.
sed.SigmaQU=sed.StokesI/100.

;=== initialize wavelengths for the EXTINCTION
use_next=100 ; number of extinction measurements
use_wave=10^(alog10(0.01)+alog10(50)*findgen(use_next)/float(use_next)) ; wavelengths at which extinction is defined
if keyword_set(next) then use_next=next
if keyword_set(wave) then $
   use_wave=range_gen(use_next,[wave[0],wave[1]],/log)

ext=dustem_initialize_ext(use_next)
ext.instru='EXTINCTION'
ext.filter='SPECTRUM'
ext.wave=use_wave

;=== initialize IQU and associated errors to avoid problems when checking EXT in dustem_set_data.pro
ext[*].EXT_I=1.e-10
ext.EXT_Q=ext.EXT_I/100.
ext.EXT_U=ext.EXT_I/100.
ext.SIGEXTII=ext.EXT_I/1000.
ext.SIGEXTQQ=ext.EXT_I/1000.
ext.SIGEXTUU=ext.EXT_I/1000.
ext.SIGEXTIQ=ext.EXT_I/1000.
ext.SIGEXTIU=ext.EXT_I/1000.
ext.SIGEXTQU=ext.EXT_I/1000.

dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext

;=== initialize at least one parameter as 'free'
;=== this is only to keep dustemwrap happy in its management of variables
;=== no fitting is going to be done
;pd = ['(*!dustem_params).gas.G0'] ; G0
;pval = [-1.e-20]
pd=parameters_description

;=== declare some other parameters as fixed
;=== this shows how to add any non-dust plugins to the observed spectra/SED
;=== for the dust-model parameters, it is not necessary, but allows the
;=== user to (i) change default values and (2) see the adopted values in the graphical output windows
;fpd=parameters_description
seds=ptrarr(Nc)
em_spectra=ptrarr(Nc)
ext_spectra=ptrarr(Nc)
yrange=[1.e-5,1.e5]
colors=long(range_gen(Nc,[0,255]))
IF keyword_set(show_seds) THEN loadct,13
FOR i=0L,Nc-1 DO BEGIN
	;dustem_init,model=model,pol=use_polarisation
	;dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext
  st=0 ;otherwise model is not recomputed !!

	;fpval = *parameter_values[i]
  pval = *parameter_values[i]
	message,'= computing model '+strtrim(i,2)+'/'+strtrim(Nc,2)+' ('+strtrim(1.*i/Nc*100.)+'%)',/continue
	message,'free parameter values=',/continue
	print,'free_parameter_values=',pval
	message,'fixed parameter values=',/continue
	print,'fixed_parameter_values=',fiv
  ;stop
	;=== initialise all this information into dustemwrap's brain
	dustem_init_params,model,pd,pval,fpd=fpd,fiv=fiv,pol=use_polarisation,/force_reset

  ;IF keyword_set(print_params) THEN BEGIN
 ; 	dustem_print_params,0
 ; 	;stop
 ; ENDIF

	;=== compute the predicted SED and extinction curve for the dust-model
	;=== and any plugins
	dustem_Ised=dustem_compute_sed(pval,st=st)   ;is in MJy/sr/1e20 H/cm2
	
	dustem_Qsed=dustem_Ised*0.
	dustem_Used=dustem_Ised*0.

	;== Now we have a predicted SED (I), we fill up a structure that we
	sed.stokesI=dustem_Ised
	sed.stokesQ=dustem_Qsed
	sed.stokesU=dustem_Used
	;=== set uncertainties for I,Q,U to 0.
	sed.sigmaII=0.
	sed.sigmaQQ=0.
	sed.sigmaUU=0.
	;=== set covariances to 0.
	sed.sigmaIQ=0.
	sed.sigmaIU=0.
	sed.sigmaQU=0.

	;==== fill in dependent columns of the SED.
	sed=dustem_fill_sed_dependent_columns(sed)
	seds[i]=ptr_new(sed)
	em_spectra[i]=ptr_new(st.sed)
	ext_spectra[i]=ptr_new(st.ext)
  IF keyword_set(show_seds) THEN BEGIN
  	  xrange=[0.1,1.e3]
	    IF i EQ 0 THEN BEGIN
		    tit='Grid SEDs model '+model
		    xtit='Wavelength [mic]'
		    ytit='SED [MJy/sr for NH=1.e20 H/cm^2]'
		    wavs=(*seds[i]).wave
		    Is=(*seds[i]).stokesI
	    	cgplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,/ylog,yr=yrange,xtit=xtit,ytit=ytit,tit=tit,/xlog,xrange=xrange,/xsty
	    ENDIF ELSE BEGIN
	    	cgoplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,color=colors[i]
	    ENDELSE
	    ;stop
	ENDIF
	;stop
ENDFOR
;stop

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


;toto=mrdfits(fits_file,1,header)
;hprint,header

;stop
the_end:

END