function spiral(r1, r2, b, point)
% SPIRAL Spiral plot
%	Usage SPIRAL(R1, R2, B, POINT)
%	R1 = Integer radius of big circle 
%	R2 = Integer radius of small circle 
%	B  = Distance between pencil hole and center of small circle
%	POINT = No. of points per rotation around center of big circle
%
%	Something to try on:
%		spiral(10, 7, 0);	% a big circle
%		spiral(10, 7, 1);	% dount
%		spiral(10, 7, 3);	% rounded big petals
%		spiral(10, 7, 4);	% small petal
%		spiral(10, 7, 7);	% close to a star
%
%		spiral(10, 13, 1);	% another donut
%		spiral(10, 13, 6);	% outside-in patal
%		spiral(10, 13, 8);	% outside-in patal
%
%		spiral(5, 3, 12);	% a fat star

%	Roger Jang, 1991

if nargin < 4, point = 100; end
if nargin < 3, b = 31; end
if nargin < 2, r2 = 65; end
if nargin < 1, r1 = 132; end

p = r1;
q = r2;

while p ~= q*floor(p/q)
	r = p - q*floor(p/q);
	p = q;
	q = r;
end

% q is the gcd of r1 and r2
rotation = r2/q;

theta1 = linspace(0, 2*pi*rotation, point*rotation);
theta2 = r1*theta1/r2;
x = b*cos(theta1 - theta2) + (r1 - r2)*cos(theta1);
y = b*sin(theta1 - theta2) + (r1 - r2)*sin(theta1);
%plot(x,y)
axis image;
mycomet(x, y, 0.01);