orbe_p.m

clear;  %%% Plot data from orbe.f, orbec.f, orbrk.f or orbrka.f %%%
load tplot.dat; load thplot.dat; load rplot.dat; 
load kinetic.dat;  load potential.dat;
% 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)