% 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];
v0 = input('Enter initial tangential velocity - '); % (au/yr)
v = [0 v0];
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;
for istep=1:nstep
rplot(istep) = norm(r); % Record orbit for polar plot
thplot(istep) = atan2(r(2),r(1));
tplot(istep) = time;
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,'+')
grid
ylabel('Distance (AU)')
title('Orbital motion')
subplot(122)
totalE = kinetic + potential;
plot(tplot,kinetic,'-.',tplot,potential,'--',tplot,totalE,'-')
xlabel('Time (yr)')
ylabel('Energy')
title('KE (Dot) PE (Dash) Total (Solid)')
subplot(111)
