matlab part of a program for nbody sim

1)

Instead of calling the circle function once for each body, it would be nice to let a new function call the circle function for all the bodies.

1. Xs: this is the vector of all the bodies' x coordinates
2. Ys: this is the vector of all the bodies' y coordinates
3. Masses: the vector of masses
4. color: the color to plot all the bodies

Create a new function called plotBodies that has these parameters:
Determine the number of bodies that are to be plotted (it's the length of either the Xs, Ys, or Masses vectors) and store it in a variable.

Start a for loop that goes from 1 to the number of bodies. Use a loop control variable like b or body, because you're iterating over the number of bodies. I'll refer to the variable as b.

Call the circle function. You need to supply 4 arguments when you call that function: the x coordinate, the y coordinate, the radius, and the color.

Recall how you tested the circle function:
circle(Xs(1), Ys(1), Masses(1), 'k');
This function call draws the circle for mass 1. Now you want to drawn the circle for mass 1, you want to draw a circle for mass b.

End the for loop

so far I got this :

Masses, Xs, and Ys are set numbers given in my program

my circle function is

function circle(x,y,Radius,color)

Theta = 0 : 0.1 : 2 * pi;
X = Radius * cos ( Theta )+ x;
Y = Radius * sin ( Theta )+ y;
plot ( X, Y, color )
end

and I tried the plotBodies but was lost:

function plotBodies(Xs,Ys,Masses,color)
n=length(Ys);

for b=(1:1:n)
   circle1(Xs(b), Ys(b), Masses(b), 'k'));

end

2)
This function has these parameters

1. A body number. I called mine b.
2. Xs: these are the x positions of all the bodies
3. Ys: these are the y positions of all the bodie
4. Dxs: these are the x velocities of all the bodies
5. Dys: these are the y velocities of all the bodies
6. Masses

This function returns a vector of two values: dx and dy, which is the new velocity of the body.
To calculate the total force on a body, you must add up all the separate forces on a body by all bodies other than itself.
Start a for loop that does this:

For each body number from 1 to the number of bodies:
.
○ If the body number is not equal to b:
Calculate the force on the body by calling calculateForce.
Calculate the x & y accelerations. The acceleration is the force divided by the mass. The mass is the mass of body b.
Add the x & y accelerations to the x & y velocities of body b.

I calculate force using this program

function calculateForce(body1,body2,Xs,Ys,Masses)
dx=Xs(body1)-Xs(body2);
dy=Ys(body2)-Ys(body1);
r=sqrt(dx^2+dy^2);

Radi=(Masses/2);
if r <(Radi*2)
r=Radi+Radi;
disp(r)
else
disp('invalid')
end

f= (6.67 * 10^-11)*(body1*body2)/r^2;
disp(f)

FA=atan2(dy,dx);
disp(FA)
fx=cos(FA);
disp(fx)
fy=sin(FA);
disp(fy)
end

I only got to

function accelerateBody(b,Xs,Ys,Dxs,Dys,Masses)=(Dy,Dx)

and now im stuck on this part

Answer 1

[arushi]

Hi John,

You want to plot circles for all bodies using the circle function. Here's how you can do that:

function plotBodies(Xs, Ys, Masses, color)
    n = length(Ys); % Determine the number of bodies
    for b = 1:n % Loop through each body
        circle(Xs(b), Ys(b), Masses(b), color); % Call the circle function
    end
end

function circle(x, y, Radius, color)
    Theta = 0 : 0.1 : 2 * pi;
    X = Radius * cos(Theta) + x;
    Y = Radius * sin(Theta) + y;
    plot(X, Y, color);
    hold on; % Ensure all circles are plotted on the same figure
end

Complete accelerateBody Function -

This function calculates the new velocity of a body based on the forces exerted by all other bodies.

function [Dx, Dy] = accelerateBody(b, Xs, Ys, Dxs, Dys, Masses)
    G = 6.67 * 10^-11; % Gravitational constant
    Dx = Dxs(b);
    Dy = Dys(b);
    for i = 1:length(Xs)
        if i ~= b
            [fx, fy] = calculateForce(b, i, Xs, Ys, Masses, G);
            ax = fx / Masses(b); % Acceleration in x direction
            ay = fy / Masses(b); % Acceleration in y direction
            Dx = Dx + ax; % Update velocity in x direction
            Dy = Dy + ay; % Update velocity in y direction
        end
    end
end

function [fx, fy] = calculateForce(body1, body2, Xs, Ys, Masses, G)
    dx = Xs(body2) - Xs(body1);
    dy = Ys(body2) - Ys(body1);
    r = sqrt(dx^2 + dy^2);
    Radi = (Masses(body1) / 2) + (Masses(body2) / 2); % Sum of radii

    if r < Radi
        r = Radi; % Adjust distance to avoid division by zero or too small distances
    end

    f = G * (Masses(body1) * Masses(body2)) / r^2; % Gravitational force
    FA = atan2(dy, dx); % Angle of force
    fx = f * cos(FA); % Force in x direction
    fy = f * sin(FA); % Force in y direction
end

Hope this helps.