FUNCTION dustem_cc,wavein,spec,filter_names,cc=cc,fluxconv=fluxconv,help=help

;+
; NAME:
;    dustem_cc
; PURPOSE:
;    Computes observed sed and color correction for a given spectrum in a given filter
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    sed=dustem_cc(wave,spec,filter_names,cc=cc,fluxconv=fluxconv,help=help)
; INPUTS:
;    wave: array of wavelengths for spec
;    spec: spectrum (must be in brightness units)
;    filter_names: names of filters for which color correction is needed
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    sed: SED as observed in filters provided in filter_names
; OPTIONAL OUTPUT PARAMETERS:
;    cc        = color correction coefficients (spec*cc=sed)
; ACCEPTED KEY-WORDS:
;    fluxconv  = if set, these are taken for flux conventions
;              possible values: 'nuInu=cste', 'IRAC', 'MIPS', 'CMB'
;    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.P. Bernard April 18th 2011
;    replaces individual cc routines for various instruments
;    see evolution details on the dustem cvs maintained at CESR
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

;For MIPS color corrections, see
;http://irsa.ipac.caltech.edu/data/SPITZER/docs/mips/mipsinstrumenthandbook/51/

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_cc'
  sed=0.
  goto,the_end
ENDIF

;eps=1.e-30
eps=0.
;=== needed for CMB flux convention
Tcmb=2.725D0

;=== needed for MIPS flux convention
h=6.62606876*1d-34            ; J s
c=2.99792458*1d8              ; m*s-1
K=1.380658*1d-23              ; J/K 
T_0=1d4                       ; calibration star (K units)
fact=h*c
;fact1=h*c
fact3=K*T_0

; sort input
order =sort(wavein)
spec2 = [eps,spec(order)]
wavein2 = [eps,wavein(order)]

IF not keyword_set(fluxconv) THEN BEGIN
  fluxconvs=dustem_filter2fluxconv(filter_names)
ENDIF ELSE BEGIN
  fluxconvs=fluxconv
ENDELSE
instrus=dustem_filter2instru(filter_names)

inst_tags=tag_names((*!dustem_filters))

;DEBUG VG

;j=where(instrus eq 'HFI',count)
;if count gt 0 then instrus(j) = 'HFI_V3P02'

Nbands=n_elements(filter_names)
sed=dblarr(Nbands)

cc=dblarr(Nbands)
cc(*)=1.d0

FOR i=0L,Nbands-1 DO BEGIN
  iinstru=where(inst_tags EQ instrus(i),countinstr)

  IF countinstr EQ 0 THEN BEGIN
    message,'Instrument '+instrus(i)+' not found',/continue
    stop
  ENDIF

  st=(*!dustem_filters).(iinstru)
  ii=(where(st.filter_names EQ filter_names(i)))(0)
  spec0=interpol(spec2,wavein2,st.central_wavelengths(ii))
  IF !dustem_do_cc NE 0 and not !dustem_never_do_cc THEN BEGIN   ;Color corrections are computed
;    print,'DOING color correction calculations'
    CASE strupcase(fluxconvs(i)) OF
      'HFI0': BEGIN  ;This is purely for test (stupid since reads filters each time)
        filter_dir=!dustem_wrap_soft_dir+'/Data/FILTERS/'
        dir_hfi=filter_dir+'HFI/OFFICIAL'+'/'
;       CAUTION: This below should be consistent with what is in dustem_read_filters.pro
;        hfi_filter_version='v101'
        hfi_filter_version='v201'
        bp = hfi_read_bandpass(hfi_filter_version, /rimo,path_rimo=dir_hfi)
        spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
        nu=*(st.use_frequencies(ii))/1.d9
        order=sort(nu)
        spec_int=spec_int(order)
        nu=nu(order)
        channel=strtrim(round(dustem_filter2freq(filter_names(i))),2)
;stop
        ccc=hfi_color_correction_dustem(bp.name,bp.freq/1.d9,bp.trans,channel, bolo_cc,nu, spec_int)
