pro power2d,bx,by,bz,kk=kk,ff=ff,asp=asp,noplo=noplo
;
if n_params(0) eq 0 then begin
  print,'pro power2d,bx,by,bz,kk=kk,ff=ff,asp=asp'
  return
endif
if n_elements(asp) eq 0 then asp=[3,3]
;
bx=reform(bx)
by=reform(by)
bz=reform(bz)
;
;  makes a 2-D power spectrum
;
s=size(bx)
mx=s(1)
my=s(2)
print,mx,my
;
;  compute k-array
;
dkx=2.*!pi/asp(0)
dky=2.*!pi/asp(1)
kx=findgen(mx)*dkx
ky=findgen(my)*dky
kk=sqrt(kx^2#(ky-ky+1.)+(kx-kx+1.)#ky^2)
kk=kk(1:mx/2,1:my/2)
;
;  work space
;
t1=complexarr(mx,my)
t2=complexarr(mx,my)
;
;  add together the contributions from the x,y,z components
;
for l=0,mx-1 do t1(l,*)=fft(bx(l,*),-1)
for m=0,my/2 do t2(*,m)=fft(t1(*,m),-1)
ff=abs(t2(1:mx/2,1:my/2))^2
;
for l=0,mx-1 do t1(l,*)=fft(by(l,*),-1)
for m=0,my/2 do t2(*,m)=fft(t1(*,m),-1)
ff=ff+abs(t2(1:mx/2,1:my/2))^2
;
for l=0,mx-1 do t1(l,*)=fft(bz(l,*),-1)
for m=0,my/2 do t2(*,m)=fft(t1(*,m),-1)
ff=ff+abs(t2(1:mx/2,1:my/2))^2
;
;  multiply by kk to get the spectrum as if it was integrated over shells
;  (although here we do not integrate, but plot individual points).
;
ff=ff*kk
if n_elements(noplo) eq 0 then plot,alog10(kk),alog10(ff),ps=3,xtitl='!6log!d10!n!8k!6',ytitl='!6log!d10!n!8E!dK!n!6(!8k!6)'
;
end