# Particle Simulation in 3D, help please!

1. Dec 3, 2009

### logix88

I am trying to simulate a 3D particle simulation, where there are a 1000 particles of random x,y,z postn varying area between 2 n 4, and mass... There are different forces acting on each particles, Force due to heat, wind, and gravity, we calc acceleration from that, after which change in distance and hence new x,y,z positions are obtained. I have written most of it (at least I thnk I did),

I need help updating these new x,y,z postion at time range (1:10), I am sort of stuck there.

I also need to lock down this 3D space to a 100x100m box, it seems to go beserk,...

and finally, these particles should get stuck to the top and bottom if they reach there... and if it goes beyond left or right, it should appear in the opposite side.

Here's the code, let me knw what you guys think.

NB : particle data file here : http://bit.ly/6qJaNb [Broken]

Code (Text):

%n = no.of particles  REPLACE 1000 with n

scatter3(data(:,1),data(:,2),data(:,3),data(:,5),'filled'), view(-60,60)

m=data(:,4);

g=9.81;

x=data(:,1);
y=data(:,2);
z=data(:,3);

F_w=zeros(3,1,1000);

F_g=zeros(3,1,1000);

F_h=zeros(3,1,1000);

q=0.08;

A=data(:,5);

for t=2:1:10
for w=1:1:1000;

F_g_z(w)=-m(w)*g*z(w);   % force due to grav

F_g(1,1,w)=0;
F_g(2,1,w)=0;
F_g(3,1,w)=F_g_z(w);

rx=50-y(w);
ry=x(w)-50;

F_w_x(w)=1.6*(exp(-(sqrt(rx^2+ry^2)-25/20)^2))*rx;    % force wind
F_w_y(w)=1.6*(exp(-(sqrt(rx^2+ry^2)-25/20)^2))*ry;
F_w(1,1,w)=F_w_x(w);
F_w(2,1,w)=F_w_y(w);
F_w(3,1,w)=0;

T(w)=1/16*(70*x(w)+(50-y(w))^2);

F_h1(w)=q*A(w)*T(w)*z(w);             %force heat

F_h(1,1,w)=0;
F_h(2,1,w)=0;
F_h(3,1,w)=F_h1(w);

F=F_g+F_w+F_h;

a=F/m(w);                     % acceleration

u(1,1,w,1)=a(1,1,w);
u(2,1,w,1)=a(2,1,w);
u(3,1,w,1)=a(3,1,w);

u(1,1,w,t)=u(1,1,w,t-1)+a(1,1,t)*t;    % to store u in x direction, of w-th particle with time t
u(2,1,w,t)=u(2,1,w,t-1)+a(2,1,t)*t;
u(3,1,w,t)=u(3,1,w,t-1)+a(3,1,t)*t;

%             end

%             fprintf('%u\n',u);

s(1,1,w,t)=u(1,1,w,t)*t+(((1/2)*a(1,1,w))*t^2);   % s=ut+1/2at^2;
s(2,1,w,t)=u(2,1,w,t)*t+(((1/2)*a(2,1,w))*t^2);
s(3,1,w,t)=u(3,1,w,t)*t+(((1/2)*a(3,1,w))*t^2);

%             ds=s_fin-s_ini;
ds(1,1,w,t)=s(1,1,w,t)-s(1,1,w,t-1);
ds(2,1,w,t)=s(2,1,w,t)-s(2,1,w,t-1);
ds(3,1,w,t)=s(3,1,w,t)-s(3,1,w,t-1);

x_change(w,t)=ds(1,1,w,t);
y_change(w,t)=ds(2,1,w,t);
z_change(w,t)=ds(3,1,w,t);

x_final(w,1)=0;
y_final(w,1)=0;
z_final(w,1)=0;

x_final(w,t)=x_final(w,t-1)+x_change(w,t);

y_final(w,t)=y_final(w,t-1)+y_change(w,t);

z_final(w,t)=z_final(w,t-1)+z_change(w,t);

x(w)=x_final(w,t);
y(w)=y_final(w,t);
z(w)=z_final(w,t);

%             fprintf('%z\n',z(w));

% %
%             scatter3(x_final(w,t),y_final(w,t),z_final(w,t),'filled'), view(-60,60)
%             pause(0.1);

%             NEED TO UPDATE particledata / 3D to change accordingly
%             HOW TO?
%             ballp(i)=(x(i),y(i),z(i))
%             dt = 0.5
%             t=0.0
%
%
%             while True:
%
%               t = t + dt
%               ballpos(i) = ballpos(i) + (ballp(i)/m(i))*dt
%
%               if not (side > ball.x > -side):   (0,Y,Z) OR (100,Y,Z)
%                 ball.p.x = -ball.p.x
%               if not (side > ball.y > -side):   (X,0,Z) OR (X,100,Z)
%                 ball.p.y = -ball.p.y
%               if not (side > ball.z > -side):
%                 ball.p.z = -ball.p.z   HERE BALL IS STUCK TO TOP N BOTTOM  for any (X,Y,0) or (X,Y,100)

end

end

