Attraction Force between 2D balls

Question

I have a simulation with multiple circles moving in 2D space, with elastic collisions between them.

I'd like to add an attraction force between particles, so that particles move towards other particles depending on mass, etc. How would I go about this?

My collision management function looks like this:

void manageCollision(Particle particleA, Particle particleB)
{
    float distanceX = particleA.Position.X - particleB.Position.X;
    float distanceY = particleA.Position.Y - particleB.Position.Y;
    double collisionAngle = Math.Atan2(distanceY, distanceX);
    double pA_magnitude = Math.Sqrt(particleA.Velocity.X * particleA.Velocity.X + particleA.Velocity.Y * particleA.Velocity.Y);
    double pB_magnitude = Math.Sqrt(particleB.Velocity.X * particleB.Velocity.X + particleB.Velocity.Y * particleB.Velocity.Y);
    double pA_direction = Math.Atan2(particleA.Velocity.Y, particleA.Velocity.X);
    double pB_direction = Math.Atan2(particleB.Velocity.Y, particleB.Velocity.X);
    double pA_newVelocityX = pA_magnitude * Math.Cos(pA_direction - collisionAngle);
    double pA_newVelocityY = pA_magnitude * Math.Sin(pA_direction - collisionAngle);
    double pB_newVelocityX = pB_magnitude * Math.Cos(pB_direction - collisionAngle);
    double pB_newVelocityY = pB_magnitude * Math.Sin(pB_direction - collisionAngle);
    double pA_finalVelocityX = ((particleA.Mass - particleB.Mass) * pA_newVelocityX + (particleB.Mass + particleB.Mass) * pB_newVelocityX) / (particleA.Mass + particleB.Mass);
    double pB_finalVelocityX = ((particleA.Mass + particleA.Mass) * pA_newVelocityX + (particleB.Mass - particleA.Mass) * pB_newVelocityX) / (particleA.Mass + particleB.Mass);
    double pA_finalVelocityY = pA_newVelocityY;
    double pB_finalVelocityY = pB_newVelocityY;
    particleA.Velocity = new Vector2((float)(Math.Cos(collisionAngle) * pA_finalVelocityX + Math.Cos(collisionAngle + Math.PI / 2) * pA_finalVelocityY), (float)(Math.Sin(collisionAngle) * pA_finalVelocityX + Math.Sin(collisionAngle + Math.PI / 2) * pA_finalVelocityY));
    particleB.Velocity = new Vector2((float)(Math.Cos(collisionAngle) * pB_finalVelocityX + Math.Cos(collisionAngle + Math.PI / 2) * pB_finalVelocityY), (float)(Math.Sin(collisionAngle) * pB_finalVelocityX + Math.Sin(collisionAngle + Math.PI / 2) * pB_finalVelocityY));
}

Each ball or particle spawns with a random mass and radius.

The function is called within an update type of method, like this:

Vector2 globalGravity = new Vector2(0f, gravityScale / 6000);

    for (int i = 0; i < particles.Count(); i++)
{
    particles[i].Update((float)updateTimer.Interval, globalGravity);
    Vector2 position = particles[i].Position;
    Vector2 velocity = particles[i].Velocity;
    collisionWallCheck(ref position, ref velocity, particles[i].Radius);
    particles[i].Position = position;
    particles[i].Velocity = velocity;


    Particle pA = particles[i];
    for (int k = i + 1; k < particles.Count(); k++)
    {
        Particle pB = particles[k];
        Vector2 delta = pA.Position - pB.Position;
        float dist = delta.Length();

        if (dist < particles[i].Radius + particles[k].Radius && !particles[i].Colliding && !particles[k].Colliding)
        {
            particles[i].Colliding = true;
            particles[k].Colliding = true;
            manageCollision(particles[i], particles[k]);
            particles[i].initColorTable(); // Upon collision, change the color
            particles[k].initColorTable();
            totalCollisions++;
        }
        else
        {
            particles[i].Colliding = false;
            particles[k].Colliding = false;
        }
    }
}

I'm storing the initial position, velocity and masses of each ball.

What I apparently need to do, and don't know how to implement, is:

Calculate the magnitude and direction of the gravitational force.
Knowing the force, you can calculate the acceleration of each body.
Knowing the acceleration you can calculate the new velocity.
Knowing the velocity you can calculate the new position.

I'm shaky with the equations for it essentially, and I'd like to start off by making an attraction force between just two balls.

Using Steven's suggestion, this is the new integrated code.

