- #1

eltrinco

- 8

- 0

## Homework Statement

I have created random mass and try to make them go around each other . somehow i couldn't find where all the code goes wrong. how do i fix this

## The Attempt at a Solution

I have coded this. but still the file can't be run. i don't know where to fix this .

% Creating initial positions/velocities (Randomly)

clear all

close all

clc

n_masses=10;

m=rand(1,n_masses) %RANDOM masses fore each body IF YOU WANT

% the same mass for all use instead :

%m=k*ones(1,n_masses) where k is the mass

x=randn(1,n_masses); %X position

y=randn(1,n_masses); %Y position

vx(1,:)=randn(1,n_masses); %X component velocity (initial)

vy(1,:)=randn(1,n_masses) %Y component velocity (initial)

scatter(x,y,round(100.*m),'filled','blue') % different masses, different sizes

title('Initial configuration')

figure

% Calculating force matrix (F{i,j} means Force exterted by the i-th

% particle on the j-th.

%when j=i -> force is zero

[F,F_totalx,F_totaly]=calculating_F_matrix(x,y,m);

%Evolution

dt=0.05; %interval

n_steps=40; %number of timesteps

%starting positions for the simulation

xx{1}=x; %Will contain the succession of the X positions

yy{1}=y; %Will contain the succession of the Y positions

for s=1:n_steps+1

disp(['Simulating at time t=' num2str(dt*(s-1))]);

for i=1:n_masses

[F,F_totalx,F_totaly]=calculating_F_matrix(xx{s},yy{s},m);

%initial accelerations (fixed)

ax(s,:)=F_totalx;

ay(s,:)=F_totaly;

%Updating velocities

vx(s+1,i)=vx(s,i)+dt*ax(s,i);

vy(s+1,i)=vy(s,i)+dt*ay(s,i);

%Updating positions

xx{s+1}(i)=xx{s}(i)+dt*xx{s}(i);

yy{s+1}(i)=yy{s}(i)+dt*yy{s}(i);

end

end

% Graphics (plot the evolution of the particles)

%find max x and y to correctly set the reference for the plot

for i=1:length(xx)

xm(i)=max((abs(xx{i})));

ym(i)=max(abs(yy{i}) );

end

axis (1.1.*[-max(xm) max(xm) -max(ym) max(ym) ]); % set axis to see all particles

for s=1:n_steps

scatter(xx{s},yy{s},round(100.*m),'filled','red')

axis (1.2.*[-max(xm) max(xm) -max(ym) max(ym) ]); % set axis to see all particles

title(['Time t=' num2str(dt*(s-1))])

pause(0.15)

end

function [F,F_totalx,F_totaly]=calculating_F_matrix(2,3,5)

n_masses=length(m);

G=6.67*10^(-11); % Universal gravity constant

for i=1:n_masses %influenced index

for j=1:n_masses %index generating field

if j~=i % excluding j=i

Xdistance=abs(x(i)-x(j));

Ydistance=abs(y(i)-y(j));

%distance=sqrt((x(i)-x(j)).^2 )%+(y(i)-y(j)).^2); %distance i-j

% dij=(x(i)-x(j)) %unitary vector dij

F{i,j}(1)=G*m(i)*m(j)/(Xdistance);

F{i,j}(2)=G*m(i)*m(j)/(Ydistance);

else

F{i,j}(1)=0;

F{i,j}(2)=0;

end

end

end

% Calculating the total force (per unit of mass m_i) acting on i-th particle (by components)

F_totalx(1)=0;

F_totaly(1)=0;

for i=1:n_masses

Z=cell2mat(F(i,:)'); %vector of all forces acting on the i-th particle

F_totalx(i)=sum(Z(:,1))./m(i);

F_totaly(i)=sum(Z(:,2))./m(i);

end