- #1
logix88
- 12
- 0
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
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
Code:
%n = no.of particles REPLACE 1000 with n
data=load('particledata.dat');
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: