function power_xy2r, fftxy, kr = kr

if (n_params() eq 0) then begin
    Message,'NAME: POWER_XY2R.PRO',/info
    Message,'',/info
    Message,'PURPOSE: This function sums power in radial annuli of',/info
    Message,'  a 2-D array of (k_x,k_y), to compute power vs k_r.',/info
    Message,'',/info
    Message,'USAGE: fftr = power_xy2r(fftxy, kr = kr)',/info
    Message,'',/info
    Message,'INPUT:',/info
    Message,'  FFTXY: complexarr(nx,ny) of FFT power in (k_x, k_y)',/info
    Message,'',/info
    Message,'OUTPUT:',/info
    Message,'  FFTR: complexarr(nkr) of FFT power in (k_r)',/info
    Message,'',/info
    Message,'KEYWORDS:',/info
    Message,'  KR: Set to return 1-D array of radial wavenumbers',/info
    Message,'',/info
    Message,'BLAME: Written 01 May 2007, B.T. Welsch, UCB/SSL',/info
    return,0
endif

nx = n_elements(fftxy(*,0))
ny = n_elements(fftxy(0,*))

; Define 1-D frequency arrays as IDL does -- for even and odd cases.
;====================================================================

; First for X
if ((nx)mod(2) eq 0) then kx1d = $
[0,(findgen(nx/2) + 1)/(nx),-reverse(findgen(nx/2-1) + 1)/(nx)] $
else kx1d = $
[0,(findgen(nx/2.-.5) + 1)/(nx),-reverse(findgen(nx/2.-.5) + 1)/(nx)]

; Next, for Y.
if ((ny)mod(2) eq 0) then ky1d = $
[0,(findgen(ny/2) + 1)/(ny),-reverse(findgen(ny/2-1) + 1)/(ny)] $
else ky1d = $
[0,(findgen(ny/2.-.5) + 1)/(ny),-reverse(findgen(ny/2.-.5) + 1)/(ny)]

; Make 1-D array of k_r values using max([nx,ny])*sqrt(2)
if (nx ge ny) then kr1d = kx1d(where(kx1d ge 0)) else $
  kr1d = ky1d(where(ky1d ge 0))
sort_kr = kr1d(sort(kr1d))
 n_kr = n_elements(sort_kr)
fftr = complexarr(n_kr)

; Define 2-D frequency arrays, in X, Y, and R
kx2d = kx1d#replicate(1,ny)
ky2d = replicate(1,nx)#ky1d
kr2d = sqrt(kx2d^2 + ky2d^2)

for i=0,n_kr-2 do begin
    inbin = where((kr2d ge sort_kr(i)) and $
                  (kr2d lt sort_kr(i+1)), $
                  n_inbin)
    if (n_inbin ne 0) then fftr(i) = total(fftxy(inbin)*conj(fftxy(inbin)))
endfor

kr = sort_kr
return,fftr

end


pro freq_filter, data_in, data_out, dt=dt, dx=dx, $
                 cutoff=cutoff, readcut=readcut, $
                 kwnoplot=kwnoplot, kwarr=kwarr, kwplot=kwplot, $
                 saveplot=saveplot, hipass=hipass, verbose=verbose

if (n_params() eq 0) then begin
    Message,'NAME: FREQ_FILTER.PRO',/info
    Message,'',/info
    Message,'PURPOSE: Zeroes-out low (default) or high freqs. in',/info
    Message,'  (nx,ny,nt) data array. Can compute k vs. omega if nx=ny.',$
      /info
    Message,'',/info
    Message,'USAGE: freq_filter, data_in, data_out, cutoff=5.,/kwplot',/info
    Message,'',/info
    Message,'ARGUMENTS: DATA_IN: fltarr(nx,ny,nt), must be regular cadence', $
      /info
    Message,'           DATA_OUT: fltarr(nx,ny,nt) filtered data',/info
    Message,'',/info
    Message,'KEYWORDS: DT: Data cadence, seconds assumed; default=1',/info
    Message,'  DX: Pixel size, default=1 (k vs. omega needs square pixels)', $
      /info
    Message,'  HIPASS: Set to 1 to override DEFAULT lo-pass filtering.',/info
    Message,'  CUTOFF: zeros freqs. above CUTOFF (below CUTOFF if /HIPASS)',$
      /info
    Message,'  READCUT: Prompts for CUTOFF from keyboard; useful',/info
    Message,'    if /KWPLOT is set, for choosing CUTOFF at run-time',/info
    Message,'  KWPLOT: Computes k-omega array, makes contour plot. ',/info
    Message,'  SAVEPLOT: Saves contour plot to file "kwcontours.ps".',/info
    Message,'  KWNOPLOT: Computes k-omega array, nothing is plotted. ',/info
    Message,'  KWARR: Returns k-omega array if either /KWNOPLOT or /KWPLOT.',$
      /info
    Message,'  VERBOSE: Set to print progress, other info.',/info
    Message,'',/info
    Message,'BLAME: Written 01 May 2007, B.T. Welsch, UCB/SSL',/info
    return
endif

