function [kx,ky,kz,ksquare]=makek(M,N,P);

% makes MxNxP matrices holding the components of the vector k
% in the corresponding position in the MxNxP grid.
% The first position in each row, in the case of kx, for example, 
% holds 0, the next 1, etc. up to N/2-1, the next holds -N/2,
% the next -N/2+1, up to -1.  These assignments correspond
% to Matlab's implementation of the FFT. 

ktemp=kron([0:N/2-1 -N/2:-1],ones(M,1));
for i=1:P
   kx(:,:,i)=ktemp;
end

ktemp=kron(ones(1,N),[0:M/2-1 -M/2:-1]');
for i=1:P
   ky(:,:,i)=ktemp;
end

ktemp=[0:P/2-1 -P/2:-1];
for i=1:P
   kz(:,:,i)=ktemp(i)*ones(M,N);
end

ksquare=kx.*kx+ky.*ky+kz.*kz;