;5/17/2000 - emi removal tools
 ;=====================================================================
func peakzone(p, ipeak, vrange, irange)
 ;specialized, assumes p is n by 3 array, returns indices of peak "zone"
 n = dimen(p,0)
 yq = p(ipeak, 1)		;our (supposed) peak
 ty,'peak =', yq
 yq = yq * vrange		;take only values larger than this
 ty,'limit =', yq
 ;a bit map for our hits
 hits = zero(intarr(n,3))
 ;go off in 6 directions
 for j=0,2 do {
  i = ipeak
  if p(i,j) ge yq then {
   hits(i,j) = 1
   ;and check further downstream
   while (1) {
    if i le 0 or i le (ipeak - irange) then break
    i = i - 1
    if p(i,j) ge yq then hits(i,j) = 1 else break
    }
   ;and further upstream
   while (1) {
    if i ge n-1 or i ge (ipeak + irange-1) then break
    i = i + 1
    if p(i,j) ge yq then hits(i,j) = 1 else break
    }
   }
  }
 ;have all our hits now
 in = sieve(hits ne 0)
 if isarray(in) then {
 in_i = in%n
 in_j = in/n
 ty,in_i,in_j
 } else ty,'no hits'
 
 
 return, in
 endfunc
 ;=====================================================================
block test1
 ;assume window 0 for initial testing and data in x
 w = 0
 nx = dimen(x,0)
 ny = dimen(x,1)
 sc2d,x,a,b,c,d
 ;we want to examine power spectrum so don't need a displayable version,
 ;also must keep the transforms for inverse
 xq=(a*a+b*b+c*c+d*d)
 aq=2.*(c*d-a*b)
 dq=xq+aq
 aq=xq-aq
 ;not concerned about the edges
 ;look for peaks along the 12.5 and 42 KHz lines
 i12k = fix(26*ny/1024)
 i42k = fix(86*ny/1024)
 ty,'power lines =', i12k, i42k
 nf = dimen(a,1)*2 -1
 ncenter = nf/2
 p12 = fltarr(nf,3)
 p42 = fltarr(nf,3)
 for j = -1, 1 do {
 zq = aq(i12k+j,*)
 zq = reverse(zq(1:*))
 p12(0,j+1) = [zq, dq(i12k+j,*)]
 zq = aq(i42k+j,*)
 zq = reverse(zq(1:*))
 p42(0,j+1) = [zq, dq(i42k+j,*)]
 }
 ;;xport,21
 ;;plot,p12
 ;;xport,22
 ;;plot,p42
 endblock
 ;=====================================================================
block test2
 ;establish smooth bases for locating peaks, use width of 21 in the slow
 ;(small freq step or CCD column direction) and 3 in the fast direction
 ps = 21.0*smooth(p12(*,+), 21)
 pp = p12(*,1)
 pbase = (ps - pp)*(1./63.)
 dp12 = pp - pbase
 dpr12 = dp12/pbase
 ;reject very small peaks
 peak12_limit = 0.02
 peak12_irange = 15
 nx = dimen(x,0)
 peak12_dcrange = 35*nx/1024	;avoid the low freqs
 ty,'peak12_dcrange =', peak12_dcrange
 peak12_vrange = 0.02
 peak12_ncomp = 5

 dpr12 = dpr12 * (pp gt peak12_limit)
 peak12 = maxloc(dpr12)
 peak12 = peak12(0)
 ty,'max for 12.5 is at ', peak12
 pv = !lastmax
 ;check if too close to DC
 if abs(peak12 - ncenter) lt peak12_dcrange then {
   pv = 0
   ty,'peak too close to DC area'
   }
 ;check if too many competitors
 nc = num_elem(sieve(dpr12 gt .5*dpr12(peak12)))
 ty,'competitors at least 0.5 of top peak =', nc
 if nc gt peak12_ncomp then pv = 0
 if pv lt peak12_limit then {
  in12 = 0	;means no points will be removed
  } else {
  ;poke around to find the size of the peak
  in12 = peakzone(p12, peak12, peak12_vrange, peak12_irange)
 }
 ;convert to coordinates in full array
 if isarray(in12) then in12 = in12 + nf*(i12k-1)
 ;the 42 KHz case
 ;establish smooth bases for locating peaks
 ps = 21.0*smooth(p42(*,+), 21)
 pp = p42(*,1)
 pbase = (ps - pp)*(1./63.)
 dp42 = pp - pbase
 dpr42 = dp42/pbase
 peak42_limit = 0.01
 peak42_irange = 10
 peak42_dcrange = 10	;avoid the low freqs
 peak42_vrange = 0.01
 peak42_ncomp = 5
 dpr42 = dpr42 * (pp gt peak42_limit)
 peak42 = maxloc(dpr42)
 peak42 = peak42(0)
 ty,'max for 42 is at ', peak42
 pv = !lastmax
 ;check if too close to DC
 if abs(peak42 - ncenter) lt peak42_dcrange then {
   pv = 0
   ty,'peak too close to DC area'
   }
 ;check if too many competitors
 nc = num_elem(sieve(dpr42 gt .5*dpr42(peak42)))
 ty,'competitors at least 0.5 of top peak =', nc
 if nc gt peak42_ncomp then pv = 0
 if pv lt peak42_limit then {
  in42 = 0	;means no points will be removed
  } else {
  ;poke around to find the size of the peak
  in42 = peakzone(p42, peak42, peak42_vrange, peak42_irange)
 }
 if isarray(in42) then in42 = in42 + nf*(i42k-1)
 if isarray(in12) then {
   if isarray(in42) then inall = [in12, in42] else inall = in12
   } else inall = in42
 xport, 21
 plot, dpr12, 2
 xport,22
 plot,dpr42, 2
 ;show these as a bit array that matches the top half of the 2-D power spectrum
 nf2 = dimen(a,0)
 p2d = zero(bytarr(nf, nf2))
 ;the inall indices are already oriented for top half of p2d
 p2d(inall) = 1
 p2d=concat(reverse(p2d),p2d(*,1:*))

 $data4 = p2d
 view, 4
 
 ;now zero these bad guys
 
 nx = dimen(a,0)	ny = dimen(a,1)
 ;note redef of nx and ny above
 iyp = inall%nf - ny + 1
 ixp = inall/nf			;this is always positive
 iq = iyp lt 0
 iyp = abs(iyp)
 ind = ixp + ny*iyp
 inp = abs(sieve(ind, iq))
 if isarray(inp) then {
 ;this is the "lower rhs" zone
 xq = .5*( a(inp) - b(inp))
 a(inp) = xq
 b(inp) = - xq
 xq = .5*( c(inp) + d(inp))
 c(inp) = xq
 d(inp) = xq
 }

 ;then the other half
 iq = not(iq)
 inp = abs(sieve(ind, iq))
 if isarray(inp) then {
 ;this is the "upper rhs" zone
 xq = .5*( a(inp) + b(inp))
 a(inp) = xq
 b(inp) = xq
 xq = .5*( c(inp) - d(inp))
 c(inp) = xq
 d(inp) = -xq
 }


 ;finally, need to check for edges, note how x and y are switched
 iq = (ixp eq 0) or (ixp eq (ny-1)) or (iyp eq 0) or (iyp eq (nx-1))
 inp = abs(sieve(ind, iq))
 if isarray(inp) then {
 a(inp) = 0.0
 c(inp) = 0.0
 d(inp) = 0.0
 } 

 sc2db, y, a, b, c, d
 
 endblock
 ;=====================================================================