PRO xvclean , tvec, dvec, gplot, flat, tclean, dclean



;; Announce programme and options

IF (gplot ne 2) THEN BEGIN

  print, 'XVCLEAN, tvec, dvec, gplot, flat, tclean, dclean'
  print
  print, 'clean variable star data by removing outliers and background trends'
  print, 'Input:'
  print, 'tvec, dvec  : arrays of time and data values '
  print, 'gplot = 0 (no graphics), 1 (X or PS), 2 (X window), >2 (more) '
  print, 'flat  = 0 (no background), 1 (sine wave), 2 (polynomial) '
  print
  print, 'Output:'
  print, 'tclean, dclean  : cleaned data '
  print 

  IF (gplot eq 0) THEN $
    print, 'gplot =', gplot, '; no graphics'
  IF (gplot eq 1) THEN $
    print, 'gplot =', gplot, '; graphics on / X or PS'
    IF (strupcase(!d.name) EQ 'X') THEN window,3
  IF (gplot eq 2) THEN $
    print, 'gplot =', gplot, '; graphics on / X window'
    IF (strupcase(!d.name) EQ 'X') THEN window,3
  IF (gplot gt 2) THEN $
    print, 'gplot =', gplot, '; diagnostic graphics'
  IF (flat eq 1) THEN $
    print, 'flat  =', flat, '; sine wave background'
  IF (flat eq 2) THEN $
    print, 'flat  =', flat, '; polynomial background'
  print

ENDIF


;;  clean the data

; Initialize
tclean = tvec
dclean = dvec

; Compute sigma, omitting far out data
  nspec = N_ELEMENTS(dvec)
  svec = STDDEV(dvec)
  mvec = TOTAL(dvec)/nspec
  FOR it=0,nspec-1 DO BEGIN
    IF (ABS(dvec(it)-mvec) gt 2*svec) THEN BEGIN
      dclean(it) = mvec
    ENDIF  
  ENDFOR
  svec = stddev(dclean)

; Select data within +/- 3 sigma
ic     = 0
FOR it=0,nspec-1 DO BEGIN
  IF (abs(dvec(it)-mvec) lt 3*svec) THEN BEGIN
     tclean(ic) = tvec(it)
     dclean(ic) = dvec(it)
     ic = ic + 1
  ENDIF  
ENDFOR

IF (ic eq 0) THEN ic = nspec
 tclean = tclean(0:ic-1)
 dclean = dclean(0:ic-1)
 mclean = TOTAL(dclean)/ic
 dclean = dclean - mclean


; plot the raw data 
IF (gplot gt 2) THEN BEGIN
        print, mvec, mclean, svec
        plot,tvec,dvec,        $
             psym=1,                                      $ 
             /xstyle, xtitle='Time',                      $
             ytitle='cleaned data',                       $
             yrange=[mclean-2*svec,mclean+2*svec],        $
             title=ptitle,                                $
             charsize=1.4
ENDIF



;; Subtract low-frequency background


; initialize background fit
bg_fit = fltarr( ic )
bg_fit = 0.


; Polynomial background 
IF ( flat eq 2 ) THEN BEGIN
  bg_coeffs = poly_fit (tclean, dclean, 2, bg_fit)
  dclean = dclean - bg_fit
  print, "Background coefficients"
  print, "bg = a0  +  a1 * x + a2 * x^2"
  print, bg_coeffs, format='(3f10.2)'
ENDIF


; Sine wave background
IF ( flat eq 1 ) THEN BEGIN
  bg_coeffs = [0.,5.,0.55,1.]
  wvec = fltarr( ic )
  wvec = [replicate(1.0, ic)]
  bg_fit = curvefit (tclean, dclean, wvec, $
                      bg_coeffs, function_name='sinefit', $
                      bg_sigma, noderivative=1  )
  dclean = dclean - bg_fit
  print, "Background coefficients"
  print, "bg = a0  +  a1 * sin ( 2.pi.(x-a2) / P )
  print, bg_coeffs
  print, "+/-", bg_sigma
ENDIF


;; Plot the data and the background fit

IF (gplot gt 0) THEN BEGIN
        plot,tclean,dclean+bg_fit,        $
             psym=+1,                                     $ 
             /xstyle, xtitle='Time',                      $
             ytitle='Velocity (km/s)',                    $
             $ ;  yrange=[-50,+50],                            $
             title=ptitle,                                $
             charsize=1.4
        IF (flat gt 0 ) THEN oplot,tclean,bg_fit,color=45
ENDIF


END