How to give trail to particles/ How to fade out plots with time (complex example)

Inspired from a webproject I have created this n-body-simulation:

First of a function to acces on:

% calculating acceleration of particles 

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

followed by the Main script:

clear 

%% setting initial conditions

% Number of particles

N = 10; 

% 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 

figure
 
grid on
 
view(3)


%% 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;
    
    hold all
    
    for j = 1:N
        plot3(pos(j,1),pos(j,2),pos(j,3),'b.','MarkerSize',4)
    end

    
    pause(0.005)
end

However I'm running into a problem here: the particles path will not disapear gradually since all positions are hold. I saw some clever ways implementing a fade but I'm not able adapting it here. Maybe someone brave enough working through all the code may give some advises.

Answer 1

[DGM]

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.

[DGM]

This is doing single-point scatter inside the loop, adjusting face alpha for every step.

% setting initial conditions
% Number of particles
N = 3; 
% 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)

% each trajectory has its own solid color
cmap = hsv(N);
% alpha step size
dalph = 1/50; 

% starting main loop
hold all
for kt = 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;
    
    for kn = 1:N
        % can't change marker alpha for line objects afaik
        %H(kn,kt) = plot3(pos(kn,1),pos(kn,2),pos(kn,3),'.','MarkerSize',4,'color',cmap(kn,:));
        % it works for scatter objects, but it's obscenely slow
        H(kn,kt) = scatter3(pos(kn,1),pos(kn,2),pos(kn,3),5,cmap(kn,:),'filled');
        H(kn,kt).AlphaDataMapping = 'none';
        
        % redo alpha for older points
        for kkt = kt:-1:1
            H(kn,kkt).MarkerFaceAlpha = max(H(kn,kkt).MarkerFaceAlpha - dalph,0);
        end
    end

    % pause(0.001) % don't even need this because it's so slow anyway
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

Can it be done better? Probably. I can't run any of this due to version differences in scatter3(), and running it here makes the forum post editor ridiculously unstable and laggy. I'm bailing out at this point.

[Niklas Kurz]

Even if damn slow you really did manage doing it! Thankful Congratulations.

[DGM]

There may be some way of addressing the individual descendant (marker) objects in a scatter plot that I don't know about (perhaps something undocumented). If that were possible, then you might only need N scatter objects instead of N*t_end objects. I don't know if that is an option, but it might be a big boost to speed if it were.

Answer 2

[Image Analyst]

Perhaps if you used scatter() to create the markers. You can create an array with the colors for all the markers so have your older ones be faded and your newer ones be brighter/bolder/more colored. You can then put hold on and call plot() if you want lines between the markers.