void updateTimer_Tick(object sender, EventArgs e)
{
    const double G = 6.67398 * 0.00000000001;

    for (int i = 0; i < particles.Count(); i++)
    {
        double sumX = 0;
        double sumY = 0;

        Particle pA = particles[i];
        for (int k = i + 1; k < particles.Count(); k++)
        {
            Particle pB = particles[k];
            Vector2 delta = pA.Position - pB.Position;
            float dist = delta.Length();

            if (dist < particles[i].Radius + particles[k].Radius && !particles[i].Colliding && !particles[k].Colliding)
            {
                particles[i].Colliding = true;
                particles[k].Colliding = true;
                manageCollision(particles[i], particles[k]);
                particles[i].initColorTable();
                particles[k].initColorTable();
                totalCollisions++;
                particles[i].Colliding = false;
                particles[k].Colliding = false;
            }
            else
            {
                double distanceX = particles[i].Position.X - particles[k].Position.X;
                double distanceY = particles[i].Position.Y - particles[k].Position.Y;
                double r = Math.Sqrt(Math.Pow(distanceX, 2) + Math.Pow(distanceY, 2));
                double force = G * particles[i].Mass * particles[k].Mass / (r * r);
                double theta = Math.Tan(distanceY / distanceX);
                sumX += force * Math.Cos(theta);
                sumY += force * Math.Sin(theta);
                particles[i].Colliding = false;
                particles[k].Colliding = false;
            }
        }
        double netForce = Math.Sqrt(Math.Pow(sumX, 2) + Math.Pow(sumY, 2));
        double a = netForce / particles[i].Mass;
        double aTheta = Math.Tan(sumY / sumX);

        // Here we get accelerations for X and Y.  You can probably figure out velocities from here.
        double aX = a * Math.Cos(aTheta);
        double aY = a * Math.Sin(aTheta);
        Vector2 accel = new Vector2((float)aX, (float)aY);

        particles[i].Update((float)updateTimer.Interval, accel);
        //particles[i].Update((float)updateTimer.Interval, globalGravity);
        Vector2 position = particles[i].Position;
        Vector2 velocity = particles[i].Velocity;
        collisionWallCheck(ref position, ref velocity, particles[i].Radius);
        particles[i].Position = position;
        particles[i].Velocity = velocity + accel;
    }
    Draw();
}

The Update function for the particles is simple, and before it used a global gravity Vector which was 0,0.

        public void Update(float timeStep, Vector2 gravity)
        {
            velocity = velocity + timeStep * gravity;
            position = position + timeStep * velocity;
        }

I'm now unsure how to deal with the cases of 0.

Steve Westbrook Steve Westbrook · Accepted Answer · 2013-03-21T01:59:35

Start by calculating the force of gravity acting on each object. This is given by

F = Gm1m2/r*r

where m1 and m2 are the masses of two objects, G is the gravitational constant, and r is the distance between the two objects.

Now, r is a vector, so you may want to split this up into separate components - Fx and Fy. You can do this as follows:

Fx = F * cos(theta)
Fy = F * sin(theta)

For each mass, calculate the force of gravity acting on it and every other object. Sum the vectors to get the net force of gravity. (Note - that link is available for your interest, but takes a long time to get to the point). At this point you will have a net force on each object, from which you can calculate acceleration. Here's the code to get to this point:

const double G = 6.67398 * 0.00000000001;

for (int i = 0; i < particles.Count(); i++)
{
    double sumX = 0;
    double sumY = 0;

    for (int j = 0; j < particles.Count(); j++)
    {
        // Don't add attraction to self
        if (i == j)
            continue;

        double distanceX = particles[i].Position.X - particles[j].Position.X;
        double distanceY = particles[i].Position.Y - particles[j].Position.Y;
        double r = Math.Sqrt(Math.Pow(distanceX, 2) + Math.Pow(distanceY, 2));
        double force = G * particles[i].Mass * particles[j].Mass / (r * r);
        double theta = Math.Tan(distanceY / distanceX);
        sumX += force * Math.Cos(theta);
        sumY += force * Math.Sin(theta);
    }

    double netForce = Math.Sqrt(Math.Pow(sumX, 2) + Math.Pow(sumY, 2));
    double a = netForce / particles[i].Mass;
    double aTheta = Math.Tan(sumY / sumX);

    // Here we get accelerations for X and Y.  You can probably figure out velocities from here.
    double aX = a * Math.Cos(aTheta);
    double aY = a * Math.Sin(aTheta);
}

NOTES

This doesn't take stuff like 0-values into account - you'll have to clean up this code to deal with special cases before it will run without crashing.

Don't update any positions until you've calculated all the forces, or else you'll be off for later elements in the list.

Another thing worth noting: This algorithm is O(n^2), so if you have more than a few bodies it's going to take a lot of crunching. That's unfortunately just the way it is; if you find a fast way to calculate gravitational attraction for a large number of bodies, you should probably call NASA.

Depending on your coordinate system, you may find the y-vectors getting reversed. This is because Euclidean geometry thinks of positive values of y as "going up" whereas programmers tend to measure y in positive units "going down" from the top of the screen. This can play havoc with your angles and things.

Attraction Force between 2D balls

2 Answers