function nbody
% NBODY  3-D motion of n bodies under mutual gravitional attraction.
% P = n-by-3 array of position coordinates
% V = n-by-3 array of velocities
% M = n-by-1 array of masses

   [P,V,M] = initialize_orbits;
   [H,stop,speed] = initialize_graphics(P);
   n = size(P,1);

   while ~get(stop,'value')
      % Advance velocities using current positions
      dt = get(speed,'value');
      for i = 1:n
         for j = i+1:n
            r = P(j,:) - P(i,:);
            g = dt*r/norm(r)^3;
            V(i,:) = V(i,:) + M(j)*g;
            V(j,:) = V(j,:) - M(i)*g;
         end
      end
     
      % Advance positions using advanced velocities
      P = P + dt*V;
      update_plot(H,P);
   end

   finalize_graphics(stop)
end 

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

function [P,V,M] = initialize_orbits
% The solar system.
% Our values for the inner planets are just reasonable guesses.

   sun.p = [0 0 0];
   sun.v = [0 0 0];
   sun.m = 1;

   me.p = [.3 0 .005];
   me.v = [0 1.75 0];
   me.m = 1.8e-7;

   ve.p = [.7 0 .01];
   ve.v = [0 1.1 0];
   ve.m = 2.4e-6;

   ea.p = [1 0 0];
   ea.v = [0 1 0];
   ea.m = 3e-6;

   ma.p = [1.5 0 -.02];
   ma.v = [0 .8 0];
   ma.m = 2e-6;

   ju.p = [4.8414 -1.1603 -0.10362];
   ju.v = [0.096500 0.44754 -0.0040136];
   ju.m = 9.5479e-04;
   
   sa.p = [8.3434 4.1248 -0.40352];
   sa.v = [-0.16087 0.29056 0.0013394]; 
   sa.m = 2.8588e-04;
   
   ur.p = [12.894 -15.111 -0.22331];
   ur.v = [0.17233 0.13826 -0.0017241];
   ur.m = 4.3662e-05;
   
   ne.p = [15.380 -25.919 0.17926];
   ne.v = [0.15583 0.094649 -0.0055316];
   ne.m = 5.1514e-05;
   
   P = [sun.p; me.p; ve.p; ea.p; ma.p; ju.p; sa.p; ur.p; ne.p];
   V = [sun.v; me.v; ve.v; ea.v; ma.v; ju.v; sa.v; ur.v; ne.v];
   M = [sun.m; me.m; ve.m; ea.m; ma.m; ju.m; sa.m; ur.m; ne.m];

end

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

function [H,stop,speed] = initialize_graphics(P)
%  color =       [ gold magenta  gray   blue   red   dkred  orange cyan dkgreen];
   color = (1/4)*[4 3 0; 2 0 2; 1 1 1; 0 0 3; 4 0 0; 2 0 0; 4 2 0; 0 3 3; 0 2 0];
   marksize = [36 12 18 18 16 30 24 20 18]';
   n = size(P,1);
   clf reset
   for i = 1:n
      H(i) = line(P(i,1),P(i,2),P(i,3),'color',color(i,:), ...
         'marker','.','markersize',marksize(i));
   end
   axis(32*[-1 1 -1 1 -1/8 1/8])
   box on
   cameratoolbar
   cb = ['axis(s*axis);' ...
      'b=findobj(''style'',''slider'');' ...
      'set(b,''max'',s*get(b,''max''));' ...
      'set(b,''value'',s*get(b,''value''))'];
   speed = uicontrol('style','slider','min',0,'max',1/2,'value',1/8, ...
      'units','normal','position',[.02 .02 .36 .03],'sliderstep',[1/80 1/40]);
   stop = uicontrol('string','stop','style','toggle', ...
     'units','normal','position',[.90 .02 .08 .03]);
   uicontrol('string','in','style','pushbutton','units','normal', ...
      'position',[.42 .02 .06 .03],'callback',['s=1/sqrt(2);' cb])
   uicontrol('string','out','style','pushbutton','units','normal', ...
      'position',[.50 .02 .06 .03],'callback',['s=sqrt(2);' cb])
   uicontrol('string','x','style','pushbutton','units','normal', ...
      'position',[.58 .02 .06 .03],'callback','view(0,0)')
   uicontrol('string','y','style','pushbutton','units','normal', ...
      'position',[.66 .02 .06 .03],'callback','view(90,0)')
   uicontrol('string','z','style','pushbutton','units','normal', ...
      'position',[.74 .02 .06 .03],'callback','view(2)')
   uicontrol('string','3d','style','pushbutton','units','normal', ...
      'position',[.82 .02 .06 .03],'callback','view(3)')
   drawnow
end

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

function update_plot(H,P)
   n = size(P,1);
   for i = 1:n
      set(H(i),'xdata',P(i,1),'ydata',P(i,2),'zdata',P(i,3))
   end
   drawnow
end

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

function finalize_graphics(stop)
   delete(findobj('style','pushbutton'))
   delete(findobj('style','slider'))
   set(stop,'string','close','value',0,'callback','close')
end