%
% for w=1:1000
%     for t=2:1:10
%
%     scatter3(x_final(w,t),y_final(w,t),z_final(w,t),'filled'), view(-60,60)
%     pause(0.1);
%     end
% end
%
%

Last edited by a moderator: May 4, 2017
2. Dec 3, 2009

### ApexOfDE

I dont understant what you means :( sry. When I glance at your code, I see that you can find x/y/z_final. You just need to draw a figure again to obtain new pos.

do you mean you want to see fig in 2d? If that, view(2) can help you.

It calls periodic condition. You can do this simply. Lets assume that z-pos ranges from 0 to 100 and there is another matrix to indicate whether particles have its z-pos out of this range (0 - no, 1 - yes). If (1), skip calculating this particle.

For left-right side, you can put a condition after calculating and modifying data for this special move. It should be more complicated.

p/s: I can only say this because I dont have time to read your code and completely understand it. Ill try it again when I am free :D

3. Dec 3, 2009

### logix88

Let me elaborate... I have x,y,z position of a 1000 particles. Now, over time these particles will be moving, due to the different forces applied on them, I want to simulate how these particles move around this 100x100x100m box. The t(1:10) was just time range going form 1 to 10, in steps of 1.

the 100mx100m was a mistake, it is supposed to be 100x100x100m

I have updated the code, to the one previously used, to make things a little clear.

Thanx!

4. Dec 4, 2009

### ApexOfDE

It is a simple code that i tried to write. In this code, I used your coordinate data and created new position by using random integers. If z-pos is out of [0 100] range, this particle will not move next time.

Code (Text):

clc
clear all

%
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
a_z = zeros(length(z), 1);

scatter3(x, y, z, data(:,5),'filled'), view(-60,60)

for t = 1 : 100
for j = 1 : length(x)
if a_z(j) == 1
continue
end

d = randi(10, 1, 3);
s = randi(2, 1, 3);

x(j) = x(j) + d(1) * (-1)^s(1);
y(j) = y(j) + d(2) * (-1)^s(2);
z(j) = z(j) + d(3) * (-1)^s(3);

if x(j) <= 0 || x(j) >= 100
C = (abs(x(j)) - x(j)) / (2 * abs(x(j)));
x(j) = 100 * C;
end

if y(j) <= 0 || y(j) >= 100
C = (abs(y(j)) - y(j)) / (2 * abs(y(j)));
y(j) = 100 * C;
end

%if p is at bottom or top, it cant move anymore.
if z(j) < 0 || z(j) > 100
C = (abs(z(j)) - z(j)) / (2 * abs(z(j)));
z(j) = 100 * C;
a_z(j) = 1;
end
end

scatter3(x, y, z, data(:,5),'filled'), view(-60,60)
pause(1)
end

Hope this helps.

5. Dec 4, 2009

### logix88

Thnx for the code, Its seems pretty straight forward with random x,y,z positions, its just when trying to implement forces and mapping them wrt to time, it seemed a bit too confusing... I will give it a ago, and come bak if I have any problems, Thanx a lot for ur help!

6. Dec 4, 2009

### D H

Staff Emeritus
logix88,

Your code has a *big* problem in the way you are propagating the particles from one time step to the next. If you want realistic physics, those particles need to have position and velocity. The acceleration changes the velocity, and the velocity in turn changes the position.

Make sure you do things in that order. If you update the position first and velocity second you will have what is called the Euler integration technique. Simply switch the order (velocity, then position) and you will have what is called symplectic Euler integration technique. The difference between the two is immense. Basic Euler integration is truly abysmal. Symplectic Euler, while not that great, represents a vast improvement in realism.

7. Dec 4, 2009

### logix88

I thought thats what I did, first I calculate Forces due to heat, wind, gravity for each particle, and I get acceleration from those values. Initial vel is 0, I then use the acceleration to obtain u=u0+at and then plug that into s=ut+1/2*at^2. Then obtain the change in distance for each x,y,z between previous time step and current time step and finally update the particles position. Do you mean, that I should update final position before distance? What other way can I do that?

Code (Text):

clc
clear all

%
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
m = data(:, 4);
a_z = zeros(length(z), 1);
g = 9.81;
q = 0.08;

A=data(:,5);

F_w=zeros(3,1,1000);

F_g=zeros(3,1,1000);

F_h=zeros(3,1,1000);

u=zeros(3,1,1000,100);

scatter3(x, y, z, data(:,5),'filled'), view(-60,60)

for t = 2 : 100
for j = 1 : length(x)
%         if a_z(j) == 1
%             continue
%         end

F_g_z(j)=-m(j)*g*z(j);

F_g(1,1,j)=0;
F_g(2,1,j)=0;
F_g(3,1,j)=F_g_z(j);

rx=50-y(j);
ry=x(j)-50;

F_w_x(j)=1.6*(exp(-(((sqrt(rx^2+ry^2))-25)/20))^2)*rx;
F_w_y(j)=1.6*(exp(-(((sqrt(rx^2+ry^2))-25)/20))^2)*ry;
F_w(1,1,j)=F_w_x(j);
F_w(2,1,j)=F_w_y(j);
F_w(3,1,j)=0;

T(j)=1/16*(70*x(j)+(50-y(j))^2);

F_h1(j)=q*A(j)*T(j)*z(j);

F_h(1,1,j)=0;
F_h(2,1,j)=0;
F_h(3,1,j)=F_h1(j);

F=F_g+F_w+F_h;
%
%     a=a_g+a_w+a_h;
a=F/m(j);

u(1,1,j,1)=a(1,1,j);
u(2,1,j,1)=a(2,1,j);
u(3,1,j,1)=a(3,1,j);

u(1,1,j,t)=u(1,1,j,t-1)+a(1,1,j)*t;
u(2,1,j,t)=u(2,1,j,t-1)+a(2,1,j)*t;
u(3,1,j,t)=u(3,1,j,t-1)+a(3,1,j)*t;

%             end

%             fprintf('%u\n',u);

s(1,1,j,t)=u(1,1,j,t)*t+(((1/2)*a(1,1,j))*t^2);
s(2,1,j,t)=u(2,1,j,t)*t+(((1/2)*a(2,1,j))*t^2);
s(3,1,j,t)=u(3,1,j,t)*t+(((1/2)*a(3,1,j))*t^2);

%             ds=s_fin-s_ini;
ds(1,1,j,t)=s(1,1,j,t)-s(1,1,j,t-1);
ds(2,1,j,t)=s(2,1,j,t)-s(2,1,j,t-1);
ds(3,1,j,t)=s(3,1,j,t)-s(3,1,j,t-1);

x_change(j,t)=ds(1,1,j,t);
y_change(j,t)=ds(2,1,j,t);
z_change(j,t)=ds(3,1,j,t);

x_final(j,1)=0;
y_final(j,1)=0;
z_final(j,1)=0;

x_final(j,t)=x_final(j,t-1)+x_change(j,t);

y_final(j,t)=y_final(j,t-1)+y_change(j,t);

z_final(j,t)=z_final(j,t-1)+z_change(j,t);

x(j)=x_final(j,t);
y(j)=y_final(j,t);
z(j)=z_final(j,t);

if x(j) <= 0 || x(j) >= 100
C = (abs(x(j)) - x(j)) / (2 * abs(x(j)));
x(j) = 100 * C;
end

if y(j) <= 0 || y(j) >= 100
C = (abs(y(j)) - y(j)) / (2 * abs(y(j)));
y(j) = 100 * C;
end

%if p is at bottom or top, it cant move anymore.
if z(j) < 0 || z(j) > 100
C = (abs(z(j)) - z(j)) / (2 * abs(z(j)));
z(j) = 100 * C;
a_z(j) = 1;
end
end

scatter3(x, y, z, data(:,5),'filled'), view(-60,60)
pause(0.1)
end

8. Dec 4, 2009

### D H

Staff Emeritus
This assumes constant acceleration. You don't have constant acceleration.

What you should be doing is something like this:

\aligned \vec u(t+\Delta t) &= \vec u(t) + \vec a(t)\Delta t \\ \vec x(t+\Delta t) &= \vec x(t) + \vec u(t+\Delta t) \Delta t \endaligned

You going from t=2 to 100 (and why 2?) in steps of 1. So in your case, $\Delta t = 1$.

Do not do this:

\aligned \vec u(t+\Delta t) &= \vec u(t) + \vec a(t)\Delta t \\ \vec x(t+\Delta t) &= \vec x(t) + \vec u(t)\Delta t \endaligned

That second one is the basic Euler method. The first one is symplectic Euler.

9. Dec 4, 2009

### logix88

u(1,1,j,1)=a(1,1,j);
u(2,1,j,1)=a(2,1,j);
u(3,1,j,1)=a(3,1,j);

u(1,1,j,t+1)=u(1,1,j,t)+a(1,1,j)*t;
u(2,1,j,t+1)=u(2,1,j,t)+a(2,1,j)*t;
u(3,1,j,t+1)=u(3,1,j,t)+a(3,1,j)*t;

x_final(j,1)=x(j);
y_final(j,1)=y(j);
z_final(j,1)=z(j);

x_final(j,t+1)=x_final(j,t)+u(1,1,j,t+1);
y_final(j,t+1)=y_final(j,t)+u(2,1,j,t+1);
z_final(j,t+1)=z_final(j,t)+u(3,1,j,t+1);

x(j)=x_final(j,t+1);
y(j)=y_final(j,t+1);
z(j)=z_final(j,t+1);

Is that what it should be like?

Last edited: Dec 4, 2009
10. Dec 4, 2009

### D H

Staff Emeritus
No. Why are you doing this?
Think of it this way: Suppose you had time going from 102 to 200? None of your forces depend on time. When time is zero is just a label. Moving that label doesn't change the physics.

11. Dec 4, 2009

### logix88

Hey, I am sry I dont understand what you said, I am not a physicist, I knw I have some prb understanding the simple physics involved in this.

TBH, I am more confused now than before...