;+
;FUNCTION: THM_LSP_FIND_SPIKES
;
;           NOT FOR GENERAL USE. CALLED BY THM_EFI_REMOVE_SPIKES
;           ONLY FOR ISOLATED PARTICLE OR WAVE BURSTS
;           NOT FOR ENTIRE ORBIT.
;
;PURPOSE:
;    Remove the non-physical spiky signals in the efw data.
;
;INPUT:
;    t             -NEEDED. Time array
;    xx            -NEEDED. EFP data, DO NOT GIVE EFW DATA. Won't work.
;    per           -NEEDED. Spin period.
;
;KEYWORD:
;    nwin          -OPTIONAL. Number of points in spike search window. DFLT = 16
;    spikesig      -OPTIONAL. Sigma of spikes. DFLT = 5
;    sigmin        -OPTIONAL. Minimun of sigma. DFLT = 0 mV/m.
;
;OUTPUT:
;    tpks          RETURN VALUE. Time of peaks (mod spin per) + t0
;    Amp           Estimated amplitude of spikes. Helps remove_spikes.
;
;HISTORY:
;   2009-05-30: REE. Broke out to run with wave burst.
;-
function thm_lsp_find_spikes, t, xx, per, nwin=nwin, $
  spikesig=spikesig, sigmin=sigmin, talk=talk, Amp=Amp


IF not keyword_set(width) then width = 1.d-3 ; 1 ms

; SET UP TIME AND DATA
tt = t-t[0]
dt = tt[1:*]-tt[*]
mdt = median(dt)
x = xx

; CHECK KEYWORDS
if not keyword_set(nwin) then nwin = 16L
if not keyword_set(spikesig) then spikesig = 5.0d
if n_elements(sigmin) EQ 0 then sigmin = 0.01d

nwin = long(nwin)

; CALCULATE AVERAGE SIGNAL OVER PERIOD

phs      = tt/per
npps     = long(round(2.0d*per/mdt))    ; SET NUMBER OF POINTS PER SPIN AT ~2X DT
nspins   = ceil(max(phs))               ; NUMBER OF SPINS
ntot     = long(nspins) * long(npps)    ; NUMBER OF PHASE POINTS
phx      = dindgen(ntot)/npps
maxphs   = max(phs)
xp       = interpol([x,0,0], [phs, maxphs+mdt/per, maxphs+2*mdt/per], phx)
xr       = reform(xp, npps, nspins)
xe       = total(xr,2)/ nspins


if keyword_set(talk) then $
  plot, phx[0:npps-1], xe, xtit = 'Spin Phase (/2pi)', $
  title = 'Spin-Averaged Value', ytit = 'Amplitude'

; FIND SPIKES
;xel  = [xe(npps-nwin/2:npps-1), xe, xe(0:nwin/2-1)]
;xel2 = abs(xel-thm_lsp_median_smooth(xel,nwin))
;sig  = abs(xel2)/(stddev(xel2(nwin/2: npps+nwin/2-1))+sigmin)
;if keyword_set(nms) then sig = thm_lsp_median_smooth(sig,nms)
;if keyword_set(nms) then sig = smooth(sig,nms)
;sdum = sig                                 ; ISOLATE MAXIMUMS IN SIG
;FOR i= nwin/2-1, npps+nwin/2-1 DO BEGIN
;  smax = max(sdum(i-nwin/2+1:i+nwin/2))
;  if sdum(i) EQ smax then sig(i)=smax else sig(i) = 0
;ENDFOR
;sig  = sig(nwin/2: npps+nwin/2-1)
;ppk  = where(sig GE spikesig, nppk)


; NEW VERSION
xe3   = [xe, xe, xe]
Amp   = xe3-thm_lsp_median_smooth(xe3,nwin)
xe3s  = abs(Amp)
xe3sq = sqrt(smooth(xe3s*xe3s, npps/4+1))
sig   = xe3s/(xe3sq+sigmin)
sdum  = sig                                 
FOR i= nwin/2-1, 3*npps-nwin/2-1 DO BEGIN
  smax = max(sdum[i-nwin/2+1:i+nwin/2])
  if sdum[i] EQ smax then sig[i]=smax else sig[i] = 0
ENDFOR
sig   = sig[npps: 2*npps-1]
Amp   = Amp[npps: 2*npps-1]
ppk   = where(sig GE spikesig, nppk)


; CHECK PPK FOR ERRORS
IF nppk EQ 0 then BEGIN
  print, 'NO SPIKES FOUND'
  IF keyword_set(talk) then BEGIN
    xyouts, 0.25, 0.75, 'NO SPIKES FOUND', charsize=2, /normal
    wait, 1.0
  ENDIF 
  return, -1
ENDIF
IF nppk GT 4 THEN BEGIN
  print, 'TOO MANY SPIKES - CANNOT REMOVE'
  return, -1
ENDIF 

; PEAKS
tpks  = ppk*per/double(npps) + t[0]

IF keyword_set(talk) then BEGIN
  oplot, phx[ppk], xe[ppk], psym=6, col = 250
  xyouts, 0.25, 0.75, 'SPIKES', charsize=2, /normal
  wait, 0.25
ENDIF 

; ESTIMATE AMPLITUDES
Amp = Amp[ppk]

return, tpks

end