nx = n_elements(data_in(*,0,0))
ny = n_elements(data_in(0,*,0))
nt = n_elements(data_in(0,0,*))

if not(keyword_set(dt)) then dt = 1. ; temporal spacing
if not(keyword_set(dx)) then dx = 1. ; grid spacing

if (keyword_set(verbose)) then print,'dt is:',dt
if (keyword_set(verbose)) then print,'dw is:',dw

data_out = fltarr(nx,ny,nt)

; Generate frequency array as IDL does.
;===================================================
if ((nt)mod(2) eq 0) then $
  warr0 = [0,(findgen(nt/2) + 1)/(nt), $ ; for nt = even case
           -reverse(findgen(nt/2-1) + 1)/(nt)] else $
  warr0 = [0,(findgen(nt/2.-.5) + 1)/(nt), $ ; for nt = odd case
           -reverse(findgen(nt/2.-.5) + 1)/(nt)]

warr = warr0/dt                 ; scale to units of dt


; If either /KWPLOT or /KWNOPLOT, will generate a k-omega diagram.
if (keyword_set(kwplot) or keyword_set(kwnoplot)) then begin 

    if keyword_set(verbose) then $
      print,'Generating k-omega diagram.'

    fft_xy = fft(data_in(*,*,0))  ; spatial FFT of 0th time slice
    fft_r = power_xy2r(fft_xy, kr = kr) ; Convert spatial (kx,ky) to (kr)
    nkr = n_elements(kr)        ; number of elements 1-D array of kr

    ktarr = complexarr(nkr,nt)  ; define 2-D array of (kr, time)
    ktarr(*,0) = fft_r          ; 1st row, for 0th time slice

    for i=1,nt-1 do begin ; compute power vs. kr for other time slices
        fft_xy = fft(data_in(*,*,i)) 
        fft_r = power_xy2r(fft_xy)    
        ktarr(*,i) = fft_r ; now fill in ith rows

        if (keyword_set(verbose)) then $
          print,'Computing k_r in slice',i

    endfor

    ; Next, FFT in time, to get k vs. omega
    kwarr = fft(ktarr,-1,dimension=2)

    ; Locate (+/-) modes.
    where_posw = where(warr gt 0, n_posw) ; if nt is even, n_posw=n_negw+1
    where_negw = reverse(where(warr lt 0, n_negw))

    ; Define array for combined power in negative & positive freqs.
    power_1side = fltarr(nkr, n_posw+1) ; zero-mode is DC, include it

    ; Sum power in (+/-) modes
    ;===============================
    osc_power = abs(kwarr)^2    ; compute modulus^2 prior to summing
    power_1side(*,1:n_posw) = osc_power(*,where_posw)
    power_1side(*,1:n_negw) = power_1side(*,1:n_negw) + $ 
      osc_power(*,where_negw) ; for nt even, negw is missing highest pos freq 
    power_1side(*,0) = osc_power(*,0)

    warr_1side = [0,warr(where_posw)] ; 1-D array of freqs. in power_1side

    if keyword_set(kwplot) then begin

        if keyword_set(saveplot) then begin
            set_plot,'ps'
            device,filename='kwcontours.ps'
        endif

        nzpower = where(power_1side ne 0, n_nzpower)
        if (n_nzpower eq 0) then print,fail_variable ; No power! data_in empty?
        max_lev = ceil(max(alog10(power_1side(nzpower))))
        levels = 10^(findgen(max_lev-1)+1) ; powers of 10, exclude lowest

        contour,power_1side,kr,warr_1side,chars=1.5, $
          /follow,levels=levels,$
          xtit='Wavenumber, |k|', $
          ytit='Frequency, !4x!x',$
          tit='1-Sided Power, k vs. !4x!x'

        if keyword_set(saveplot) then set_plot,'x'

        if (keyword_set(verbose) and keyword_set(saveplot)) then $
          print,'Saving k-omega plot.'

    endif ; ends if-KWPLOT check for 
            
endif ; ends if-KWNOPLOT or KWPLOT check

if keyword_set(readcut) then read,cutoff,prompt= $
  'Please enter cutoff frequency:'

data_xyw = fft(data_in,-1,dimension=3) ; temporal FFT of input data

if not(keyword_set(hipass)) then begin ; Low Pass Filter

    if keyword_set(verbose) then print,'Low-pass filtering.'

    hi_w = where(abs(warr) gt cutoff, n_hiw)
    if (n_hiw eq 0) then print,'No frequencies above chosen cutoff.' $
    else data_xyw(*,*,hi_w) = 0.    ; zeros out both Im() & Re() parts!
    data_out_complex = fft(data_xyw,1,dimension=3)

endif else begin ; High Pass Filter

    if keyword_set(verbose) then print,'High-pass filtering.'

    lo_w = where(abs(warr) lt cutoff, n_low)
    if (n_low eq 0) then print,'No frequencies below chosen cutoff.' $
    else data_xyw(*,*,lo_w) = 0.         ; zeros out both Im() & Re() parts!
    data_out_complex = fft(data_xyw,1,dimension=3)

endelse

data_out = real_part(data_out_complex)

end