function Demo(t1,t2,P1,theta1,P2,r2,theta2,P3,r3,theta3)
% Displays the epicycle predictions of Mars from time t1 to t2 and plots
% the error.
% t1 an integer that indicates the starting day.
% t2 an integer that specifies the ending day. Assume t1<t2.
% Circle 1 has period P1 (days), radius 1, and initial angle theta1 (degrees).
% Circle 2 has period P2 (days), radius r2, and initial angle theta2 (degrees).
% Circle 3 has period P3 (days), radius r3, and initial angle theta3 (degrees).

% Useful values..
c = cos(linspace(0,2*pi));
s = sin(linspace(0,2*pi));
del = .05;

hold on
% The first circle...
plot(c,s,'w')
% Its center...
fill(del*c,del*s,'w')
% The celestial sphere...
r = 2.5;
plot(r*c,r*s,'b')
% The window...
d = max(1.5*(1+r2+r3),3);
hold on
axis([-d d -d d])
axis square
set(gca,'color','k')
set(gcf,'position',[100 100 800 800],'color','k')
shg
sinErr = zeros(t2-t1+1,1);
k = 0;
for t=t1:t2
% Display the situation at time t.
% First compute the location of the three "dots"
% (x1,y1) is the location of the dot on the first circle...
mu1 = 2*pi*t/P1 + theta1*pi/180; x1 = cos(mu1); y1 = sin(mu1);
% (x2,y2) is the location of the dot on the second circle...
mu2 = 2*pi*t/P2 + theta2*pi/180; x2 = x1+r2*cos(mu2); y2 = y1+r2*sin(mu2);
% (x3,y3) is the location of the dot on the third circle...
mu3 = 2*pi*t/P3 + theta3*pi/180; x3 = x2+r3*cos(mu3); y3 = y2+r3*sin(mu3);
% Compute their projections (x1p,y1P), (x2P,y2P), and (x3P,y3P) onto the "celestial sphere"...
d1 = r/sqrt(x1^2+y1^2); x1P = x1*d1; y1P = y1*d1;
d2 = r/sqrt(x2^2+y2^2); x2P = x2*d2; y2P = y2*d2;
d3 = r/sqrt(x3^2+y3^2); x3P = x3*d3; y3P = y3*d3;
% (x,y)= the position Mars on the "celestial sphere" as seen from Earth...
[xE,yE] = Earth(t);
[xM,yM] = Mars(t);
x = xM-xE;
y = yM-yE;
d = r/sqrt(x^2+y^2);
x = x*d;
y = y*d;
% Compute the sine of the discrepancy angle...
u = [x;y]; u = u/norm(u); v = [x3;y3]; v = v/norm(v);
k = k+1;
sineErr(k) = abs(u(1)*v(2)-u(2)*v(1));

% Now the graphics...
h1 = plot(x1+r2*c,y1+r2*s,'w'); % circle 2
h2 = plot(x2+r3*c,y2+r3*s,'w'); % circle 3
h3 = plot([0 x1 x2 x3],[0 y1 y2 y3],'w'); % The radii
h4 = fill(x1+del*c,y1+del*s,'y'); % 1st circle dot
h5 = fill(x2+del*c,y2+del*s,'c'); % 2nd circle dot
h6 = fill(x3+del*c,y3+del*s,'m'); % 3rd circle dot
h10 = fill(x+1.5*del*c,y+1.5*del*s,'r');
h7 = fill(x1P+del*c,y1P+del*s,'y'); % 1st circle dot projection
h8 = fill(x2P+del*c,y2P+del*s,'c'); % 2nd circle dot projection
h9 = fill(x3P+del*c,y3P+del*s,'m'); % 3rd circle dot projection
h11 = text(-2.7,-2.7,sprintf('t = %4d',t),'fontsize',26,'color','y');
h12 = text(-2.7,-3.2,sprintf('|sin(errorAngle)| = %5.3f ',sineErr(k)),'fontsize',26,'color','y');
if t==t1
% take a look at the situation on the first day...
pause(2)
else
% Subsequent situations displayed briefly so we get a "movie"
pause(.001)
end
if t<t2
% Erase the "moving parts" unless the simulation has ended...
delete(h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11,h12)
end
end