function [cx,cy,cz]=nonlinear(bx,by,bz,kx,ky,kz,ksquare);

% function [cx,cy,cz]=nonlinear(bx,by,bz,kx,ky,kz,ksquare);
% All arguments are MxNxP matrices. 
% One must first have run [kx,ky,kz,ksquare]=makek(M,N,P);

[M,N,P]=size(bx);
ksquare(1,1,1)=1;   % temporary, to avoid a warning message
i=sqrt(-1);
ax=i*(ky.*bz-kz.*by)./ksquare;
ay=i*(kz.*bx-kx.*bz)./ksquare;
az=i*(kx.*by-ky.*bx)./ksquare;
ksquare(1,1,1)=0;   % undo temporary fix above
ax(1,1,1)=0;        % average x-velocity = 0
ay(1,1,1)=0;        % average y-velocity = 0
az(1,1,1)=0;        % average z-velocity = 0
dx=fft3(ifft3(ay).*ifft3(bz)-ifft3(az).*ifft3(by))*M*N*P;
dy=fft3(ifft3(az).*ifft3(bx)-ifft3(ax).*ifft3(bz))*M*N*P;
dz=fft3(ifft3(ax).*ifft3(by)-ifft3(ay).*ifft3(bx))*M*N*P;
cx=i*(ky.*dz-kz.*dy);
cy=i*(kz.*dx-kx.*dz);
cz=i*(kx.*dy-ky.*dx);