% 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]; % Initial position is on x-axis
v0 = input('Enter initial tangential velocity - '); % (au/yr)
v = [0 v0]; % Initial velocity is in 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; % Total number of time steps
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; % Record time for plot
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)]; % Unravel position and velocity
v = [state(3) state(4)]; % from state vector
time = time + tau;
fprintf('Finished %g steps out of %g \n',istep,nstep);
end
% Graph the trajectory of the comet
subplot(121)
polar(thplot,rplot,'+') % Polar plot of orbit
grid % Put 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)
