% orbit.m: Demo of numerical modeling issues from orbital mechanics.
%
% Numerically integrate  x''(t) = -x(t)/||x||^3 in two dimensions
% x = position of a satellite about a planet, which is at position (0,0).
% The orbit should be a perfect circle; deviations are due to errors in 
% the method used to numerically integrate the orbit.
% Compare the errors observed for the Euler and the Stoermer/Verlet methods.

 fprime = inline('[x(3:4); -x(1:2)/norm(x(1:2))^3]','x');    % define the RHS of the ODE

 r0 = [1;0];                  % initial position of satellite
 v0 = [0;1];                  % initial velocity of satellite
 Tfinal = 25;                 % integrate from 0 to Tfinal

  
 method = input(' (1) Euler method or (2) Stoermer/Verlet? ');
 h  = input('time step, h (e.g, 0.1) ? ');

 t  = [0:h:Tfinal];           % set up time grid

 figure(1), clf
 set(gcf,'doublebuffer', 'on')                % prevents screen flashing during loop

% solve the differential equation using Euler's method or the Stoermer/Verlet method

 r  = r0;
 v  = v0;

 E = zeros(size(t));       % energy
 E(1) = .5*(v(:,1)'*v(:,1))-1/norm(r(:,1));
 
 for k=2:length(t)
    if method==1           % Euler
       r(:,k) = r(:,k-1) + h*v(:,k-1);
       v(:,k) = v(:,k-1) - h*r(:,k-1)/norm(r(:,k-1))^3;
    else                   % Stoermer/Verlet
       rhalf  =  r(:,k-1) + (h/2)*v(:,k-1);
       v(:,k) =  v(:,k-1) - h*rhalf/norm(rhalf)^3;
       r(:,k) =  rhalf    + (h/2)*v(:,k);
    end
    E(k) = .5*(v(:,k)'*v(:,k))-1/norm(r(:,k));

% plotting commands follow
    figure(1), clf   
    axes('position', [0 .33 1 .62])
    plot(r(1,:),r(2,:),'b-','linewidth',2), hold on
    plot(0,0,'r.','markersize',40)
    plot(r(1,end),r(2,end),'g.','markersize',30)
    axis equal, axis([-2.5 2.5 -2.5 2.5])
    xlabel('r_1', 'fontsize', 20) 
    ylabel('r_2', 'fontsize', 20) 
    text(3,2,sprintf('t = %5.3f', t(k)),'fontsize',18)

    axes('position', [.25 .05 .5 .15])
    plot(t(1:k), E(1:k), 'b-','linewidth',3)
    axis([0 t(end) -.55 -.2]) 
    text(26.5,-.3,'energy','fontsize',18)
    drawnow
 end