Here's one way using scatter3()
% setting initial conditions
% Number of particles
N = 3; % i'm using fewer points for ease of viewing
% max distance between partices at start
rmax = 4;
pos = rmax * rand(N,3);
% max velocity for particles
vmax = 1;
vel = vmax * rand(N,3);
% max mass of particles
mmax = 5;
mass = mmax * ones(N,1);
% minimal distance parameter
softening = 0.5;
% calculating acceleration of particles
acc = getAcc(pos,mass,softening);
% selecting start time
t_start = 0;
% selecting timestep
dt = 0.04;
% selecting end time
t_end = 250;
% initializing figure
grid on
view(3)
% preallocate outputs
xpos = zeros(N,t_end);
ypos = zeros(N,t_end);
zpos = zeros(N,t_end);
% starting main loop
for i = 1:t_end
vel = vel + acc * dt/2;
[pos] = pos + vel * dt;
acc = getAcc(pos, mass, softening);
vel = vel + acc * dt/2;
t_start = t_start + dt;
% store outputs
xpos(:,i) = pos(:,1);
ypos(:,i) = pos(:,2);
zpos(:,i) = pos(:,3);
end
% each trajectory has its own solid color
cmap = hsv(N);
% alpha is a linear gradient from 1 to 0
alph = fliplr(linspace(0,1,t_end));
for ks = 1:N
hs = scatter3(xpos(ks,:),ypos(ks,:),zpos(ks,:),10,cmap(ks,:),'filled');
hold on;
hs.AlphaData = alph;
hs.AlphaDataMapping = 'none';
hs.MarkerFaceAlpha = 'flat';
end
function [a] = getAcc(pos,mass,softening)
% gravitation constant
G = 1;
% Number of particles
N = size(pos,1);
% acceleration matrix
a = zeros(N,3);
% loop to fill in acceleration by interaction with each particle
for i = 1:N
for j = 1:N
if i ~= j
dx = pos(j,1)-pos(i,1);
dy = pos(j,2)-pos(i,2);
dz = pos(j,3)-pos(i,3);
dr = (dx.^2 + dy.^2 + dz.^2+ softening.^2).^(-3/2);
a(i,1) = a(i,1) + G * mass(j,1) * (dr * dx);
a(i,2) = a(i,2) + G * mass(j,1) * (dr * dy);
a(i,3) = a(i,3) + G * mass(j,1) * (dr * dz);
end
end
end
end
I imagine this isn't really what you want though. I imagine you want something that fades points dynamically as new points are plotted. That can be done, but it'll be slow.
