function ydot = astronomy(t, y)
%
%   ydot = astronomy(t, y)
%
% Implements the ODE for the motion of the third negligibly sized body in
% a three body gravitational problem.  The details are taken from Ascher
% & Greif, chapter 16, exercise 8.
%
% Parameters:
%
%   t    Current time.  Ignored.
%   y    Current state.  Stored as u_1, u_2, u_1', u_2'.
%
%   ydot Derivative of current state.

% Ian Mitchell for CPSC 303, 3/26/06.

  % Normalized masses.
  mu = 0.012277471;  
  mu_hat = 1 - mu;

  % Unpack input vector.
  % Positions.
  u1 = y(1);
  u2 = y(2);
  % Velocities.
  u1p = y(3);
  u2p = y(4);

  % More parameters.
  D1 = ((u1 + mu).^2 + u2.^2) .^ 1.5;
  D2 = ((u1 - mu_hat).^2 + u2.^2) .^ 1.5;
  
  % Output vector.
  ydot = zeros(4, 1);
    
  % Derivatives of position are velocities.
  ydot(1) = u1p;
  ydot(2) = u2p;
  
  % Derivatives of velocities are accelerations.
  ydot(3) = u1 + 2 * u2p - mu_hat * (u1 + mu) ./ D1 - mu * (u1 - mu_hat) ./ D2;
  ydot(4) = u2 - 2 * u1p - mu_hat * u2 ./ D1 - mu * u2 ./ D2;