%   Example code for "Generating Rectangular Coordinates 
%   in Polar Coordinate Order" article in Nov. 2005 issue
%   of "DSP Tips & Tricks" colunm of the IEEE Signal 
%   Processing magazine.
%
%   Article Author:  Prof. Charles M. Rader

clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ****** FIRST ROUTINE   ************** 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = 5;
L = 5;

a=0;b=1;c=1;d=L; %initialization
m=[b d];
n=[a c]; 
while ~((c==M)&(d==1)) % end test
   Z = floor((L+b)/d);
   e = Z*c-a;
   f= Z*d-b;
   if e<=M
      % report pair (e,f)
      m = [m  f];
      n = [n  e];
   end
   b=d;d=f;a=c;c=e;
end 
figure(1)
plot(m,n,'-.o','markersize',10)
axis([0 L 0 M]), grid on
set(gca, 'FontSize', 14); set(gca, 'FontName', 'Arial');  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ****** SECOND ROUTINE   ************** 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=0; b=1; c=1; d=L; %initialization
% report pair a,b
% report pair c,d 
while ~((c==1)&(d==1)) %end test
   Z = floor((L+b)/d);
   e = Z*c-a;
   f= Z*d-b; 
   % report pair e,f

         e,f

   b=d; d=f; a=c; c=e;
end