FUNCTION dustem_compute_gb_sed_fast,p_dim,_extra=extra,waves=waves,spec=spec

;+
; NAME:
;    dustem_compute_gb_sed
; PURPOSE:
;    Computes an SED from a given Grey Body spectrum
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    sed=dustem_compute_sed(p_dim[,st=][,cont=][,_extra=][,/help])
; INPUTS:
;    p_dim      = parameter values
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    sed       = computed SED for filters in !dustem_data
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

;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)
;w0=350.
w0=100.

;GET THE OBSERVATIONS AND ERRORS
obs_sed = (*!dustem_data.sed).values
err_sed = (*!dustem_data.sed).sigma

;COMPUTE THE MODEL SED
;dustem_sed = obs_sed*0.D0
;ind_sed=where((*!dustem_data).filt_names NE 'SPECTRUM',count_sed)

;stop
print,'!dustem_do_cc=',!dustem_do_cc
print,'!dustem_never_do_cc=',!dustem_never_do_cc
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
  ;stop
  NNw=200L
  wwmin=50. & wwmax=10000.
  waves=dindgen(NNw)/(1.*NNw-1)*(alog10(wwmax)-alog10(wwmin))+alog10(wwmin)
  waves=10.D0^waves
  ;waves=dindgen(NNw)/(1.*NNw-1)*(wwmax-wwmin)+wwmin
  spec=pp(0)*(waves/w0)^(-1.*pp(2))*dustem_planck_function(pp(1),waves)
;    dustem_sed(ind_sed)=dustem_cc(waves,spec,((*!dustem_data).filt_names)(ind_sed),cc=cc)
  dustem_sed=dustem_cc(waves,spec,((*!dustem_data.sed).filt_names),cc=cc)
  spec0=interpol(spec,waves,dustem_filter2wav((*!dustem_data.sed).filt_names))
  ;dustem_sed=spec0*cc
  print,dustem_sed,cc,pp,spec0
;stop
;  ENDIF
ENDIF ELSE BEGIN ;(!dustem_do_cc EQ 0 OR !dustem_never_do_cc EQ 1)
  message,'NOT DOING color correction calculations',/info
;  waves=1.d0*dustem_filter2wav(((*!dustem_data).filt_names)(ind_sed))
  waves=1.d0*dustem_filter2wav(((*!dustem_data.sed).filt_names))
  spec=pp(0)*(waves/w0)^(-1.*pp(2))*dustem_planck_function(pp(1),waves)
  dustem_sed=spec*(*!dustem_previous_cc)
  print,dustem_sed,*!dustem_previous_cc,pp,spec
;  print,dustem_sed,spec,waves,*!dustem_previous_cc
;  stop
;  dustem_sed(ind_sed)=spec*(*!dustem_previous_cc)
ENDELSE

;!dustem_do_cc=       1
;!dustem_never_do_cc=       0
;% DUSTEM_COMPUTE_GB_SED_FAST: DOING color correction calculations
;       27736.339       11914.539       3501.9079       850.50556
;      0.97084945       1.0627671       1.0643585       1.1851495
;       1.0000000       20.000000       2.0000000
;       28569.145       11210.865       3290.1583       717.63569
;IDL>   print,dustem_sed,*!dustem_previous_cc,pp,spec
;       28987.803       11820.542       3111.8034       799.60098
;      0.97084945       1.0627671       1.0643585       1.1851495
;       1.0000000       20.000000       2.0000000
;       29858.186       11122.420       2923.6422       674.68365

;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).filt_names EQ 'SPECTRUM',count_spec)
;IF count_spec NE 0 THEN BEGIN
;  dustem_sed(ind_spec)=interpol(spec,ww,(((*!dustem_data).wav)(ind_spec)))
;;  dustem_sed(ind_spec)=10^interpol(alog10(spec),alog10(st.dustem.wav),alog10((((*!dustem_data).wav)(ind_spec))))
;ENDIF

;clean pointers
;heap_gc

the_end:

RETURN,dustem_sed

END