  function [xcStar,phiStar,maxErr] = BestEccentric(k)
% k is an integer that satisfies 1<=k<=20
% xcStar and phiStar are the optimal eccentric model parameters
% from the standpoint of predicting the first 4k season events.
% maxErr returns the max(|s(1)-shat(1)|,...,|s(4k)-shat(4k)|}
% where s(i) is the observed start time of the i-th event and
% shat(i) is the predicted start time of the i-th event obtained
% by the eccentric model with parameters xcStar and phiStar.

% Assemble the relevant start time data in s
S = StartTimes(); 
S = S(1:4*k);
% A rendom start vector works just fine...
x0 = rand(2,1);
% Turn off lsqnonlin output and set termination criteria...
options = optimset('Display','off','TolFun',10^-6,'TolX',10^-6);
% The call...
[x,resNorm,resVec] = lsqnonlin(@(x) MyF(x,S),x0,[],[],options);
% Deciphering x and looking at the maximum error among all the
% predicted start times...
xcStar =  x(1);
phiStar = x(2);
maxErr = max(abs(resVec));

  function R = MyF(x,S)
% The vector of discrepancies between S and the start time predictions
% based on the current choice of eccentric model parameters...
xc  = x(1);
phi = x(2);
%aveYear = 365.2421875;
aveYear = 365.25;       % AH

% Compute the season season angles according to the eccentric model...
c = cos(phi); s = sin(phi); tau = -c*xc + sqrt(1-(s*xc)^2);
thetaSpring = MyAtan(tau*s,xc+tau*c);
c = cos(phi+pi/2); s = sin(phi+pi/2); tau = -c*xc + sqrt(1-(s*xc)^2);
thetaSummer = MyAtan(tau*s,xc+tau*c);
c = cos(phi+pi); s = sin(phi+pi); tau = -c*xc + sqrt(1-(s*xc)^2);
thetaAutumn = MyAtan(tau*s,xc+tau*c);
c = cos(phi+3*pi/2); s = sin(phi+3*pi/2); tau = -c*xc + sqrt(1-(s*xc)^2);
thetaWinter = MyAtan(tau*s,xc+tau*c);
% Compute the season lengths based on the season angles..
delSpring = thetaSummer-thetaSpring;
if delSpring>=0
    SpringLength = delSpring*aveYear/(2*pi);
else
    SpringLength = (delSpring+2*pi)*aveYear/(2*pi);
end
delSummer = thetaAutumn-thetaSummer;
if delSummer>=0
    SummerLength = delSummer*aveYear/(2*pi);
else
    SummerLength = (delSummer+2*pi)*aveYear/(2*pi);
end
delAutumn = thetaWinter-thetaAutumn;
if delAutumn>=0
    AutumnLength = delAutumn*aveYear/(2*pi);
else
    AutumnLength = (delAutumn+2*pi)*aveYear/(2*pi);
end
delWinter = thetaSpring-thetaWinter;
if delWinter>=0
    WinterLength = delWinter*aveYear/(2*pi);
else
    WinterLength = (delWinter+2*pi)*aveYear/(2*pi);
end
% Vector of start times as predicted by the eccentric model...
k = length(S)/4;
R = zeros(4*k,1);
R(1) = SpringLength;
for i=2:4*k
    if mod(i,4)==2
        R(i) = R(i-1)+SummerLength;
    elseif mod(i,4)==3
        R(i) = R(i-1)+AutumnLength;
    elseif mod(i,4)==0
        R(i) = R(i-1)+WinterLength;
    else
        R(i) = R(i-1)+SpringLength;
    end
end
% The vector of discrepancies...
R = R-S;

  function theta = MyAtan(y,x)
% Determines theta (0<=theta<2pi) so that (x,y)= (rcos(theta),rsin(theta))
% where r = sqrt(x^2+y^2).
theta = atan2(y,x);
if theta<0
    theta = 2*pi+theta;
end

% -2 for comparing season lengths instead of cumulative start times
% -1 for phi off by pi/2
