% orbrk - Program to compute the orbit of a comet
% using the 4th order Runge-Kutta method
clear; help orbrk; % 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;
state = [ r(1) r(2) v(1) v(2) ];
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
state = rk4(state,time,tau,'gravrk',GM);
r = [state(1) state(2)];
v = [state(3) state(4)];
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)
