function astronomy
% Script file to demonstrate 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.
%
% Ian Mitchell for CPSC 303, 3/26/06.
  
  % Problem parameters.
  y0 = [ 0.994; 0.0; 0.0; -2.00158510637908252240537862224 ];
  
  error_tolerance = 1e-6;
  tspan = [ 0, 17.1 ];
  
  % First, run with Matlab's ode45 to demonstrate the "accurate" trajectory.
  % Note that we must modify the default error tolerance threshold, and
  % that we also turn off output refinement to get a more accurate
  % measure of the actual steps.
  options = odeset('RelTol', error_tolerance, ...
                   'AbsTol', 0.01 * error_tolerance, 'Refine', 1);
  tic;
  [ ts, ys ] = ode45(@astroODE, tspan, y0, options);
  time = toc;

  % Plot the result.
  fprintf([ '\nMatlab ode45 in %d steps ' ...
            '(smallest: %5.3e, biggest: %5.3e) in %f seconds' ], ...
          length(ts), min(diff(ts)), max(diff(ts)), time);
  h = plot(ys(:,1), ys(:,2), 'b+-');
  set(h, 'LineWidth', 2);
  hold on;
  xlabel('u_1');  ylabel('u_2');

  legend_handles = h(1);
  legend_labels = { 'ode45, error tol 1e-6' };
  legend(legend_handles, legend_labels);
  
  % Now try running with classic fixed timestep RK4.
  steps =  [ 100; 200; 1000; 2000; 10000; 20000 ];
  hs = 17.1 ./ steps;
  styles = { 'r:', 'y:', 'm.-', 'c.-', 'k--', 'g--' };
  for i = 1:length(hs);
    % Wait to digest previous result.
    pause;

    % Get the next result.
    fprintf('\nRunning fixed timestep RK4 with %d steps (h = %5.3e)', ...
            steps(i), hs(i));
    tic
    [ ts, ys ] = rk4(@astroODE, tspan, y0, hs(i));
    time = toc;
    h = plot(ys(:,1), ys(:,2), styles{i});
    set(h, 'LineWidth', 2);

    fprintf(' in %f seconds', time);
    
    % Set the axis to see the latest results.
    axis([ min(ys(:,1)), max(ys(:,1)), min(ys(:,2)), max(ys(:,2)) ]);
    
    % Replace the legend with a bigger legend.
    new_label = sprintf('rk4, %d steps', steps(i));
    legend_handles = [ legend_handles, h(1) ];
    legend_labels = { legend_labels{:}, new_label };
    legend(legend_handles, legend_labels);
  
  end

  fprintf('\n');