function IdealGasSimulation
%IDEALGASSIMULATION Physics engine for multi-sphere collision.
%
%% PARAMETER
%Constants
G_newton = 0;%6.6738480e-11; %Newtonian constant of gravitation
g_grav = 0;%9.80665; %Gravitational acceleration
k_coulomb = 0;%8.987551787e9; %Coulomb's constant
Cr = 1; %Coefficient of restitution
dflt_axis_lim = [-10,10;-10,10;-10,10];
balls = repmat([0.2 1 0],10,1);
N = size(balls,1);
ball_pos = zeros(N,3);
rmax = max(balls(:,1))+0.5;
lx = dflt_axis_lim(1,1)+1+rmax:2*rmax:dflt_axis_lim(1,2)-1-rmax;
ly = dflt_axis_lim(2,1)+1+rmax:2*rmax:dflt_axis_lim(2,2)-1-rmax;
lz = dflt_axis_lim(3,1)+1+rmax:2*rmax:dflt_axis_lim(3,2)-1-rmax;
[mx,my,mz] = meshgrid(lx,ly,lz);
for kk = 1:N
ball_pos(kk,:) = [mx(kk),my(kk),mz(kk)];
end
ball_pos(:,3) = 0;
ball_v = 100*rand(N,3);
%ball_v(:,3) = 0;
%Properties
ball_color = repmat([0 0.5 1 0.75],N,1);
v = 0.5;
H = @heaviside;
w_pln_v = {...
{@(x) 0,@(x) 0,@(x) -v*(H(x)-1.5*H(x-20)+0.5*H(x-30)+H(x-30)-H(x-40))},...
{@(x) 0,@(x) 0,@(x) v*(H(x)-1.5*H(x-20)+0.5*H(x-30)+H(x-30)-H(x-40))},...
{@(x) v*(H(x)-1.5*H(x-20)+0.5*H(x-30)+H(x-30)-H(x-40)),@(x) 0,@(x) 0},...
{@(x) -v*(H(x)-1.5*H(x-20)+0.5*H(x-30)+H(x-30)-H(x-40)),@(x) 0,@(x) 0},...
{@(x) 0,@(x) -v*(H(x)-1.5*H(x-20)+0.5*H(x-30)+H(x-30)-H(x-40)),@(x) 0},...
{@(x) 0,@(x) v*(H(x)-1.5*H(x-20)+0.5*H(x-30)+H(x-30)-H(x-40)),@(x) 0}};
% w_sph_v = {...
% {@(x) 0,@(x) 0,@(x) 0}};
w_sph_v = {};
pln_pnt = cell(1,6); %{[pnt1;pnt2;pnt3;pnt4]}
pln_pnt{1} = [-10 -10 10;-10 10 10;10 10 10;10 -10 10];
pln_pnt{2} = [-10 -10 -10;-10 10 -10;10 10 -10;10 -10 -10];
pln_pnt{3} = [-10 -10 -10;-10 -10 10;-10 10 10;-10 10 -10];
pln_pnt{4} = [10 -10 -10;10 -10 10;10 10 10;10 10 -10];
pln_pnt{5} = [-10 10 -10;-10 10 10;10 10 10;10 10 -10];
pln_pnt{6} = [-10 -10 -10;-10 -10 10;10 -10 10;10 -10 -10];
pln_area = [100 100 100 100 100 100];
% wall_sph = [3,0,0,0];
wall_sph = [];
GravityVec = [0 0 -1];
GravitySrc = [];
CoulombSrc = [];
GroundPoint = [0 0 0];
friction_factor = 0;
axis_lim = [-10,10;-10,10;-10,10];
dt = 0.1;
tmax = 200;
debug_mode = 1;
%Variables
tlen = tmax/dt;
tlen_fig = 50/dt;
cldmax = 25000;
n_grav = size(GravitySrc,1);
n_coulomb = size(CoulombSrc,1);
n_pln = length(pln_pnt);
n_sph = size(wall_sph,1);
wall_pln = generate_walls;
impulse_frame = zeros(1,n_pln+n_sph); %[plane sphere]
if debug_mode
dn_press = 50; %measurement time interval for pressure
impulse_wall = zeros(dn_press,n_pln+n_sph); %[plane sphere]
press_factor = dn_press*dt*sum(pln_area,2); %not including spherical walls
n_v_dist = 10;
v_max = max(sum(sqrt(ball_v.^2),2));
v_int = linspace(0,v_max,n_v_dist);
v_int(end+1) = Inf;
v_dist = NaN(1,n_v_dist);
end
%% MODIFY
ball_pos(:,4) = 1;
ball_pos_new = ball_pos;
ball_v_new = ball_v;
GravityVec = GravityVec/sqrt(GravityVec*GravityVec'); %normalize
GravityField = g_grav*GravityVec;
GravityPlane = [GravityVec -GravityVec*GroundPoint'];
%CHECK
[~,kc,~] = check_overlap;
if kc > 0
fprintf(['ERROR: There is an overlap in the current sphere layout.'...
'\nPlease check ''ball_pos'' and try again.\n\n'])
return
end
%% PLOT and SIMULATION
%Generate Figure and Axes
t_dura = max([min([10^(round(log10(dt/N^2))),0.01]),1e-4]);
h_f = figure('name','Sphere Collider','menubar','none','numbertitle','off',...
'position',[45 280 600 450],'color',[1 1 1]);
h_a = axes('parent',h_f);
set(h_a,'box','on','projection','perspective','dataaspectratio',[1 1 1])
set(h_a,'xlim',axis_lim(1,:),'ylim',axis_lim(2,:),'zlim',axis_lim(3,:))
set(h_a,'xlimmode','manual','ylimmode','manual','zlimmode','manual')
set(h_a,'xgrid','on','ygrid','on','zgrid','on')
xlabel(h_a,'x')
ylabel(h_a,'y')
zlabel(h_a,'z')
rotate3d(h_a)
%Debugging Figure
if debug_mode
%Momentum & Energy
lin_color = [1 0 0;0 1 0;0 0 1;0 1 0;0 0 1;0 0 0];
lin_style = {'-','-','-',':',':','-'};
h_lin1 = zeros(1,6);
p_total = NaN(tlen_fig,3);
E_total = NaN(tlen_fig,3);
% figure1: Total Momentum and Energy
% subplot(2,1,1): Momentum, subplot(2,1,2): Energy;
h_fd1 = figure('name','Momentum and Energy','menubar','none',...
'numbertitle','off',...
'position',[728 380 600 350]); %[728 415 600 305]
h_ad1 = zeros(1,2);
E_total(1,:) = total_energy;
h_ad1(2) = subplot(2,1,2,'parent',h_fd1,'xlim',[0 dt*tlen_fig]);
title(h_ad1(2),'Energy in Total');
xlabel(h_ad1(2),'time interval')
ylabel(h_ad1(2),'Energy')
for dk = 4:6
h_lin1(dk) = line('xdata',dt*(1:tlen_fig),'ydata',E_total(:,dk-3),...
'parent',h_ad1(2),'color',lin_color(dk,:),'linestyle',lin_style{dk});
end %energy
p_total(1,:) = total_momentum;
h_ad1(1) = subplot(2,1,1,'parent',h_fd1,'xlim',[0 dt*tlen_fig]);
title(h_ad1(1),'Momentum in Total')
xlabel(h_ad1(1),'time interval')
ylabel(h_ad1(1),'Momentum')
for dk = 1:3
h_lin1(dk) = line('xdata',dt*(1:tlen_fig),'ydata',p_total(:,dk),...
'parent',h_ad1(1),'color',lin_color(dk,:),'linestyle',lin_style{dk});
end %momentum
%Pressure & Velocity Distribuion
press_total = zeros(1,tlen_fig); %for computing press_average
press_average = NaN(1,tlen_fig);
%pressure
h_lin2 = zeros(1,4); %1: current pressure; 2: average pressure;
%3: theoratical pressure;
%4: velocity distribution;
h_fd2 = figure('name','Pressure & Velocity Distribution',...
'menubar','none','numbertitle','off',...
'position',[728 55 600 285]); %[728 415 600 305]
h_ad2 = zeros(1,2);
h_ad2(1) = subplot(1,2,1,'parent',h_fd2,'xlim',[0 dt*tlen_fig]);
title(h_ad2(1),'Pressure Measurement');
xlabel(h_ad2(1),'time interval')
ylabel(h_ad2(1),'Pressure')
h_lin2(1) = line('xdata',dt*(1:tlen_fig),'ydata',press_total,...
'parent',h_ad2(1),'color',[0 0 1],'linestyle','-');
h_lin2(2) = line('xdata',dt*(1:tlen_fig),'ydata',press_average,...
'parent',h_ad2(1),'color',[0 0 0],'linestyle',':');
%velocity distribution
h_ad2(2) = subplot(1,2,2,'parent',h_fd2,'xlim',[0 v_max],'ylim',[0 1]);
title(h_ad2(2),'Velocity Distribution');
xlabel(h_ad2(2),'velocity')
ylabel(h_ad2(2),'Probability')
h_lin2(4) = line('xdata',v_int(1:n_v_dist),'ydata',v_dist,...
'parent',h_ad2(2),'color',[0 0 0],...
'linestyle',':','linewidth',1,...
'marker','o','markerfacecolor',[1 0 0.25],...
'markeredgecolor','none','markersize',4.5);
end
%Generate Balls
sph_n = 25;
[bx,by,bz] = sphere(sph_n);
ball_srf_data = cell(1,N); %{[(sph_n+1)x(sph_n+1)x3 matrix]}
sf_h = zeros(1,N);
%Initial xyz_data of balls
for ikb = 1:N
ball_srf_data{ikb} = zeros(sph_n+1,sph_n+1,3);
ball_srf_data{ikb}(:,:,1) = balls(ikb)*bx;
ball_srf_data{ikb}(:,:,2) = balls(ikb)*by;
ball_srf_data{ikb}(:,:,3) = balls(ikb)*bz;
sf_h(ikb) = surface(...
ball_srf_data{ikb}(:,:,1),...
ball_srf_data{ikb}(:,:,2),...
ball_srf_data{ikb}(:,:,3),...
'parent',h_a,...
'FaceColor',ball_color(ikb,1:3),'FaceAlpha',ball_color(ikb,4),...
'linestyle','none','visible','off');
end %Initial xyz_data of balls
%Initial plot
for ikb = 1:N
set(sf_h(ikb),...
'xdata',ball_srf_data{ikb}(:,:,1)+ball_pos(ikb,1),...
'ydata',ball_srf_data{ikb}(:,:,2)+ball_pos(ikb,2),...
'zdata',ball_srf_data{ikb}(:,:,3)+ball_pos(ikb,3),...
'visible','on')
end %Initial plot
%Plot Walls
%plane
h_pln = zeros(1,n_pln);
for ikw = 1:n_pln
h_pln(ikw) = patch(...
pln_pnt{ikw}(:,1),...
pln_pnt{ikw}(:,2),...
pln_pnt{ikw}(:,3),...
'k','parent',h_a,...
'FaceColor',[0.85 0.4 0.24],'facealpha',0.2,...
'linestyle','none');
end
%sphere
h_sph = zeros(1,n_sph);
sph_srf_data = cell(1,n_sph); %{[(sph_n+1)x(sph_n+1)x3 matrix]}
for ikw = 1:n_sph
sph_srf_data{ikw} = zeros(sph_n+1,sph_n+1,3);
sph_srf_data{ikw}(:,:,1) = wall_sph(ikw,1)*bx;
sph_srf_data{ikw}(:,:,2) = wall_sph(ikw,1)*by;
sph_srf_data{ikw}(:,:,3) = wall_sph(ikw,1)*bz;
h_sph(ikw) = surface(...
sph_srf_data{ikw}(:,:,1)+wall_sph(ikw,2),...
sph_srf_data{ikw}(:,:,2)+wall_sph(ikw,3),...
sph_srf_data{ikw}(:,:,3)+wall_sph(ikw,4),...
'parent',h_a,...
'FaceColor',[0.75 0.2 0.85],'facealpha',0.2,...
'linestyle','none');
end %plot wall_sphere
%% SIMULATION
%Start Simulation!
dt_pre = dt; % dt_p: time interval for prediction
dt_c = 0; % dt_c: cumulative time of collision
t = dt; % time
cnt_t = 1; % counter for time
cnt_c = 0; % counter for collision
cnt_error = 0; % counter for errors
impulse_sum = 0;
press_total_sum = 0;
while cnt_c < cldmax && cnt_t < tlen
cnt_c = cnt_c+1;
%predict and update position
ball_pos_new(:,1:3) = ball_pos(:,1:3)+dt_pre*ball_v;
%Collision Detection
[dt_int,kc,kc_w] = check_overlap;
if dt_int < inf %collide!
%processing
%update position and velocity
dt_c = dt_c+dt_int;
dt_pre = dt-dt_c;
ball_pos(:,1:3) = ball_pos(:,1:3)+dt_int*ball_v;
ball_v(kc,:) = ball_v_new(kc,:);
if kc_w > 0 %collide with walls
impulse_frame(kc_w) = impulse_frame(kc_w)+impulse_collide;
end
%ACCELERATION NOT UPDATED
else %No overlapping
%#(ONE FRAME)#
cnt_t = cnt_t+1;
t = t+dt;
for kb1 = 1:N
%update acceleration and velocity
a = [0 0 0];
for kb2 = [1:kb1-1,kb1+1:N]
norm_balls = ball_pos(kb2,1:3)-ball_pos(kb1,1:3); %kb1 -> kb2
dist_balls = sqrt(norm_balls*norm_balls');
a_coulomb_magn = -k_coulomb*balls(kb1,3)*balls(kb2,3)/dist_balls^2; %negative
norm_balls = norm_balls/dist_balls; %normalize
a = a+a_coulomb_magn*norm_balls;
end %Interaction between balls
for ksg = 1:n_grav
norm_src = ball_pos(kb1,1:3)-GravitySrc(ksg,2:4);
dist_src = sqrt(norm_src*norm_src');
a_grav_magn = G_newton*balls(kb1,2)*GravitySrc(ksg,1)/dist_src^2;
norm_src = norm_src/dist_src; %normalize
a = a+a_grav_magn*norm_src;
end %Gravity
for ksc = 1:n_coulomb
norm_src = ball_pos(kb1,1:3)-CoulombSrc(ksc,2:4);
dist_src = sqrt(norm_src*norm_src');
a_coulomb_magn = k_coulomb*balls(kb1,3)*CoulombSrc(ksc,1)/dist_src^2;
norm_src = norm_src/dist_src; %normalize
a = a+a_coulomb_magn*norm_src;
end %Coulomb force
%Gravity Field
a = a+GravityField-...
friction_factor/balls(kb1,2)*ball_v(kb1,:);
ball_v(kb1,:) = ball_v(kb1,:)+dt*a;
%DRAW
set(sf_h(kb1),...
'xdata',ball_srf_data{kb1}(:,:,1)+ball_pos_new(kb1,1),...
'ydata',ball_srf_data{kb1}(:,:,2)+ball_pos_new(kb1,2),...
'zdata',ball_srf_data{kb1}(:,:,3)+ball_pos_new(kb1,3))
end %update acceleration and velocity
%update position
ball_pos = ball_pos_new;
%move walls
for kw = 1:n_pln
vw = [w_pln_v{kw}{1}(t),w_pln_v{kw}{2}(t),w_pln_v{kw}{3}(t)];
pln_pnt{kw} = pln_pnt{kw}+([1;1;1;1]*vw)*dt;
set(h_pln(kw),...
'xdata',pln_pnt{kw}(:,1),...
'ydata',pln_pnt{kw}(:,2),...
'zdata',pln_pnt{kw}(:,3))
wall_pln(kw,4) = wall_pln(kw,4)-(vw*wall_pln(kw,1:3)')*dt;
end
for kw = 1:n_sph
vw = [w_sph_v{kw}{1}(t),w_sph_v{kw}{2}(t),w_sph_v{kw}{3}(t)];
wall_sph(kw,2:4) = wall_sph(kw,2:4)+vw*dt;
set(h_sph(kw),...
'xdata',sph_srf_data{kw}(:,:,1)+wall_sph(kw,2),...
'ydata',sph_srf_data{kw}(:,:,2)+wall_sph(kw,3),...
'zdata',sph_srf_data{kw}(:,:,3)+wall_sph(kw,4))
end
pause(t_dura);
%RESET dt_p and dt_c
dt_pre = dt; %this is important!
dt_c = 0;
%DEBUG
if debug_mode
%test for conservation of momentum
idx = rem(cnt_t-1,tlen_fig)+1;
p_total(idx,:) = total_momentum;
E_total(idx,:) = total_energy;
for dk = 1:3
set(h_lin1(dk),'ydata',p_total([idx+1:tlen_fig,1:idx],dk))
end %momentum
for dk = 4:6
set(h_lin1(dk),'ydata',E_total([idx+1:tlen_fig,1:idx],dk-3))
end %energy
%measure pressure & velocity
%RESET pressure container
idx_press = rem(cnt_t-1,dn_press)+1;
impulse_sum = impulse_sum-sum(impulse_wall(idx_press,:),2)+...
sum(impulse_frame,2);
impulse_wall(idx_press,:) = impulse_frame;
impulse_frame(:) = 0;
press_total_sum = press_total_sum-press_total(idx);
press_total(idx) = impulse_sum/press_factor;
press_total_sum = press_total_sum+press_total(idx);
press_average(idx) = press_total_sum/tlen_fig;
set(h_lin2(1),'ydata',press_total([idx+1:tlen_fig,1:idx]))
set(h_lin2(2),'ydata',press_average([idx+1:tlen_fig,1:idx]))
%compute velocity distribution
velocity_distribution
set(h_lin2(4),'ydata',v_dist)
end
end
end %Simulation
%% REPORT
fprintf('Number of ERRORs: %d\n',cnt_error)
%% FUNCTIONS
function p_total = total_momentum
p_total = zeros(1,3);
p_total(1) = balls(:,2)'*ball_v(:,1);
p_total(2) = balls(:,2)'*ball_v(:,2);
p_total(3) = balls(:,2)'*ball_v(:,3);
end %compute total momentum
function E_total = total_energy
E_total = zeros(1,3);
E_total(1) = balls(:,2)'*sum(ball_v.^2,2)/2;
E_total(2) = balls(:,2)'*abs(ball_pos*GravityPlane')*g_grav;
E_total(3) = E_total(1)+E_total(2);
end %compute total energy
function [dt_int,kc,kc_w] = check_overlap
dt_int = inf;
kc = -1; %unnecessary, actually :)
kc_w = -1;
for fk1 = 1:N
%update acceleration
%detect possible collision between ball and walls
for fkw = 1:n_pln
dist_plane = abs(wall_pln(fkw,:)*ball_pos_new(fk1,:)');
iscross = (wall_pln(fkw,:)*ball_pos_new(fk1,:)')*...
(wall_pln(fkw,:)*ball_pos(fk1,:)') < 0; %bug fixed
if dist_plane < balls(fk1,1) || iscross %!!
vbn = ball_v(fk1,:)*wall_pln(fkw,1:3)'*wall_pln(fkw,1:3);
vwn = ([w_pln_v{fkw}{1}(t),w_pln_v{fkw}{2}(t),w_pln_v{fkw}{3}(t)]*...
wall_pln(fkw,1:3)')*wall_pln(fkw,1:3);
dist_plane_old = abs(wall_pln(fkw,:)*ball_pos(fk1,:)');
dt_int_c = (dist_plane_old-balls(fk1,1))/...
sqrt((vbn-vwn)*(vbn-vwn)');
if dt_int_c < dt_int
dt_int = dt_int_c;
ball_v_new(fk1,:) = ball_v(fk1,:)-(1+Cr)*vbn+2*vwn;
kc = fk1;
impulse_collide = abs(ball_v_new(fk1,:)-ball_v(fk1,:))*...
balls(fk1,2);
impulse_collide = sqrt(impulse_collide*impulse_collide'); %scalarize
kc_w = fkw;
end
end
end %detect possible collision between ball and walls
for fkw = 1:n_sph
sph_norm = ball_pos_new(fk1,1:3)-wall_sph(fkw,2:4);
sph_norm_old = ball_pos(fk1,1:3)-wall_sph(fkw,2:4);
dist_sphere = sqrt(sph_norm*sph_norm');
dist_sphere_old = sqrt(sph_norm_old*sph_norm_old');
iscross = ((dist_sphere-wall_sph(fkw,1))*...
(dist_sphere_old-wall_sph(fkw,1))) < 0; %bug fixed
if abs(dist_sphere-wall_sph(fkw,1)) < balls(fk1,1) || iscross %!!
sph_norm = sph_norm/dist_sphere; %normalize
vbn = (ball_v(fk1,:)*sph_norm')*sph_norm;
vwn = ([w_sph_v{fkw}{1}(t),w_sph_v{fkw}{2}(t),w_sph_v{fkw}{3}(t)]*...
sph_norm')*sph_norm;
dt_int_c = (abs(dist_sphere_old-wall_sph(fkw,1))...
-balls(fk1,1))/sqrt((vbn-vwn)*(vbn-vwn)');
if dt_int_c < dt_int
dt_int = dt_int_c;
ball_v_new(fk1,:) = ball_v(fk1,:)-(1+Cr)*vbn+2*vwn;
kc = fk1;
impulse_collide = abs(ball_v_new(fk1,:)-ball_v(fk1,:))*...
balls(fk1,2);
impulse_collide = sqrt(impulse_collide*impulse_collide'); %scalarize
kc_w = fkw+n_pln; %!!
end
end
end %detect possible collision between ball and walls
for fk2 = fk1+1:N
%detect possible collision between balls
norm_balls = ball_pos_new(fk2,1:3)-ball_pos_new(fk1,1:3); %kb1 -> kb2
if norm_balls*norm_balls' < (balls(fk1,1)+balls(fk2,1))^2
v2_1 = ball_v(fk2,:)-ball_v(fk1,:); %relative velocity
r2_1 = ball_pos(fk2,1:3)-ball_pos(fk1,1:3); %relative position
d_sq = (det([r2_1([1,2]);v2_1([1,2])])^2+det([r2_1([2,3]);v2_1([2,3])])^2+...
det([r2_1([3,1]);v2_1([3,1])])^2)/(v2_1*v2_1');
d_c = sqrt(r2_1*r2_1'-d_sq)-sqrt((balls(fk1,1)+balls(fk2,1))^2-d_sq);
dt_int_c = abs(d_c)/sqrt(v2_1*v2_1');
%DEBUG
if d_sq >= (balls(fk1,1)+balls(fk2,1))^2
set(sf_h([fk1,fk2]),'facecolor',[0.75 0.25 0.75])
disp('ERROR: "d" >= R1+R2')
cnt_error = cnt_error+1;
end
if dt_int_c < 0
set(sf_h([fk1,fk2]),'facecolor',[0.75 0.25 0.75])
disp('ERROR: "dt_int_c" < 0')
cnt_error = cnt_error+1;
end
ball_pos_new(fk1,1:3) = ball_pos(fk1,1:3)+dt_int_c*ball_v(fk1,1:3);
ball_pos_new(fk2,1:3) = ball_pos(fk2,1:3)+dt_int_c*ball_v(fk2,1:3);
norm_balls_c = (ball_pos_new(fk2,1:3)-ball_pos_new(fk1,1:3))/...
(balls(fk1,1)+balls(fk2,1)); %normalize
if dt_int_c < dt_int
dt_int = dt_int_c;
vn1 = (ball_v(fk1,:)*norm_balls_c')*norm_balls_c;
vn2 = (ball_v(fk2,:)*norm_balls_c')*norm_balls_c;
va = (balls(fk1,2)*vn1+balls(fk2,2)*vn2+balls(fk2,2)*...
Cr*(vn2-vn1))/(balls(fk1,2)+balls(fk2,2));
vb = (balls(fk1,2)*vn1+balls(fk2,2)*vn2+balls(fk1,2)*...
Cr*(vn1-vn2))/(balls(fk1,2)+balls(fk2,2));
ball_v_new(fk1,:) = ball_v(fk1,:)-vn1+va;
ball_v_new(fk2,:) = ball_v(fk2,:)-vn2+vb;
kc = [fk1 fk2];
end
end
end %detect possible collision between balls
end %loop contains every ball
end %detect possible collisions
function wall_pln = generate_walls
wall_pln = zeros(n_pln,4);
for fkw = 1:n_pln
vec1 = pln_pnt{fkw}(2,:)-pln_pnt{fkw}(1,:);
vec2 = pln_pnt{fkw}(3,:)-pln_pnt{fkw}(1,:);
vecn = cross(vec1,vec2);
vecn = vecn/sqrt(vecn*vecn'); %normalize
wall_pln(fkw,:) = [vecn -vecn*pln_pnt{fkw}(1,:)'];
end
end %generate walls(wall_pln) from plane_pnt
function velocity_distribution
v_dist(:) = 0;
ball_v_magn = sort(sum(sqrt(ball_v.^2),2));
fkint = 2;
for fk = 1:N
while fkint <= n_v_dist+1
if ball_v_magn(fk) < v_int(fkint)
v_dist(fkint-1) = v_dist(fkint-1)+1;
break
else
fkint = fkint+1;
end
end
end
v_dist = v_dist/N;
end
end
