  function P2Sol()
% Epicycle fit
  close all

% Get Earth and Mars Data...
  N = 2000;
  tVals = linspace(0,N,N+1)';
  [xE,yE] = Earth(tVals);
  [xM,yM] = Mars(tVals);
  
% A typical plot to build intuition...
  figure
  del = zeros(N,1);
  for k=1:N
      del(k) = Sine([xM(k)-xE(k);yM(k)-yE(k)],[xM(k+1)-xE(k+1);yM(k+1)-yE(k+1)]);
  end
  plot(del)
  title('How much Mars Advances Each Night (radians)')
 
% A typical pencil-and-paper reduction...
% Suppose the two planets moved uniformly in circular orbits. Then
% the Earth-to-Mars vector [x(t);y(t)] is defined by:
%           x(t) = 228*cos((2pi/687)t) - 149.5*cos((2pi/365)*t)
%           y(t) = 228*sin((2pi/687)t) - 149.5*sin((2pi/365)*t)
% So here is a pretty good guess for a 2-circle epicycle;

  r = [ 1 -149.5/228]; P = [687 365];
  
% Let's look at the error...
  [x,y] = Epicycle(r,P,tVals);
  figure
  v = SineVec(xE,yE,xM,yM,x,y);
  plot(v)
  title(sprintf('phi = %10.5f',v'*v))
  
% Let's use the third circle to compensate for the shortcomings
% of our current (roughly determined) 2-circle model.
% Since we're not sure if r(3) should be negative, initialize it
% to zero. And reason that the period of the 3rd disk is half that of
% the second disk, i.e., P(1)/P(2) \approx P(2)/P(3)

  tau0 = [P(1) r(2) P(2) 0 P(2)/2];
  options = optimset('MaxFunEvals',10000,'MaxIter',10000);
  tauStar = fminsearch(@(tau) MyF(tau,xE,yE,xM,yM,tVals),tau0);
% Display the error
  r = [1 tauStar(2) tauStar(4)];
  P = [tauStar(1) tauStar(3) tauStar(5)];
  
% Some solutions and the corresponding value of phi...
  % r = [1  -.65083  -.02266];  P = [684.495  366.711 184.551];   %3.03
  % r = [1  .9079  -1.5902];  P = [369.6371  684.3998  368.4409];  %3.74
  % r = [1  -.6575 -.0810];  P = [686.9254  365.2598  1.12*10^14];  %.007
  % r = [ 1 -.021976 -1.531658]; P = [366.596562 317.040105 684.508767]; %3.46
  % r = [1 -.65581  1.04*10^-5]; P = [ 684.461  366.486  17.173]; %3.78
  % r = [1 -116.1  176.1]; P = [ 384.851452  366.59212  684.43]; %3.77
  % r = [ 1 -.021976  -1.531]; P = [ 366.596  317.04  684.50]; % 3.46
  r = [1 .31 .8]; P = [ 1.1  1.89  3.22]
   
  [x,y] = Epicycle(r,P,tVals);
  figure
  v = SineVec(xE,yE,xM,yM,x,y);
  plot(v)
  title(sprintf('phi = %10.5f',v'*v))
  xlabel(sprintf('P(1) = %6.3f  r(2) = %6.3f  P(2) = %6.3f  r(3) = %6.3f  P(3) = %6.3f',tauStar))
  

  function mu = MyF(tau,xE,yE,xM,yM,tVals)
% The function phi...
  r = [1 tau(2) tau(4)];
  P = [tau(1) tau(3) tau(5)];
  [x,y] = Epicycle(r,P,tVals);
  z = SineVec(xE,yE,xM,yM,x,y);
  mu = z'*z;


  function s = Sine(u,v)
% u and v are 2-vectors.
% s = crossproduct(u,v)/(norm(u)*norm(v))  
  s = (u(1)*v(2) - u(2)*v(1))/(norm(u)*norm(v));
  
  function [xVals,yVals] = Epicycle(r,P,tVals)
% r and P are m-vectors.
% tVals is an N-vector
% Returns in the column vectors xVals and yVals the
% value of the functions
%
%      x(t) = r(1)*cos(2*pi*t/P(1)) + ... + r(m)*cos(2*pi*t/P(m))
%      x(t) = r(1)*cos(2*pi*t/P(1)) + ... + r(m)*cos(2*pi*t/P(m))
%
% at t = tVals(1),...,tVals(N).

  N = length(tVals);
  xVals = zeros(N,1);
  yVals = zeros(N,1);
  m = length(P);
  for k=1:m
      xVals = xVals + r(k)*cos((2*pi/P(k))*tVals);
      yVals = yVals + r(k)*sin((2*pi/P(k))*tVals);
  end
  
  function z = SineVec(xE,yE,xM,yM,x,y)
% The input vectors are column N-vectors
% z is a column n-vector with the property that z(i) is the sine of
% the angle between [x(i);y(i)] and [xM(i)-xE(i);yM(i)-yE(i)].
  N = length(xE);
  z = zeros(N,1);
  for k=1:N
      u = [xM(k)-xE(k);yM(k)-yE(k)];
      v = [x(k);y(k)];
      z(k) = Sine(u,v);
  end
  
  function [xE,yE] = Earth(t)
% Earth is located at (xE,yE) at time t.
  PE = 147.098; AE = 152.097; rhoE = 365.24;
  theta = 2*pi*t/rhoE;
  xE = (PE-AE)/2 + ((PE+AE)/2)*cos(theta);
  yE = sqrt(PE*AE)*sin(theta);
     
  function [xM,yM] = Mars(t)
% Mars is located at (xM,yM) at time t.
  PM = 206.669; AM = 249.209; rhoM = 686.971;
  theta = 2*pi*t/rhoM;
  xM = (PM-AM)/2 + ((PM+AM)/2)*cos(theta);
  yM = sqrt(PM*AM)*sin(theta);
 

  
  
 
  
  

  
  
  
   