function orbits
% ORBITS  Motion of n bodies under mutual gravitional attraction.
%    ORBITS The solar system with n = 9 for one sun and 8 planets.
%    Graphics user interface similar to BROWNIAN and FIREWORKS.

   % P = n-by-d array of position coordinates, d = 2 or 3.
   % V = n-by-d array of velocities
   % M = n-by-1 array of masses
   % H = graphics and user interface handles

   [P,V,M] = initialize_orbits;
   H = initialize_graphics(P);

   n = size(P,1);  % Number of bodies
   steps = 20;     % Number of steps between plots
   t = 0;          % time

   while ~get(H.stop,'value')

      % Obtain step size from slider.
      delta = get(H.speed,'value')/steps;
      
      for k = 1:steps

         % Compute current gravitational forces.
         G = zeros(size(P));
         for i = 1:n
            for j = [1:i-1 i+1:n];
               r = P(j,:) - P(i,:);
               G(i,:) = G(i,:) + M(j)*r/norm(r)^3;
            end
         end
 
         % Update velocities using current gravitational forces.
         V = V + delta*G;
        
         % Update positions using updated velocities.
         P = P + delta*V;

      end

      t = t + steps*delta;
      update_plot(P,H,t)
   end
   finalize_graphics
end 

% ------------------------------------------------------------

function [P,V,M] = initialize_orbits
% The solar system.
% Obtain data from Jet Propulsion Laboratory HORIZONS.
% http://ssd.jpl.nasa.gov/horizons.cgi  
% Ephemeris Type: VECTORS
% Coordinate Orgin: Sun (body center)
% Time Span: 2008-7-24 to 2008-7-25

   sol.p = [0 0 0];
   sol.v = [0 0 0];
   sol.m = 1.9891e+30;

   mer.p = [-1.02050180e-2  3.07938393e-1  2.60947941e-2];
   mer.v = [-3.37623365e-2  9.23226497e-5  3.10568978e-3];
   mer.m = 3.302e+23;
   
   ven.p = [-6.29244070e-1  3.44860019e-1  4.10363705e-2];
   ven.v = [-9.80593982e-3 -1.78349270e-2  3.21808697e-4];
   ven.m = 4.8685e+24;
   
   ear.p = [ 5.28609710e-1 -8.67456608e-1  1.28811732e-5];
   ear.v = [ 1.44124476e-2  8.88154404e-3 -6.00575229e-7];
   ear.m = 5.9736e+24;
   
   mar.p = [-1.62489742e+0 -2.24489575e-1  3.52032835e-2];
   mar.v = [ 2.43693131e-3 -1.26669231e-2 -3.25240784e-4];
   mar.m = 6.4185e+23;
   
   jup.p = [ 1.64800250e+0 -4.90287752e+0 -1.65248109e-2];
   jup.v = [ 7.06576969e-3  2.76492888e-3 -1.69566833e-4];
   jup.m = 1.8986e+27;
   
   sat.p = [-8.77327303e+0  3.13579422e+0  2.94573194e-1];
   sat.v = [-2.17081741e-3 -5.26328586e-3  1.77789483e-4];
   sat.m = 5.6846e+26;
   
   ura.p = [ 1.97907257e+1 -3.48999512e+0 -2.69289277e-1];
   ura.v = [ 6.59740515e-4  3.69157117e-3  5.11221503e-6];
   ura.m = 8.6832e+25;
   
   nep.p = [ 2.38591173e+1 -1.82478542e+1 -1.74095745e-1];
   nep.v = [ 1.89195404e-3  2.51313400e-3 -9.54022068e-5];
   nep.m = 1.0243e+26;

   P = [sol.p; mer.p; ven.p; ear.p; mar.p; jup.p; sat.p; ura.p; nep.p];
   V = [sol.v; mer.v; ven.v; ear.v; mar.v; jup.v; sat.v; ura.v; nep.v];
   M = [sol.m; mer.m; ven.m; ear.m; mar.m; jup.m; sat.m; ura.m; nep.m];

   % Scale mass by solar mass.
   M = M/sol.m;

   % Scale velocity to radians per year.
   V = V*365.25/(2*pi);

   % Adjust sun's initial velocity so that system total momentum is zero.
   V(1,:) = -sum(diag(M)*V);

