function A8()
% Illustrates a 3-circle Epicycle model.
close all
% First circle...
P1 = 300;
theta1 = 120;
% Second circle...
P2 = 200;
r2 = .5;
theta2 = 200;
% Third circle...
P3 = 130;
r3 = .2;
theta3 = 10;

t1 = 0;
t2 = 2200;

Demo(t1,t2,P1,theta1,P2,r2,theta2,P3,r3,theta3)

function Demo(t1,t2,P1,theta1,P2,r2,theta2,P3,r3,theta3)
% Displays the epicycle predictions of Mars from time t1 to t2.
% 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 = 1.5*(1+r2+r3);
hold on
axis([-d d -d d])
axis square
set(gca,'color','k')
set(gcf,'position',[100 100 800 800],'color','k')
shg
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;
% 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
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
h10 = fill(x+1.5*del*c,y+1.5*del*s,'r');
if t==t1
% take a look at the situation on the first day...
pause
else
% Subsequent situations displayed briefly so we get a "movie"
pause(.01)
end
if t<t2
% Erase the "moving parts" unless the simulation has ended...
delete(h1,h2,h3,h4,h5,h6,h7,h8,h9,h10)
end
end
hold off

function [xE,yE] = Earth(t)
% (xE,yE) is the location of Earth at time t (days).
P = 147.1; A = 152.1; phi = 102*pi/180; Y = 365.24;
x = (P-A)/2 + ((P+A)/2)*cos((2*pi/Y)*t);
y = sqrt(P*A)*sin((2*pi/Y)*t);
c = cos(phi); s = sin(phi);
xE = c*x - s*y;
yE = s*x + c*y;

function [xM,yM] = Mars(t)
% (xM,yM) is the location of Mars at time t (days).
P = 206.6; A = 249.22; phi = 336*pi/180; Y = 686.93;
x = (P-A)/2 + ((P+A)/2)*cos((2*pi/Y)*t);
y = sqrt(P*A)*sin((2*pi/Y)*t);
c = cos(phi); s = sin(phi);
xM = c*x - s*y;
yM = s*x + c*y;

function [P1,theta1,P2,r2,theta2,P3,r3,theta3,err] = ThreeCircle(t1,t2)
% t1 and t2 are integers that satisfy 0<=t1 and t1+7<=t2.
% Let sigma(t) be the sine of the angle between the vector [xM(t)-xE(t); yM(t)-yE(t)]
% where
% [xM(t);yM(t)] is the location of Mars at time t
% [xE(t);yE(t)] is the location of Earth at time t.
% and the vector [x(t);y(t)] where
% x(t) = cos(2*pi*t/P1 + theta1*pi/180) +
% r2*cos(2*pi*t/P2 + theta2*pi/180) +
% r3*cos(2*pi*t/P3 + theta3*pi/180)
% y(t) = sin(2*pi*t/P1 + theta1*pi/180) +
% r2*sin(2*pi*t/P2 + theta2*pi/180) +
% r3*sin(2*pi*t/P3 + theta3*pi/180)
% The epicycle parameters P1,theta1,P2,r2,theta2,P3,r3,theta3 are chosen to
% minimize sigma(t1)^2 + sigma(t1+1)^2 + ... + sigma(t2)^2.
% err is the maximum of |sigma(t1)|,...,|sigma(t2)|