The following gets the code working, but I've no idea if the results are meaningful!!
n_part=4;
kn=5;
kt=2/7*kn;
m=0.3;
g=9.81;
rad(1:n_part)=0.5;
v_init=0.5;
y=0;
% testing i_particle=1:n_part;
% testing j_particle=n_part+1:2*n_part;
position_x(1,1:n_part)=0.05;
theta=2*pi*rand(size(position_x));
velocity_x=v_init*sin(theta);
position_y(1,1:n_part)=0.05;
velocity_y=v_init*cos(theta);
timestep=100;
dt=0.001;
acceleration_x=zeros(timestep,n_part); %%%%%%%%%%%%%%
acceleration_y=zeros(timestep,n_part); %%%%%%%%%%%%%%
Fn=0;
Fn_i=0;
Fn_j=0;
v_half_x=zeros(timestep,n_part); %%%%%%%%%%%%%%
v_half_y=zeros(timestep,n_part); %%%%%%%%%%%%%%
for n=1:n_part
for k=2:timestep
% position_x
v_half_x(k,n)=velocity_x(k-1,n)+0.5*dt*acceleration_x(k-1,n);
position_x(k,n)=position_x(k-1,n)+v_half_x(k-1,n)*dt;
% position_y
v_half_y(k,n)=velocity_y(k-1,n)+0.5*dt*acceleration_y(k-1,n);
position_y(k,n)=position_y(k-1,n)+v_half_y(k-1,n)*dt;
for i=1:n_part
for j=i+1:n_part
%if i>j %%%%% i CANNOT be greater than j as you %%%%%%
%%%%% set j to be i+1 upwards! %%%%%%
lx=position_x(k-1,i)-position_x(k-1,j);
ly=position_y(k-1,i)-position_y(k-1,j);
root_xy=sqrt(ly^2+ly^2);
% force calculation
Fn=kn*root_xy^1.5;
Fn_i=Fn_i+Fn;
Fn_j=Fn_j+Fn;
% acceleration term
acceleration_x(k,:)=Fn_i./m;
acceleration_y(k,:)=Fn_j./m;
% end
end
end
velocity_y(k,n)=v_half_y(k-1,n)+0.5*dt*acceleration_y(k-1,n);
velocity_x(k,n)=v_half_x(k-1,n)+0.5*dt*acceleration_x(k-1,n);
end
end
