restore,'data/bbk.sav'

s=size(tb1) & nx=s[1] & ny=s[2] & nz=s[3]
if n_elements(size) eq 0 then size=2.*!pi
k0=2.*!pi/size

rad=sqrt(x1^2+y1^2+z1^2)                        ; wave vector radius

KK=0
ta1=shift(tb1,0,0,KK)  ; shift by 1
ta2=shift(tb2,0,0,KK)  ; shift by 1
ta3=shift(tb3,0,0,KK)  ; shift by 1

imax=min([nx,ny,nz]/2)                       ; cut corners
;if keyword_set(nocorners) then $
; imax=min([nx,ny,nz]/2) $                      ; cut corners
;else $
; imax=fix(max(rad))                            ; include corners
spectrum=fltarr(imax)

for i=0,imax-1 do begin                         ; cut the corners?
  ii=where((rad ge i-.5) and (rad lt i+.5))     ; indices into shell
  tmp=float(ta1(ii)*conj(tb1(ii))+$
            ta2(ii)*conj(tb2(ii))+$
            ta3(ii)*conj(tb3(ii)))
  spectrum(i)=total(tmp)                      ; sum over shell
endfor

k0=1.
wavenumbers=indgen(imax)*k0                     ; dimensional wavenos
spectrum=spectrum/k0                            ; normalize to k0

plot_oo,wavenumbers(1:*),spectrum(1:*)

end