  function ShowMyError
% Take a look at the error of your solution...
  close all

% Get Earth and Mars Data...
  N = 1000;
  tVals = linspace(0,N,N+1)';
  [xE,yE,xM,yM] = SetUp(tVals);
  
% Here are some student solutions. (Some had the same optimum values for r and P).
% Remove the "%" next to yours to see error plot.
%                                                                         Phi   
   r = [1 -.63648  -.12816];    P = [691.946  367.4544  3.7*10^17]; %   1.476 
  % r = [1 -.6365  -.1282];      P = [691.946  367.4548  2.55*10^10];%   1.476 
  % r = [1 -.63735  -.1395];     P = [694.8273  368.309  3*10^8];    %   1.78
  % r = [1 7.2  -4.4];           P = [ 1221.7  682.6  367.7];        %   2.471
  % r = [1  -.5968  -.0681];     P = [ 685.7895 372.1769  241.2704]; %   2.788
  % r = [1  -.2818 -.2985];      P = [687.9  373.1  373.1];          %   8.394 
  % r = [1 -.61714  -.09164];    P = [.666026 .6654798 .667846];     %  20.192
  % r = [1 .4948  .48226];       P = [ 757.18  500.284  2093.45];    %  47.9
  % r = [1  2.1  1.0];           P = [ 2093.5  757.2  500.3];        %  49.3 
  % r = [1 -.6282  -.1426];      P = [694.0522  367.0190  1.001];    %  59.88 
  % r = [1   .4915     .0671];   P = [ 693.88  3821      493.78];    %  70.72
  % r = [1 .45789 .1828];        P = [743.5629  535.1864  535.1863]; % 162.41 
  % r = [1 .56 .29];             P = [688.44  4170.15  1];           % 166.127
  % r = [1 .2939  .241];         P = [ 693.7    298.76  339];        % 272.96
  % r = [1  -.61723    -.0915];  P = [.666    .285      .667];       % 331.6 
  % r = [139.49  229.83 0.4496]; P = [373.13  697  126.3427];        % 533.1 
  
% 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))
 


  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
  
  
     
 

  
  
 
  
  

  
  
  
   