function P1()
% Checks out Next Eclipse
% Compute the first 100 eclipses and report on your code's accuracy

clc
rand('seed',0)
disp('Eclipse   Start Time          Duration      Start Error    Duration Error')
disp('--------------------------------------------------------------------------')
t0 = 5;
for k=1:100
  [tNextCVL,dCVL] = NextEclipseCVL(t0);
  [tNext,d] = NextEclipse(t0);
  tNextError = abs(tNextCVL-tNext);
  dError = abs(dCVL-d);
  fprintf('%3d   %15.8f   %15.8f   %15.8f  %15.8f\n',k,tNextCVL,dCVL,tNextError,dError)
  % Pick a somewhat random t0 for the next eclipse...
  t0 = tNextCVL + 2*dCVL + 25*rand(1);
end


function [tNext,d] = NextEclipseCVL(t0)
% tNext is the time (days) when the next eclipse starts
% d is its duration

% Sample fEnter at four times...
z0  = fEnter(t0);
z10 = fEnter(t0+10);
z20 = fEnter(t0+20);
z30 = fEnter(t0+30);
% Given the properties of fEnter, the eclipse start is in one
% of the intervals [t0,t0+10], [t0+10,t0+20], [t0+20,t0+30].
% Sign values athe endpoints tell you which...
if z0>=0 && z10<=0
    tNext = fzero(@fEnter,[t0,t0+10]);
elseif z10>=0 && z20<=0
    tNext = fzero(@fEnter,[t0+10,t0+20]);
else
    tNext = fzero(@fEnter,[t0+20,t0+30]);
end
% Experimenting, the eclipse is never more than .3 days in duration...
tEnd = fzero(@fExit,[tNext,tNext+.3]);
d = tEnd-tNext;



function z = fEnter(t)
% Let RTP be the right tangent point as seen from the sun
% This function returns 
%              sine(angle Moon-Sun-Earth) - sine(angle RTP-Sun-Earth)
% This function is positive just before the eclipse starts and goes negative
% after the moon crosses the right shadow line
[xE,yE] = Earth(t); u = [xE;yE]/sqrt(xE^2+yE^2);
[xM,yM] = Moon(t);  v = [xM;yM]/sqrt(xM^2+yM^2);
r = .006371;
sigma = r/sqrt(xE^2+yE^2);
z = v(1)*u(2)-v(2)*u(1) - sigma;

function z = fExit(t)
% Let LTP be the left tangent point as seen from the sun
% This function returns 
%              sine(angle Earth-Sun-Moon) - sine(angle LTP-Sun-Earth)
% This function is negative just before the eclipse ends and goes positive
% after the moon crosses the left shadow line
[xE,yE] = Earth(t); u = [xE;yE]/sqrt(xE^2+yE^2);
[xM,yM] = Moon(t);  v = [xM;yM]/sqrt(xM^2+yM^2);
r = .006371;
sigma = r/sqrt(xE^2+yE^2);
z = u(1)*v(2)-u(2)*v(1) - sigma;

function [xE,yE] = Earth(t)
% Earth location at time t
pE = 147.098;
aE = 152.098;
rhoE = 365.256;
tau = 2*pi*t/rhoE;
xE = (pE-aE)/2 + ((pE+aE)/2)*cos(tau);
yE = sqrt(pE*aE)*sin(tau);

function [xM,yM] = Moon(t)
% Moon location at time t
pM = .363104;
aM = .405696;
rhoM = 27.321;
tau = 2*pi*t/rhoM;
[xE,yE] = Earth(t);
phi = atan2(yE,xE);
c = cos(phi);
s = sin(phi);
xM0 = (pM-aM)/2 + ((pM+aM)/2)*cos(tau);
yM0 = sqrt(pM*aM)*sin(tau);
xM = xE + c*xM0 - s*yM0;
yM = yE + s*xM0 + c*yM0;


% 7 points.
% 5 for all but a handful of the examples
% 2-4 if major starting interval problems/inaccuracy.

