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:

;    DustEMWrap, Mid-Level, Distributed

; 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',
;              'FLAMBDA=cste', 'HFI', 'LABOCA'

;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:

;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed

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

;    Written JPB Apr-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_cc'
  sed=0.
  goto,the_end
ENDIF

;eps=1.e-30
eps=0.

;=== constants needed for CMB flux convention

Tcmb=2.725D0

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

  IF finite(cc[i]) THEN sed[i]=spec0*cc[i] ELSE sed[i]=spec0

ENDFOR

the_end:
RETURN,sed

END