FUNCTION dustem_cc,wavein,spec,filter_names,cc=cc,fluxconv=fluxconv,no_sort=no_sort,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, Color corrections
;
; 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', 'FLAMBDA=cste', 'FLAMBDA=cste'
;              , 'IRAC', 'MIPS', 'CMB', 'HFI', 'LABOCA', 'NIKA2'
;    no_sort   = if set, input arrauys are assumed to already be orderd by increasing wavelengths
;    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
fact3=K*T_0

;sort input, unles /no_sort
IF not keyword_set(no_sort) THEN BEGIN
  order =sort(wavein)
  spec2 = [eps,spec[order]]
  wavein2 = [eps,wavein[order]]
ENDIF ELSE BEGIN
  wavein2=wavein
  spec2=spec
ENDELSE

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))

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)
;  stop
  ii=(where(st.filter_names EQ filter_names[i]))[0]
  spec0=interpol(spec2,wavein2,st.central_wavelengths[ii])
  IF spec0 EQ 0 THEN BEGIN
    sed[i]=spec0
    cc[i]=1
    goto,next_filter
  ENDIF
  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
      'FNU=CSTE': BEGIN
         ;beta=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]))^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+beta))),st.use_wmin(ii),st.use_wmax(ii),/double)               
         ;cc[i]=num/(spec0*den*(st.central_wavelengths[ii])^beta)
         cc[i]=num/(spec0*st.filters_integral[ii])
      END
      'NUINU=CSTE': BEGIN
         ;beta=-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)))^(2+beta))),st.use_wmin(ii),st.use_wmax(ii),/double)               
         ;cc[i]=num/(spec0*den*(st.central_wavelengths[ii])^beta)
         cc[i]=num/(spec0*st.filters_integral[ii])
      END
      'FLAMBDA=CSTE': BEGIN
         ;beta=-2.
         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+beta))),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*(st.central_wavelengths[ii])^2)/(spec0*den)
         ;cc[i]=num/(spec0*den*(st.central_wavelengths[ii])^beta)
         cc[i]=num/(spec0*st.filters_integral[ii])
      END
      'MIPS': BEGIN ;Not used ??
        ;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.)
        cc[i]=num/(spec0*st.filters_integral[ii])
      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))
        cc[i]=num/(spec0*st.filters_integral[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
        cc[i]=num/(spec0*st.filters_integral[ii])
     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.))
        cc[i]=num/(spec0*st.filters_integral[ii])
      END
      'NIKA2': BEGIN
         ;beta=1.6 
         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+beta))),st.use_wmin(ii),st.use_wmax(ii),/double)               
         ;cc[i]=num/(spec0*den*(st.central_wavelengths[ii])^beta)
         cc[i]=num/(spec0*st.filters_integral[ii])
      END
      ELSE:BEGIN
        message,strupcase(fluxconvs[i])+' not recognized, assuming NUINU=CSTE',/info
        ;cc[i]=1. ; changed default behaviour in version 4.3
         ;beta=-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)))^(2+beta))),st.use_wmin(ii),st.use_wmax(ii),/double)               
         ;cc[i]=num/(spec0*den*(st.central_wavelengths[ii])^beta)
         cc[i]=num/(spec0*st.filters_integral[ii])
      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
  next_filter:
ENDFOR

the_end:
RETURN,sed

END