FUNCTION dustem_compute_gb_sed,p_dim $
                               ,_extra=extra $
                               ,waves=waves $
                               ,spec=spec $
                               ,help=help

;+
; NAME:
;    dustem_compute_gb_sed
; PURPOSE:

;    Computes an SED for a modified blackbody 

; CATEGORY:
;    Dustem
; CALLING SEQUENCE:

;    sed=dustem_compute_sed(p_dim[,st=][,_extra=][,/help])

; INPUTS:

;    p_dim      = [a0,a1,a2] modified blackbody parameter values
;                 a0 -- normalisation
;                 a1 -- dust temperature
;                 a2 -- emissivity

; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    sed       = computed SED for filters in !dustem_data.sed
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:

;    help      = if set, print this help

; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:

;    The DustEMWrap IDL code must be installed

; PROCEDURE:
;    None
; EXAMPLES
;    
; MODIFICATION HISTORY:

;    Written JP Bernard  2011
;    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_compute_gb_sed'
  dustem_sed=0.
  goto,the_end
ENDIF

;==== Compute Black body spectrum
pp=double(p_dim)

NNw=20L
wwmin=50. & wwmax=10000.

ww=dindgen(NNw)/(1.*NNw-1)*(alog10(wwmax)-alog10(wwmin))+alog10(wwmin)
ww=10.^ww
waves=ww

;MAKE A COPY OF THE (EMPTY) OBSERVED VALUES

;stop
obs_sed = (*(*!dustem_data).sed).values

;obs_sed = (*!dustem_data.sed).values

;COMPUTE THE MODEL SED
dustem_sed = obs_sed*0.D0

ind_sed=where((*(*!dustem_data).sed).filt_names NE 'SPECTRUM',count_sed)

IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
  message,'DOING color correction calculations',/info
ENDIF ELSE BEGIN
  message,'SKIPPING color correction calculations',/info
ENDELSE

;stop
IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
;  message,'DOING color correction calculations',/info
  IF count_sed NE 0 THEN BEGIN
    spec=pp(0)*(ww)^(-1.*pp(2))*dustem_planck_function(pp(1),ww)
;    filter_names=((*!dustem_data).filt_names)(ind_sed)
;    FOR ii=0L,n_elements(filter_names)-1 DO BEGIN

      ssed=dustem_cc(ww,spec,((*(*!dustem_data).sed).filt_names)(ind_sed),cc=cc)

;print,cc
;    ENDFOR

    dustem_sed[ind_sed]=ssed

  ENDIF
ENDIF ELSE BEGIN

  www=dustem_filter2wav(((*(*!dustem_data).sed).filt_names)[ind_sed])

  dustem_sed[ind_sed]=pp(0)*(www)^(-1.*pp(2))*dustem_planck_function(pp(1),www)*(*!dustem_previous_cc)

  message,'SKIPPING color correction calculations',/info
;  dustem_sed(ind_sed)=interpol(spec,ww,1.D0*(((*!dustem_data).wav)(ind_sed)))*(*!dustem_previous_cc)

spec= dustem_sed[ind_sed]

ENDELSE

;stop
;For spectrum data points, interpolate in log-log
;Linear interpolation leads to wrong values, in particular where few
;wavelengths points exist in the model (long wavelengths).

ind_spec=where((*(*!dustem_data).sed).filt_names EQ 'SPECTRUM',count_spec)

IF count_spec NE 0 THEN BEGIN

  dustem_sed[ind_spec]=interpol(spec,ww,(((*!dustem_data.sed).wav)[ind_spec]))

;  dustem_sed(ind_spec)=10^interpol(alog10(spec),alog10(st.sed.wav),alog10((((*!dustem_data.sed).wav)(ind_spec))))
ENDIF

;clean pointers
heap_gc

the_end:

RETURN,dustem_sed

END