% Demonstrates use of nbodyplot, which plots particle trajectories,
% assuming motion is restricted to the xy plane.
% 
% Data represents particles moving on circular Keplerian
% orbits at various distances from a central mass (i.e. a mock 
% "solar system").

% Number of timesteps ...
nt = 1501;

% Minimum and maximum times ...
tmin = 0;
tmax = 8.0;

% Vector of (output) times ...
t = linspace(tmin, tmax, nt);

% Number of particles ...
np = 9;

% Minimum and maximum particle masses ...
mmin = 1.5;
mmax = 3.0;

% Vector of particle masses ...
m = linspace(mmin, mmax, np);

% Particle color array ... use primary colors cyclically ...
rgb = zeros(np,3);
for ip = 1 : np
   if rem(ip,7) == 1
      rgb(ip,:) = [1 0 0];
   elseif rem(ip,7) == 2
      rgb(ip,:) = [0 1 0];
   elseif rem(ip,7) == 3
      rgb(ip,:) = [0 0 1];
   elseif rem(ip,7) == 4
      rgb(ip,:) = [1 1 0];
   elseif rem(ip,7) == 5
      rgb(ip,:) = [1 0 1];
   elseif rem(ip,7) == 6
      rgb(ip,:) = [0 1 1];
   else
      rgb(ip,:) = [1 1 1];
   end
end

% Minimum and maximum distances for particle orbits ...
dmin = 0.5;
dmax = 1.5;

% Vector of orbital distances ...
d = linspace(dmin, dmax, np);

% Initial angular positions of particles ...
theta0 = 2 * pi * rand(1,np);

% Angular frequencies of particles, corresponding to
% Keplerian orbits ...
omega = d.^(-3/2);

% Create and initialize coordinate array ...
r = zeros(np, 3, nt);

% Generate coordinate data---all particles move in circular
% orbits in the xy plane ...
for it = 1 : nt
   for ip = 1 : np
      r(ip, 1, it) = d(ip) * cos(theta0(ip) + omega(ip) * t(it));
      r(ip, 2, it) = d(ip) * sin(theta0(ip) + omega(ip) * t(it));
      r(ip, 3, it) = 0;
   end
end

%{
Note: This is the beginning of a multiline comment ...

The 'for ip = 1 : np' loop above could be eliminated by
using 

   1) The ':' notation to select all elements along a specific
      dimension of an array and 
   2) The element-wise multiplication operator

for it = 1 : nt
      r(:, 1, it) = d(:) .* cos(theta0(:) + omega(:) .* t(it));
      r(:, 2, it) = d(:) .* sin(theta0(:) + omega(:) .* t(it));
      r(:, 3, it) = 0;
end

End of multiline comment
%}

% Plot the trajectories ...
nbodyplot(t, r);
% Most general form of nbodyplot invocation ...
% nbodyplot(t, r, plotoptions, markersize)
