Elastic Collision of Circles in XNA

[Tharwen]

I've spent too many days trying to make this work now, so I've come to this forum for some much-needed aid.

Basically, I've been experimenting with some simple physics in XNA 4.0, and I've become completely stuck while trying to make objects rebound correctly. I've tried to use my own knowledge of mechanics and failed, so I eventually resorted to looking up formulae on Wikipedia. These aren't working. My circular objects seem to transfer their momentum only if they feel like it, and often add a few radians to their final directions just to make it look cooler (or something... I didn't ask).

Anyway, here's the code for the collision:

Code: Select all

        public void Collide(Item Item1, Item Item2)
        {
            Vector2 TotalSpeed = Item1.Speed - Item2.Speed;
            float u1 = GetHypoteneuse(TotalSpeed);

            float m1 = Item1.Mass;
            float m2 = Item2.Mass;

            //The location of the overall centre of mass before(1) and after(2) the collision.
            Vector2 MassCentre1 = ((Item1.Position * m1) + (Item2.Position * m2)) / (m1 + m2);
            Vector2 MassCentre2 = (((Item1.Position + Item1.Speed) * m1) + ((Item2.Position + Item2.Speed) * m2)) / (m1 + m2);

            //The angle through which the overall centre of mass was deflected
            float MassDeflection = XYToBearing(MassCentre2 - MassCentre1);
           
            //The angles from Item1's initial direction through which Item1 and Item2 have been deflected
            float Deflection1 = (float)Math.Atan((m2 * Math.Sin(MassDeflection))/(m1 + (m2 * Math.Cos(MassDeflection))));
            float Deflection2 = (float)(-1 * MassDeflection);

            //The final velocities of Item1 and Item2
            float v1 = (float)(u1 * ((Math.Sqrt(Math.Pow(m1, 2) + Math.Pow(m2, 2) + (2 * m1 * m2 * Math.Cos(MassDeflection)))) / (m1 + m2)));
            float v2 = (float)(u1 * ((2 * m1) / (m1 + m2)) * Math.Sin(MassDeflection / 2));

            Item1.Speed = new Force(Deflection1, v1).ForcetoXY() - Item2.Speed;
            Item2.Speed += new Force(Deflection2, v2).ForcetoXY();

        }


Those lines that start "float Deflection1/2" and "float v1/2" are the formulae copied from here.

Some clarity, if needed:

GetHypoteneuse(Vector2 Value) returns Sqrt(Value.X^2 + Value.Y^2).
ForceToXY() converts a Force object (with magnitude and angular direction) into a Vector2 (I've used my Force type here for convenience, not because I'm simulating actual forces)
XYToBearing(Vector2 Value) converts Value into a bearing, with 0 pointing at the top of the screen and 2(PI) pointing at the bottom (I'm using radians because that's XNA's default).

Thanks to anyone who helps

[Robert'); DROP TABLE *;]

I'm still reading through your code, but I managed to find the following in one of my own projects, which AFAIK, works and does what you're trying to do.

Code: Select all

struct Vector2
{ // XNA provides one for you, but I was using standard Forms.
   public static Vector2 Rotate(Vector2 V, double Radians)
   {
      double newX = V.X * Math.Cos(Radians) - V.Y * Math.Sin(Radians);
      double newY = V.X * Math.Sin(Radians) + V.Y * Math.Cos(Radians);

      return new Vector2(newX, newY);
   }

   public double X { get; set; }
   public double Y { get; set; }
   public double Length
   {
      get
      {
         return Math.Sqrt(X * X + Y * Y);
      }
   }
   public double Angle
   {
      get
      {
         return Math.Atan(Y / X);
      }
   }
   
   /* Vector arithmetic overloads */
}   

private void HandleCollision( /* ... */)
{
/*
   Establishes Item and Other as Particles, which have mass, position and velocity. Position and velocity are vectors, mass is a double. Momentum is a Vector2, defined as (mass * velocity)
*/
   Vector2 Offset = Other.Position - Item.Position;

   double Angle = Offset.Angle;
   
   // Momentum in terms of perpendicular/parrarel to collision normal.
   Vector2 firstMomentum, secondMomentum;

   // Spin the coodinate system, so that (X,Y) corresponds to perpendicular, parrarel to collision.
   firstMomentum = Vector2.Rotate(Item.Momentum, Angle);
   secondMomentum = Vector2.Rotate(Other.Momentum, Angle);

   double temp = firstMomentum.Y;
   firstMomentum.Y = secondMomentum.Y;
   secondMomentum.Y = temp;

   Item.Momentum = Vector2.Rotate(firstMomentum, -Angle);
   Other.Momentum = Vector2.Rotate(secondMomentum, -Angle);   
}

Dunno whether that helps or not.

EDIT: Are you sure that the line is not supposed to be

Code: Select all

float Deflection2 = (float)(-1 * Deflection1);

[Tharwen]

Sorry about ignoring your post; it was quite helpful, so thanks

Although, when I tried directly implementing your code, it just made them reverse their directions (as if every collision was head-on) while still transferring momentum perfectly. If it worked in your program, that seems to suggest a strange problem somewhere else in my code...

(And the suggestion to change MassDeflection to Deflection1 made almost no difference, unfortunately)

[Robert'); DROP TABLE *;]

That's what I would expect to happen if you were swapping the momentum vectors, rather than just one component of each of them. (Are XNA vectors mutable?) Can you post what you've got at the moment?

[Tharwen]

Hi again,

Sorry about not responding; I've had loads of work to do and stuff, so I haven't done much on this. The problem is that, whenever I come back to it, I decide that my code is completely broken and try something new. Anyway, the code I have right now is this:

Code: Select all


            Vector2 Offset = Item2.Position - Item1.Position;
            float OffsetAngle = XYToBearing(Offset);

            float Angle1 = XYToBearing(Offset) - XYToBearing(Item1.Velocity);
            float Angle2 = XYToBearing(Offset) - XYToBearing(Item2.Velocity);

            Vector2 Momentum1 = AngleComponent(Item1.Momentum, OffsetAngle);
            Vector2 Momentum2 = AngleComponent(Item2.Momentum, OffsetAngle);

            Item1.Momentum += Momentum2;
            Item2.Momentum += Momentum1;
            Item1.Momentum -= Momentum1;
            Item2.Momentum -= Momentum2;


, and the XYToBearing function (where I suspect the problem may be) is this:

Code: Select all

        public float XYToBearing(Vector2 Vector)
        {
            float Temp = (float)(Math.Atan2(Vector.Y, Vector.X));
                return Temp + (float)(Math.PI / 2);
        }


Angle1 and Angle2 are supposed to be the angles between each item's direction and the bearing of 1 from 2.

Another problem with this code is that, for some reason, the objects sometimes lock together inside each other, with 1 slightly above and to the left of 2.

XNA vectors are mutable, and also, why should it just swap one component of each vector?

[Robert'); DROP TABLE *;]

Code: Select all

        public float XYToBearing(Vector2 Vector)
        {
            float Temp = (float)(Math.Atan2(Vector.Y, Vector.X));
                return Temp + (float)(Math.PI / 2);
        }


I think adding 90 degrees is messing up your results here, but I don't know if AngleComponent compensates for it.

XNA vectors are mutable, and also, why should it just swap one component of each vector?

You want to swap only one component because this is what your collision looks like: (excuse the atrocious diagramming)

In my code, I rotate v and u by the angle of y'. (y'=\vec{AB}) This means that the X components represent velocity parallel to the plane of collision, and the Y components represent velocity perpendicular to the plane of collision. Because the force imparted by the collision is only perpendicular to the plane of collision, I only swap the Y components, and leave the X components unaffected. Then rotate back by the same angle, so that the velocities are in terms of the absolute XY axes again.

[Tharwen]

It turns out I had been compensating for it in AngleComponent, so I just took out the extra 90 degrees in both functions since it seems fairly pointless to have it in there.

That explanation now makes sense, too. I was actually doing something slightly different; I was trying to find the portion of each object's momentum that was in the direction of the collision, then subtracting that from the object and adding it to the other one. Here's AngleComponent:

Code: Select all

        public Vector2 AngleComponent(Vector2 Vector, float Angle)
        {
            float newX = (float)(Vector.X * Math.Cos(Angle));
            float newY = (float)(Vector.Y * Math.Sin(Angle));

            return new Vector2(newX, newY);
        }


I hope that makes sense.

[Cosmologicon]

Not sure if you still can use this, but here's a previous post with my (IMHO complete) solution, in Javascript. It allows for non-elastic collisions as well, and has no trig calls and only one sqrt.

The sqrt is only used in the calculation that pushes the masses apart from each other. I don't see anything in your code that does this, but this (or some other measure) is necessary to prevent them from registering as colliding in the following frame when they're already moving apart.

[Tynach]

I'm not sure how relavent this topic still is, but I recently was creating an elastic collision system in Pygame. What I ended up with was this:

Code: Select all

def elastic(self, ball1, ball2):
   coll_norm = ball1.pos - ball2.pos
   distance = coll_norm.length()

   if distance == 0:
      print "Error: Balls occupy the same position in space. Avoiding Universal collapse."
      return

   coll_norm.normalize()

   coll_tan = Vector(coll_norm.y, -coll_norm.x)

   ball1norm = coll_norm * dot(ball1.speed, coll_norm)
   ball1tan = coll_tan * dot(ball1.speed, coll_tan)

   ball2norm = coll_norm * dot(ball2.speed, coll_norm)
   ball2tan = coll_tan * dot(ball2.speed, coll_tan)

   ball1.speed = ball1tan + coll_norm * ball2norm.length()
   ball2.speed = ball2tan - coll_norm * ball1norm.length()


for the elastic collisions, where 'Vector' is a class I made that is roughly equal to the XNA Vector2 class. That, of course, is only called if two balls are colliding, and it checks every ball against every other ball, with this code:

Code: Select all

for u, item1 in enumerate(self.items):
   for v, item2 in enumerate(self.items):
      if v <= u:
         continue
      elif item1.collide(item2):
         self.elastic(item1, item2)


So basically, it checks it kinda like this, assuming 3 balls:

ball1 vs ball1 -> Move to next pair (don't compute)
ball1 vs ball2 -> Check for collision
ball1 vs ball3 -> Check for collision

ball2 vs ball1 -> Move on
ball2 vs ball2 -> Move on
ball2 vs ball3 -> Check for collision

ball3 vs ball1 -> Move on
ball3 vs ball2 -> Move on
ball3 vs ball3 -> Move on

Not the most optimal solution, I know, but it thus only checks for colliding balls for each unique combination of balls - and then of course it calls elastic(ballx, bally) if there IS a collision - and doesn't check a ball against itself either. Note that 'items' are balls.

Also note that I treat collision with the 'walls' (window borders) entirely separate, calculating those when I move the balls themselves.

Attached I have put a .zip file with my current source code (after several revisions). I'm posting it here because I hope it may be helpful for other people looking for solutions, and I didn't see a rule for bumping old threads (which surprised me - and I did look at the rules, though I may have missed it).

Do note, I made my variables somewhat more... Descriptive in name. I also cut out some parts of the code that wikipedia has (like moving the balls apart as part of the collision code - I see no point in that, and for me, it just messed up my dirty rectangle updating. Or it did in my original draft, but I think it wouldn't with this version. Don't want to put it in though - my collisions look nice to me as they are). But it is based on the wikipedia code, and is essentially the same algorithm.

Edit: By the way, my collisions at one time were WAY wonky because I had these lines wrong:

Code: Select all

ball1.speed = ball1tan + coll_norm * ball2norm.length()
ball2.speed = ball2tan - coll_norm * ball1norm.length()


I had both of them as 'ballxtan + coll_norm', instead of having the second one subtract. You may want to look at that, and other things related to signs and so forth.