  function P3(yp)
% yp = dy/dt at t = 0. Assume 1<=yp<=1.4.
% We will use this global variable to count f-evals...
global fCount
% This will be used for absolute error control...
tol = .0000000001;
% Will integrate across the interval [t0,tFinal] understanding that the
% integration will stop (roughly) after one revolution. We set tFinal to a
% very large value to insure that we complete one orbit3
t0 = 0;
tFinal =  10000;
% The initial con
u0 = [ 1; 0 ; 0 ; yp];

% Determine the approximate period with a call to ode45 and an appropriate
% event function...
fCount = 0;
options = odeset('AbsTol',tol,'events',@g);
[t, u] = ode45(@Kepler, [0 tFinal], u0,options);
% Because of the event function used, the last t-value will approximately
% equal the period...
m = length(t);
Period = t(m);
% Plot the orbit and display the orbit points that are computed by ode45
plot(u(:,1),u(:,3),u(:,1),u(:,3),'.r')
axis equal
shg
% Report the total number of f-evaluations, the period, and the yp value
fCount = fCount;
title(sprintf('yp = %5.3f  Period = %6.4f    f-evals = %3d',yp,Period,fCount),'fontsize',14)

% Kepler Check
% Fit splines to the discrete snapshots of x(t), dx(t)/dt, y(t), and
% dy(t)/dt...
Sx  = spline(t,u(:,1));
Sxp = spline(t,u(:,2));
Sy  = spline(t,u(:,3));
Syp = spline(t,u(:,4));
% Compute 100 different sector areas
tVals = linspace(0,Period,101);
Avals = zeros(100,1);
for j=1:100
    % The area swept out by the radius vector from 
    % t = (j-1)P/100 to t = jP/100
    Avals(j) = integral(@(t) AreaIntegrand(t,Sx,Sxp,Sy,Syp),tVals(j),tVals(j+1),'Abstol',.0000001);
end
mu = mean(Avals);
KeplerCheck = max(abs(Avals-mu))/mu;
xlabel(sprintf('KeplerCheck = %7.3e',KeplerCheck),'fontsize',14)
shg


function z = AreaIntegrand(t,Sx,Sxp,Sy,Syp)
% The integrand function for computing the area swept out by the
% radius vector.
z = ppval(Sx,t).*ppval(Syp,t) - ppval(Sxp,t).*ppval(Sy,t);


function [gstop,isterminal,direction] = g(t,u)
% Event function for turning off ode45 after (approximately) one
% revolution. We do this when dx/dt is zero and decreasing.
gstop = u(2) ;
isterminal = 1;
direction = -1;

  function up = Kepler(t,u)
% up = Kepler(t,u)
% t (time) is a scalar and u is a 4-vector whose components satisfy
%
%            u(1) = x(t)      u(2) = (d/dt)x(t)
%            u(3) = y(t)      u(4) = (d/dt)y(t)
%
% where (x(t),y(t)) are the equations of motion in the 2-body problem.
%
% up is a 4-vector that is the derivative of u at time t.
global fCount
r3 = (u(1)^2 + u(3)^2)^1.5;
up = [  u(2)     ;...
       -u(1)/r3  ;...
        u(4)     ;...
       -u(3)/r3] ;
fCount = fCount+1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5 pts  Reasonable event function for getting the period 
% 3 pts setting up the integrand function for kepler Check
% 2 pts for display of orbit, evaluation points, etc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%