function ErrorPlot(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).
% If wantSimulation == 0, then the simulation is suppressed

r = 2.5;
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));
end
plot((t1:t2)',sineErr)
xlabel('t','fontsize',14)
ylabel('| sin( error angle at t ) |','fontsize',14)
title(sprintf('t1 = %4d t2 = %4d',t1,t2),'fontsize',14)
shg