FUNCTION dustem_BAND_CC, w0, xs, ys, xt, yt, y0=y0, cc=cc, rr=rr, nint=nint
; returns flux in given instrument band (erg/cm2/sr)
; W0    (I): centroid of filter (any unit)
; XS    (I): input x-grid for SED (same unit as W0)
; YS    (I): input SED in erg/cm2/sr
; XT  (I/O): x-grid of filter (same unit as W0)
; YT  (I/O): filter transmission (arb. unit)
; RR    (O): undersampling ratio: dx(SED)/dx(filt), usually 1. 
;            If RR<1 filter oversampled in flux integration
; Y0    (O): SED at band center
; CC    (O): color correction
; NINT  (O): nr of points in band integration

 vlow = 1d-5

; remove bad and multiple values in filter
  xtt=xt & ytt=yt
  iq = UNIQ(xtt,sort(xtt))
  xtt = xtt(iq) & ytt = ytt(iq)
  in = where( (xtt GT 0d) AND (ytt GT 0d), cn )
  if cn GT 0 then begin
     xtt = xtt(in)
     ytt = ytt(in)
  endif else begin 
; if no good point in filter, band flux and the rest is 0 
     sed=0d & rr=0d & y0=0d & cc=0d & nint=n_elements(xt)
     return,sed
  endelse
  nfilt = n_elements(xtt)

; sort arrays  
  ix = SORT(xtt)
  xtt = xtt(ix) & ytt = ytt(ix)
  ytt = ytt / MAX(ytt)
  ix = SORT(xs)
  xss = xs(ix) & yss = ys(ix)
  xt=xtt & yt=ytt

; shrink SED array
  ix = WHERE( xss GE MIN(xtt) AND xss LE MAX(xtt), cx)
  if cx GT 0 then begin
; add points at edges for interpolation
     i1 = 0 
     if ix(0) GT i1 then i1 = ix(0)-1
     i2 = n_elements(xss)-1
     if ix(cx-1) LT i2 then i2 = ix(cx-1)+1
     xss = xss(i1:i2)
     yss = yss(i1:i2)
  endif else begin 
; if no good point in SED, band flux and the rest is 0 
     sed=0d & rr=0d & y0=0d & cc=0d & nint=n_elements(xt)
     return,sed
  endelse
  nsed = n_elements(xss)

; get filter sampling f_res = dx(filt)
  ix = WHERE( ytt/MAX(ytt) GT vlow, cx )
  f_res = ABS( MEDIAN( xtt(ix(1:cx-1))-xtt(ix(0:cx-2)) ) )

; Interpolate to compute integration
  IF NFILT GE NSED THEN BEGIN
; finer filter grid: integrate band flux on filter x-grid
     xi = xtt
     rr = 1d                    ; SED sampling is filter sampling
     yi = INTERPOL(yss,xss,xi)
     filt = ytt
  ENDIF ELSE BEGIN
; finer SED grid: integrate band flux on SED x-grid
     xi = xss & yi = yss 
; get SED sampling: s_res = dx(SED) and sampling ratio
     s_res = ABS( MEDIAN( xi(1:nsed-1)-xi(0:nsed-2) ) )
     rr = s_res / f_res
     filt = INTERPOL(ytt, xtt, xi)
     io = WHERE( xi GT MAX(xtt) OR xi LT MIN(xtt), co)
     if co GT 0 then filt(io) = 0d
     in = where( filt LT 0d, cn) 
     if cn GT 0 then filt(in) = 0d
  ENDELSE
  xt=xi & yt=filt

; Integrate on filter bandpass (log-grid)
; trapezium safe and accurate enough - 
  np = n_elements(xi)
  nint = np
  dx = xi(1:np-1) - xi(0:np-2)
  yint = yi*filt/xi
;  a1 = INT_TABULATED(xi,yint, /double)
  yint = (yint(0:np-2) + yint(1:np-1)) / 2d0
  a1 = TOTAL( yint*dx ) 
  yint = filt/xi
;  a2 = INT_TABULATED(xi,yint, /double)
  yint = (yint(0:np-2) + yint(1:np-1)) / 2d0
  a2 = TOTAL( yint*dx )
  
; get color correction, color corrected band flux and sed
  y0 = INTERPOL( yi, xi, w0)
  cc = a1/a2/y0
  sed = y0*cc 

  RETURN, sed
END