function [F,u,v] = rdft2(f,x,y)
%
%  RFFT2
%    Roger's slow 2-D DFT function.
%  [F,u,v] = RFFT2(f,x,y) returns in F the 2-D Discrete
% Fourier transform of the data in f defined over the
% domain in x and y which are matrices created by the 
% MESHDOM command.  The matrices u and v are spatial
% frequency variables and may be omitted.
%
%  Note that this command is written for academic purposes
%  and is much slower than the FFT and FFT2 commands 
%  which use a fast radix-2 type transform.
%                                                  RJM  3/22/95

[n,m]=size(x);

[p,q]=meshdom(0:1:n-1,0:1:m-1);

i=sqrt(-1);
pi=3.141592653589793284626433;

dy=(y(n,1)-y(1,1))/n;
dx=(x(1,m)-x(1,1))/m;

for j = 1:m,
  for k= 1:n,
    
  Fmat=f.*exp(-i*2*pi*((k-1)/n*ones(n,m).*p+(j-1)/m*ones(n,m).*q));    
    % Above line performs all multiplications
  F(j,k) = sum(sum(Fmat));
    % Above line adds up all the multiplications
end
end

% Now F contains the 2-D DFT.

%Here we perform Quadrant-Flipping

a1=ceil(n/2);
b1=ceil(m/2);
a2=ceil(n/2+1);
b2=ceil(m/2+1);

F = [ F(a2:n,b2:m)   F(a2:n,1:b1);
      F(1:a1,b2:m)   F(1:a1,1:b1)];