function orbits3
% ORBITS3  Motion of n bodies under mutual gravitional attraction.
%          Artificial 3 body problem used in the text.

   % 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 = 100;    % 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 gravitational forces.
         V = V + delta*G;
        
         % Update positions using new 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
   % Initial position, velocity, and mass for three bodies.
   P = [ 0  0  0
         1  0  0
         0  1  0]; 
   V = [-0.12 -0.36  0
          0    0.72  0
         0.36 -0.36  0];
   M = [1/2 1/3 1/6]';
end

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

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

   dotsize = [18 14 10]';
   color = [4 3 0     % gold
            0 0 3     % blue
            4 0 0]/4; % red

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

   H.clock = title('0','fontweight','normal');
   H.speed = uicontrol('style','slider','min',0,'value',1/50,'max',1/10, ...
      'units','normal','position',[.02 .02 .30 .04],'sliderstep',[1/100 1/50]);
   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 = 4;
         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',t))
   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