;        ccc=hfi_color_correction(bp,channel, bolo_cc,nu, spec_int)
        cc(i)=ccc.cc
      END
      'HFI': BEGIN
        spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
        nu=*(st.use_frequencies(ii))/1.d9
        trans=*(st.use_transmissions(ii))
        order=sort(nu)
        spec_int=spec_int(order)
        nu=nu(order)
        trans=trans(order)
        channel=strtrim(round(dustem_filter2freq(filter_names(i))),2)
        ccc=hfi_color_correction_dustem(channel,nu,trans,channel, bolo_cc,nu, spec_int)
        cc(i)=ccc.cc
      END
      'NUINU=CSTE': BEGIN

         spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
         num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii)))^2.,st.use_wmin(ii),st.use_wmax(ii),/double)
         den=integral(*(st.use_wavelengths(ii)),(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)               
         cc(i)=num/den/spec0*(st.central_wavelengths(ii))
      END
      'FLAMBDA=CSTE': BEGIN
         spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
         num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii)))^2.,st.use_wmin(ii),st.use_wmax(ii),/double)
         den=integral(*(st.use_wavelengths(ii)),(*(st.use_transmissions(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)               
         cc(i)=num/den/spec0*(st.central_wavelengths(ii))^2

      END
      'MIPS': BEGIN
        expon=exp(fact/(*(st.use_wavelengths(ii))*1e-6*fact3))-1.
        spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
        num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii)))^2.,st.use_wmin(ii),st.use_wmax(ii),/double)
        den=integral(*(st.use_wavelengths(ii)),*(st.use_transmissions(ii))/(*(st.use_wavelengths(ii)))^5/expon,st.use_wmin(ii),st.use_wmax(ii),/double)
        ffact=st.central_wavelengths(ii)*1e-6*K*T_0
        cc(i)=num/den/spec0/(st.central_wavelengths(ii))^3/(exp(fact/ffact)-1.)
      END
      'IRAC': BEGIN
        spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
        num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
        den=integral(*(st.use_wavelengths(ii)),*(st.use_transmissions(ii)),st.use_wmin(ii),st.use_wmax(ii),/double)
        cc(i)=num/den/spec0*(st.central_wavelengths(ii))
      END
      'CMB': BEGIN
        spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
        num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii)))^2.,st.use_wmin(ii),st.use_wmax(ii),/double)
        den=integral(*(st.use_wavelengths(ii)),(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))^2)*dustem_planck_function(Tcmb,*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
        cc(i)=num/den*dustem_planck_function(Tcmb,st.central_wavelengths(ii))/spec0
     END
      'LABOCA': BEGIN
;        spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
;        num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
;        den=integral(*(st.use_wavelengths(ii)),*(st.use_transmissions(ii))/(*(st.use_wavelengths(ii))^2.5),st.use_wmin(ii),st.use_wmax(ii),/double)
;        cc(i)=num/den/spec0/(st.central_wavelengths(ii))^2.5

        alpha_new=0.

        spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
        num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
        den=integral(*(st.use_wavelengths(ii)),*(st.use_transmissions(ii))/(*(st.use_wavelengths(ii))^(alpha_new+2.)),st.use_wmin(ii),st.use_wmax(ii),/double)
        cc(i)=num/den/spec0*(st.central_wavelengths(ii))^(-(alpha_new+1.))

      END

      ELSE:BEGIN
        message,strupcase(fluxconvs(i))+' not recognized',/info
        cc[i]=1.
        ;stop
      END

    ENDCASE
    !dustem_previous_cc=ptr_new(cc)
  ENDIF ELSE BEGIN ;When no new cc calculation is required
    IF !dustem_never_do_cc EQ 0 THEN BEGIN
      cc=(*!dustem_previous_cc)
    ENDIF
  ENDELSE
  sed(i)=spec0*cc(i)
ENDFOR

the_end:
RETURN,sed

END