end

% ------------------------------------------------------------

function  H = initialize_graphics(P)
% Initialize graphics and user interface controls

   dotsize = [36 12 18 18 16 30 24 20 18]';
   color = [4 3 0     % gold
            2 0 2     % magenta
            1 1 1     % gray
            0 0 3     % blue
            4 0 0     % red
            3 0 0     % dark red
            4 2 0     % orange
            0 3 3     % cyan
            0 2 0]/4; % dark green

   clf reset
   set(gcf,'numbertitle','off','name',' n-body')
   s = max(sqrt(diag(P*P')));
   axis([-s s -s s -s/4 s/4])
   axis square
   box on
   n = size(P,1);
   for i = 1:n
      H.planets(i) = line(P(i,1),P(i,2),P(i,3),'color',color(i,:), ...
         'marker','.','markersize',dotsize(i),'userdata',dotsize(i));
   end

   H.clock = title('0 years','fontweight','normal');
   H.speed = uicontrol('style','slider','min',0,'value',1/2,'max',1, ...
      'units','normal','position',[.02 .02 .30 .04],'sliderstep',[1/80 1/40]);
   H.stop = uicontrol('string','stop','style','toggle', ...
     'units','normal','position',[.90 .02 .08 .04]);
   set(gcf,'userdata',H)

   uicontrol('string','trace','style','toggle','units','normal', ...
      'position',[.34 .02 .06 .04],'callback','tracer');
   uicontrol('string','in','style','pushbutton','units','normal', ...
      'position',[.42 .02 .06 .04],'callback','zoomer(1/sqrt(2))')
   uicontrol('string','out','style','pushbutton','units','normal', ...
      'position',[.50 .02 .06 .04],'callback','zoomer(sqrt(2))')
   uicontrol('string','x','style','pushbutton','units','normal', ...
      'position',[.58 .02 .06 .04],'callback','view(0,0)')
   uicontrol('string','y','style','pushbutton','units','normal', ...
      'position',[.66 .02 .06 .04],'callback','view(90,0)')
   uicontrol('string','z','style','pushbutton','units','normal', ...
      'position',[.74 .02 .06 .04],'callback','view(0,90)')
   uicontrol('string','3d','style','pushbutton','units','normal', ...
      'position',[.82 .02 .06 .04],'callback','view(-37.5,30)')
   drawnow
end

% ------------------------------------------------------------

function tracer
% Callback for trace button
   H = get(gcf,'userdata');
   planets = flipud(H.planets);
   trace = get(gcbo,'value');
   n = length(planets);
   for i = 1:n
      if trace
         ms = 6 + 18*(i==1);
         set(planets(i),'markersize',ms,'erasemode','none')
      else
         ms = get(planets(i),'userdata');
         set(planets(i),'markersize',ms,'erasemode','normal')
      end
   end
   if trace
      set(H.clock,'erasemode','xor')
   else
      set(H.clock,'erasemode','normal')
   end
end

% ------------------------------------------------------------

function zoomer(zoom)
% Callback for in and out buttons
   H = get(gcf,'userdata');
   T = view;
   view(3);
   axis(zoom*axis);
   view(T);
   set(H.speed,'max',zoom*get(H.speed,'max'), ...
      'value',zoom*get(H.speed,'value'));
end

% ------------------------------------------------------------

function update_plot(P,H,t)
   set(H.clock,'string',sprintf('%10.2f years',t/(2*pi)))
   for i = 1:size(P,1)
      set(H.planets(i),'xdata',P(i,1),'ydata',P(i,2),'zdata',P(i,3))
   end
   drawnow
end

% ------------------------------------------------------------

function finalize_graphics
   delete(findobj('type','uicontrol'))
   uicontrol('string','close','style','pushbutton', ...
     'units','normal','position',[.90 .02 .08 .04],'callback','close');
end