function [P1,theta1,P2,r2,theta2,P3,r3,theta3,err] = ThreeCircleCVL(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)|

options = optimset('maxiter',2000,'tolX',.0001,'tolFun',.0001);

% Compute the R unit vectors once and for all. Each f-evaluation required
% by lsqnonlin uses these.
R = zeros(2,5001);
tVals = 0:5000;
for k=1:5001
% R(:,k) = unit vector that points to Mars as seen from Earth at time
% k-1...
[xM,yM] = Mars(tVals(k));
[xE,yE] = Earth(tVals(k));
R(:,k) = [xM-xE;yM-yE];
R(:,k) = R(:,k)/norm(R(:,k));
end

% Generates a good starting guess by solving the t1A = 500, t2A = 5000
% problem.
% From the plot on the handout, the "down spikes", i.e., the instances of
% maximum retrograde motion are about 782 days apart and the first one is roughly
% at t = 500. (Use ShowDailyAngle(t1,t2) for various t1 and t2 to see this.)
t1A = 500;
t2A = 5000;
% Location of Earth and Mars at t = 500...
[xM,yM] = Mars(t1A);
[xE,yE] = Earth(t1A);
% Compute the average radian advance of the "red dot" over [500,5000]...
psi = DailyAngle(t1A,t2A);
psiAve = mean(psi);
% We use the first circle to roughly have the same period of revolution,
% i.e., solve 2*pi/P1 = psiAve for P1
P1 = 2*pi/psiAve;
% Require theta1 to satisfy
% (2*Pi/P1)*t1A + theta1*pi/180 = atan2(yM-yE,xM-xE)
% Thus,
theta1 = (180/pi)*(atan2(yM-yE,xM-xE)- (2*pi/P1)*t1A);
% Assuming P2<P1 and that the (0,0), the yellow dot and cyan dot are colinear,
% what must P2 be if it takes 782 for the dotsto line up again, i.e.,
% (2*pi/P2)*782 = (2*pi/P1)*782 + 2*pi, tus...
P2 = 1/(1/782 + 1/P1);
% At t = 500 (roughly) we have the first down spike. To realize the
% required backward motion we choose theta2 so that second circle dot
% (cyan) is on the ray from (0,0) to the yellow dot, but is INSIDE the
% first circle
theta2 = (180/pi)*(atan2(yM-yE,xM-xE)- (2*pi/P2)*t1A + pi);
% Since the 2*pi/P1 is the forward angular rate of the big circle (about
% .012 radians per day) we need the second circle radius to be big enough
% so that we get a reverse rate of about .006 radians per day.
% Thus we require 2*pi/P1 to approximately equal to 1.5*((2*pi/P2)*r2), i.e,
r2 = (1.5)*(P2/P1);
% We experimented and it pays to have the third circle spin clockwise
% with roughly the same period as the second circle.
P3 = -P2;
r3 = .01;
theta3 = 0;
z0 = [P1 theta1 P2 r2 theta2 P3 r3 theta3]';
% Let's run lsqnonline to improve the fit..
LB = [P1-100 -inf P2-50 r2-.1 -inf -50-P2 0 -inf ]';
UB = [P1+100 inf P2+50 r2+.1 inf +50-P2 r2 inf ]';
z0 = lsqnonlin(@(z) f(z,R,t1,t2),z0,LB,UB,options);
P1 = z0(1);
theta1 = z0(2);
P2 = z0(3);
r2 = z0(4);
theta2 = z0(5);
P3 = z0(6);
r3 = z0(7);
theta3 = z0(8);

LB = [P1-100 -inf P2-50 r2-.1 -inf -50-P2 0 -inf ]';
UB = [P1+100 inf P2+50 r2+.1 inf +50-P2 r2 inf ]';
[zStar,resnorm,residual] = lsqnonlin(@(z) f(z,R,t1,t2),z0,LB,UB,options);
P1 = zStar(1);
theta1 = zStar(2);
P2 = zStar(3);
r2 = zStar(4);
theta2 = zStar(5);
P3 = zStar(6);
r3 = zStar(7);
theta3 = zStar(8);
err = max(abs(residual));

function F = f(z,R,t1,t2)
% Let N = t2-t1+1. This function returns a column N-vector F
% where F(k) is the sin of the angle between Q(:,k) and R(:,k)

% Unpack z...
P1 = z(1); theta1 = z(2);
P2 = z(3); r2 = z(4); theta2 = z(5);
P3 = z(6); r3 = z(7); theta3 = z(8);
% Generate all the Q-vectors and normalize...
tau1 = (2*pi/P1)*(t1:t2) + theta1*pi/180;
tau2 = (2*pi/P2)*(t1:t2) + theta2*pi/180;
tau3 = (2*pi/P3)*(t1:t2) + theta3*pi/180;
Q = [cos(tau1)+r2*cos(tau2)+r3*cos(tau3);...
sin(tau1)+r2*sin(tau2)+r3*sin(tau3)];
for k=1:t2-t1+1
Q(:,k) = Q(:,k)/norm(Q(:,k));
end
F = (Q(1,:).*R(2,t1+1:t2+1) - Q(2,:).*R(1,t1+1:t2+1))';