orbe.m

% orbe - Program to compute the orbit of a comet
% using the Euler method.
clear;  help orbe;  % Clear memory and print header
r0 = input('Enter initial radial distance - ');  % (au)
r = [r0 0];       % Initial position is on x-axis
v0 = input('Enter initial tangential velocity - '); % (au/yr)
v = [0 v0];       % Initial velocity is in the y-direction
tau = input('Enter time step, tau - ');  % (yr)
GM = 4*pi^2;      % Grav. const. * Mass of Sun (au^3/yr^2)
mass = 1.;        % Mass of projectile 
%%%%% MAIN LOOP %%%%%%
time = 0;
nstep = 200;      % Number of time steps
for istep=1:nstep  
  rplot(istep) = norm(r);       % Record orbit for polar plot
  thplot(istep) = atan2(r(2),r(1));
  tplot(istep) = time;          % Record time for plot
  kinetic(istep) = .5*mass*norm(v)^2;  % Record energies
  potential(istep) = - GM*mass/norm(r);
  % Calculate new position and velocity
  accel = -GM*r/norm(r)^3;   % Gravity
  r = r + tau*v;       % Euler step
  v = v + tau*accel; 
  time = time + tau;    
end
% Graph the trajectory of the comet
subplot(121)
  polar(thplot,rplot,'+')  % Polar plot of trajectory
  grid                     % Include a grid on the plot
  ylabel('Distance (AU)')
  title('Orbital motion')
subplot(122)
  totalE = kinetic + potential;  % Total energy
  plot(tplot,kinetic,'-.',tplot,potential,'--',tplot,totalE,'-')
  xlabel('Time (yr)')
  ylabel('Energy')
  title('KE(dot); PE(dash); Total(solid)')
subplot(111)