PRO dustem_make_sed_table,model, $
											    parameters_description, $
											    iv_min, $
											    iv_max, $
											    iv_Nvalues, $
						  						fpd=fpd, $
						  						fiv=fiv, $
						  						filename=filename, $
						  						filters=filters, $
						  						log=log, $
						  						show_seds=show_seds, $
						  						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 = IRAS filters)
;    fpd                    : fixed parameter description
;    fiv                    : fixed parameter values
;    filename               : output fits file name for the table (default ='/tmp/dustem_seds.fits')
;    log                    : 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:
;    None
; 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',log=[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_instru2filters('NIRCAM'),dustem_instru2filters('MIRI'),dustem_instru2filters('IRAS'),dustem_instru2filters('PACS'),dustem_instru2filters('SPIRE'),dustem_instru2filters('HFI')]
IF keyword_set(filters) THEN use_filters=filters
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=log)

;=== 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)
yrange=[1.e-5,1.e5]
colors=long(range_gen(Nc,[0,255]))
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
  dummy=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
	;print,'parameter_values=',pval
    
	;=== initialise all this information into dustemwrap's brain
	;stop
	;dustem_init_params,model,pd,pval,fpd=fpd,fiv=fpval,pol=use_polarisation
	dustem_init_params,model,pd,pval,fpd=fpd,fiv=fiv,pol=use_polarisation

	;=== compute the predicted SED and extinction curve for the dust-model
	;=== and any plugins 
	dustem_Ised=dustem_compute_sed(pval,st=dummy)
	dustem_Qsed=dustem_Ised*0.
	dustem_Used=dustem_Ised*0.
	; uncomment following lines if running a polarisation model and you
	;want predictions of QU
	;toto=dustem_compute_stokes(pval,st=dummy) ;this procedure also allows for the extraction of the spectra
	;dustem_qsed = toto[0]
	;dustem_used = toto[1]

    ;This is for extinction. Not used at the moment.
	;dustem_iext=dustem_compute_ext(pval,st=dummy)
	;dustem_Qext=dustem_Iext*0.
	;dustem_Uext=dustem_Iext*0.
	; uncomment following lines if running a polarisation model and you
	;want predictions of QU
	;toto=dustem_compute_stokext(pval,st=dummy) ;this procedure also allows for the extraction of the spectra
	;dustem_qext = toto[0]
	;dustem_uext = toto[1]

	;== Now we have a predicted SED (I), we fill up a structure that we
	;== will write to the sed_outfile
	;== if you 
	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)
  IF keyword_set(show_seds) THEN BEGIN
	    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]'
	    	cgplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,/ylog,yr=yrange,xtit=xtit,ytit=ytit,tit=tit,/xlog
	    ENDIF ELSE BEGIN
	    	cgoplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,color=colors[i]
	    ENDELSE
	ENDIF
	;stop
ENDFOR
;stop

str='one_sed={'
FOR i=0L,Nparams-1 DO BEGIN
  str=str+'P'+strtrim(i+1,2)+':0.d0,'
ENDFOR
FOR i=0L,Nfilters-1 DO BEGIN
	;istr=strtrim(i+1,2)
	istr=strtrim(use_filters[i],2)
	str=str+'I'+istr+':0.d0'
	IF use_polarization EQ 1 THEN str=str+',Q'+istr+':0.d0,U'+istr+':0.d0'
	;IF i NE Nfilters-1 THEN str=str+',' ELSE str=str+'}'
	str=str+','
ENDFOR
str=str+'total:0.d0}'
toto=execute(str)

;stop
Noffset=1
IF use_polarization EQ 1 THEN Noffset=3
seds_array=replicate(one_sed,Nc)
FOR i=0L,Nc-1 DO BEGIN
	FOR j=0L,Nparams-1 DO BEGIN
		seds_array[i].(j)=(*parameter_values[i])[j]
	ENDFOR
	sum=0.d0
	FOR j=0L,Nfilters-1 DO BEGIN
		pos=Nparams+j*Noffset
		seds_array[i].(pos)=((*seds[i]).stokesI)[j]
		sum=sum+((*seds[i]).stokesI)[j] ;======== compute the sum of SEDs
	ENDFOR
	seds_array[i].total=sum
ENDFOR

;======== save SED tables
;use_sed_outfile='/tmp/dustem_seds.xcat'
;write_xcat,seds_array,use_sed_outfile
;message,'Wrote '+use_sed_outfile,/continue

fits_file='/tmp/dustem_seds.fits'
IF keyword_set(filename) THEN fits_file=filename

mwrfits,seds_array,fits_file,header,/create

st=mrdfits(fits_file,0,h)
st1=mrdfits(fits_file,1,h1)

sxaddpar,h1,'MODEL',model,'dustemwrap model used',after='TFORM'+strtrim(Nparams+Nfilters,2)
sxaddpar,h1,'POLAR',use_polarization,'dustemwrap polarization ?',after='MODEL'
FOR j=0L,Nparams-1 DO BEGIN
	value=sxpar(h1,'TTYPE'+strtrim(j+1,2))
	sxaddpar,h1,'TTYPE'+strtrim(j+1,2),value,'dustemwrap '+parameters_description[j]
	sxaddpar,h1,'PNAME'+strtrim(j+1,2),parameters_description[j],'dustemwrap param name'
	;iv_min,iv_max,iv_Nvalues,Nc=Nc,log=log
	sxaddpar,h1,'PMIN'+strtrim(j+1,2),iv_min[j],'minimum param value'
	sxaddpar,h1,'PMAX'+strtrim(j+1,2),iv_max[j],'maximum param value'
	sxaddpar,h1,'PNVAL'+strtrim(j+1,2),iv_Nvalues[j],'Number of param values'
	sxaddpar,h1,'PLOG'+strtrim(j+1,2),log[j],'1 if log scale sampling'
	last_pname='PLOG'+strtrim(j+1,2)
	IF j EQ 0 THEN first_pname='PNAME'+strtrim(j+1,2)
ENDFOR
j0=j
sxaddpar,h1,'COMMENT','**** Dustemwrap PARAMETERS ***',before=first_pname
sxaddpar,h1,'COMMENT','**** Dustemwrap FILTERS ***',after=last_pname
FOR j=0L,Nfilters-1 DO BEGIN
	jj=j+j0
	value=sxpar(h1,'TTYPE'+strtrim(jj+1,2))
	sxaddpar,h1,'TTYPE'+strtrim(jj+1,2),value,'Intensity in '+use_filters[j]
	sxaddpar,h1,'FNAME'+strtrim(j+1,2),use_filters[j],'dustemwrap filter name'
ENDFOR

;stop

mwrfits,seds_array,fits_file,h1,/create
message,'Wrote '+fits_file,/continue

toto=mrdfits(fits_file,1,header)

;hprint,header

;stop
the_end:

END