function epicycloid(a,b,rev)                    %<last updated 8/18/00>
%To generate an epicycloid using a circle of radius a as a base
%and movng a circle of radius b <= a tangent to the base circle
%in a counter clockwise direction.
%The point (a,0) is tracked as the smaller circle rotates.
%
%INPUT:  a = radius of base circle
%        b = radius of outside circle (b <= a)
%        rev = number of revolutions to use (default is 1)
%
%        If either a or b is not specified then the routine sets a = 8, b = 8, rev = 1.
%
%OUTPUT: Graphical construction of an epicycloid.
%
%  By: David R. Hill, Mathematics Dept., Temple Univ.
%      Philadelphia, PA. 19122    Email: dhill001@temple.edu

if nargin==0,a=8;b=8;rev=1;end
if nargin<=2,rev=1;end
angleend=rev*2*pi;
if b > a
   disp('ERROR: b must be less than or equal to a.')
   disp('Try again.')
   return
end

%Showing fixed circle of radius and point of tangency.
ta=0:.01:2*pi;
xa=a*cos(ta);ya=a*sin(ta);
figure
plot(xa,ya,'-k','erasemode','none')
axis((2*(a+b)+2)*[-1 1 -1 1]),axis(axis),axis('square'),hold on
p=a*[1,0];
plot(p(1),p(2),'og','erasemode','none')
bangles=0:.05:2*pi;
C=(a+b)*[cos(0),sin(0)];
bcir=plot(C(1)+b*cos(bangles),C(2)+b*sin(bangles),'.b','markersize',1);
pause(3)
delete(bcir);
%angle increment
%
inc=pi/15;
for t=0:inc:angleend;
   C=(a+b)*[cos(t),sin(t)];
   bcir=plot(C(1)+b*cos(bangles),C(2)+b*sin(bangles),'.b','markersize',1);
   x=(a+b)*cos(t)-b*cos((a+b)/b*(t));
   y=(a+b)*sin(t)-b*sin((a+b)/b*(t));
   plot(x,y,'.r','erasemode','none')
   brad=plot([C(1) x],[C(2) y],'-b');
   pause(1)
   if t<angleend
      tt=t:.01:t+inc;
      xx=(a+b)*cos(tt)-b*cos((a+b)/b*(tt));
      yy=(a+b)*sin(tt)-b*sin((a+b)/b*(tt));
      plot(xx,yy,'-r','erasemode','none')
   end
   delete(bcir),delete(brad)
end
hold off