; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

      PRO thm_SCM_casinus_vec, xi,fe,fs,ns,amp,pha,n,nbi, xo
 
;     ------------------------------------------------------------------
; *   Object : Compute a pure sine from a given signal
; *   Class  : Data processing for GEOS/UBF
; *   Author : P. Robert, CRPE, 1977-1984
;            : Conversion Fortran to IDL, P. Robert, CETP, May 2007
;            : Vectorized by K. Bromund, SPSystems/GSFC, May 2007
;     ------------------------------------------------------------------
 
 
;     The input signal contains a high sine signal, with a known
;     frequency, superimposed to a useful signal.
 
;     Input  :
;          xi:  array containing the signal
;          fe:  sampling frequency 
;          fs:  frequency of the sine signal to be removed
;          ns:  number of spins to use to comput sine wave fit
;     Output :
;         amp:  amplitude of the compute sine signal
;         pha:  phase (d) of the compute sine signal
;           n:  number of points in output
;         nbi:  number of point taken for sine computation
;          xo:  sine-subtracted signal
 
;     ------------------------------------------------------------------
 
; *** computation of the number of points coresponding to an integer
;     number of the sine period to remove
 
      n= N_ELEMENTS(xi)
      nbsp=LONG(FLOAT(n)*fs/fe)  ; number of spins in input data
      nbi =LONG((FLOAT(ns)*fe/fs)+0.5) ;; nuber of points required to 
                                       ;; fit requested number of spins
      
      if not (nbi mod 2) then nbi+=1  ;; algorithm seems to work better w/ odd
; *** test if one have enough points

      IF (nbsp EQ 0 OR nbi GT n) THEN BEGIN
               PRINT, '*** desinus: no enough points to remove the sine'
               PRINT, '    at least one sine period is required'
               PRINT, '    sampling frequency=',fe
               PRINT, '    sine frequency    =',fs
               PRINT, '    Number of spins requested for fit: ', ns
               PRINT, '    number of points required:',nbi, $
                      '    available=',n
               amp=0.
               pha=0.
               RETURN
               ENDIF
 
      ; compute the sine waves for fitting
      phase_s = dindgen(nbi)*2*!dpi*fs/fe  ;; phase of each sample in kernel
      ss=sin(phase_s) * hanning(nbi) * 2.0
      cc=cos(phase_s) * hanning(nbi) * 2.0

      center = nbi/2  ;integer division intentional

      zs = convol(xi, ss)
      zc = convol(xi, cc)

      amp=2.*sqrt(zs*zs+zc*zc)/FLOAT(nbi)
      pha=ATAN(zc,zs)  ;; this is correct, because we subtract sin rather than
                       ;; cosine  wave.


      
; *   Fill out phase and amp results for first half of first window and
;     last half of last window

      ;; Note: this method gets  a phase discontinuity when nbi is odd, but
      ;; gets amplitude discontinuity at extremes when nbi is even, even
      ;; for a constant sine input.  That is why I now force nbi to be odd.

      ;; Consider truncating data by 1/2 of sliding window at beginning and end.

      ;; some data:
      ;;  0.03% residual for 10 spins
      ;;  2% residual for 1 spin (nbi even) (extremities worse)
      ;;  1% residual for 1 spin, forced to odd nbi=49
      ;;  0.03% residual for 2 spins (nbi odd)
      ;;  0.3% residual fo 3 spins (nbi odd)
      ;;  0.03% for 4 spins (nbi even) phase discont only at begin,
      ;;                         but at end, .6/5 = 0.1
      ;;  0.54% for 4 spins (nbi forced to odd: 195), 
      ;;     phase discontinuities only
      ;; 
      ;;  2 spins seems the best choice.  Also works well when amplitude varies

      pha[0:center-1] = phase_s[0:center-1] + pha[center]
      amp[0:center-1] = amp[center]

      pha[n-center:n-1] = phase_s[nbi-center:nbi-1] + pha[n-1-center]     
      amp[n-center:n-1] = amp[n-center-1]

      pha[center:n-center-1] += phase_s[center] 

; *  subtract sine wave while pha is still in pha radians
      xo = xi - amp * sin(pha)

; *  convert phase in degree for output
 
      pha *= 180./!dpi

      END

; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX