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

;dustem_make_sed_table,'A',['P1','P2','P3'],[0.,2.,4.],[1.,3.,5.],[3,3,3]
Nparams=n_elements(parameters_description)

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

use_filters=['IRAS1','IRAS2','IRAS3','IRAS4']
IF keyword_set(filters) THEN use_filters=filters
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.]

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

	message,'================= computing model '+strtrim(i,2)+'/'+strtrim(Nc,2)+' ('+strtrim(1.*i/Nc*100.)+'%)',/continue
	print,'parameter_values=',fpval
    
	;=== initialise all this information into dustemwrap's brain
	dustem_init_params,model,pd,pval,fpd=fpd,fiv=fpval,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
	    	cgplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,/ylog,yr=yrange
	    ENDIF ELSE BEGIN
	    	cgoplot,(*seds[i]).wave,(*seds[i]).stokesI,psym=-4,/ylog,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(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+'}'
ENDFOR
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
	FOR j=0L,Nfilters-1 DO BEGIN
		pos=Nparams+j*Noffset
		seds_array[i].(pos)=((*seds[i]).stokesI)[j]
	ENDFOR
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'
	last_pname='PNAME'+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

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

toto=mrdfits(fits_file,1,header)

;hprint,header

;stop

END