A month ago, as part of my first foray into LWJGL, I decided to make an OrbCollider: a small physics simulator where you take a group of "Orbs" (2d circles, though I plan to expand the physics into 3 dimensions soon) and apply "Elastic Rigid Body Collision Physics" to them. So I actually know quite a bit about programming a (rudimentary, mostly naive) physics engine.

The first thing you need to do is think about how you want your physics engine to work. Here were my considerations:

**Which is more important: [Computational] Speed or [Physical] Accuracy?** I ended up going with Accuracy for my implementation, but if you were doing collisions for a game you could very reasonably choose speed instead. It's all a matter of what your needs are.**What kind of physics are you modelling?** As I mentioned, I chose to model Elastic Rigid Body Collisions using disks. The only thing I don't simulate is rotation, either by friction with other Orbs or, if I had used non-circular objects, by glancing collisions. You need to research how these physics work and how to calculate their effects.**What is the logic flow of your Simulation?**

That last one will require expansion. For the remainder of this post, I'll be going into substantial detail on the implementation I created, and some of the rhetoric behind the decisions I made.

In my implementation, I worked out that the logic flow should work like this:

- First, we are given an "Amount of Time to Simulate, in Seconds" (dT). In my implementation, this usually comes in the form of fractions of a second (most runs of the simulate(dT) function take an input of ~.01). The actual units used don't matter so long as it's consistent.
- We then query every single Orb, against every single other Orb, and determine when the "soonest" collision occurs. In other words, if you're given 20 orbs, then you find out which two orbs collide first. This also includes the boundaries of the simulation; for all intents and purposes, a collision with a boundary counts as a collision with an orb. The physics is just different. the amount of time until the soonest collision we'll call mT.

*It's worth noting that this step is, by far, the most computationally expensive part of my program. I've been looking into an OpenCL solution to offload some of the work onto the GPU.*

- Then, we do a check against the time we're supposed to simulate:

**dT <= mT**: The amount of time we're simulating is less than the time it takes to reach the nearest collision. In this case, we advance ALL Orbs by dT seconds (they have velocity and direction vectors, so it's easy to calculate their new positions) and we're done simulating.**dT > mT && mT >= 0**: A collision will happen before we finish simulating. In this case, we advance ALL Orbs up to the point that a collision occurs (mT seconds). Then, we apply the collision algorithm to the two orbs that collided. Then, we subtract mT from dT and return to step 2. We will repeat this until dT == 0 or the first condition becomes true.

So the only things left unexplained are the physics calculations themselves. You'll almost certainly need to derive your own formulas, but I can give you the formulas I used (which I specifically developed for my implementation), which I give you permission to copy/modify however you wish.

To start with: we need to know how Orbs and Barriers are defined. In my program, Orbs have the following properties, all of which are Doubles:

**Position**: 2 variables: x and y**Velocity**: v**Angle**: a**Mass**: m**Radius**: r**... Other Variables which aren't relevant to this discussion** (Color, texture, etc)

Barriers:

**Position**: 2 variables: x and y**Angle**: a

In my implementation, Barriers are defined as Infinite Lines that pass through the specified Position. A smarter implementation would give the lines length and only check collisions within that length, which would let me build Finite Obstacles; I haven't done that yet. The Angle is the Angle that it "repels" orbs, which means that the "line" describing the barrier will be perpendicular to the Angle. This lets me build a "Pen" to keep the orbs contained by defining 3 or more barriers and pointing them inwards.

*Implementation note: Unless specified otherwise, angles are in Radians. Also, this is where things get super Math-heavy, so if you can't follow what's happening here, a Physics engine probably won't be for you.*So, in my implementation, I need to first know, given two Orbs (o

_{1}, o

_{2}), when do they collide? For the future, any references to Orb properties will be subscripted, i.e. the mass of the second orb will be m

_{2}, the angle of the first orb will be a

_{1}.

When do two Orbs Collide? When the distance between their centers (sqrt((x

_{1}-x

_{2})

^{2} + (y

_{1}-y

_{2})

^{2})) is less than the sum of their respective Radii (r

_{1}+r

_{2}). Measured relative to time (t), the formula we have to solve for t is:

r

_{1}+r

_{2} = sqrt((x

_{1}(t)-x

_{2}(t))

^{2} + (y

_{1}(t)-y

_{2}(t))

^{2})

So, how are x

_{n}(t) or y

_{n}(t) defined? Well, that's actually pretty simple:

x

_{n}(t) = x

_{n} + v

_{n}*t*cos(a

_{n})

y

_{n}(t) is calculated exactly the same, except with sin instead of cos. So, filling in the formulas, you get the next equation:

r

_{1}+r

_{2} = sqrt((x

_{1}+t*v

_{1}*cos(a

_{1}) - x

_{2}-t*v

_{2}*cos(a

_{2}))

^{2} + (y

_{1}+t*v

_{1}*sin(a

_{1}) - y

_{2}-t*v

_{2}*sin(a

_{2}))

^{2})

We can simplify this a little by combining similar terms: dx = x

_{1}-x

_{2}, dvx = vx

_{1}*[(t*v*_{1}*cos(a_{1}))] - vx

_{2}, dy = y

_{1}-y

_{2}, dvy = vy

_{1}-vy

_{2}. The formula then looks a little more manageable:

r

_{1}+r

_{2} = sqrt((dx + t*dvx)

^{2} + (dy + t*dvy)

^{2})

Solving for t, we eventually have to apply the quadratic formula (I'm not showing the intermediate steps) and the formula will look like this:

a = dvx

^{2} + dvy

^{2}b = 2*(dx*dvx + dy*dvy)

c = dx

^{2} + dy

^{2} - (r

_{1}+r

_{2})

^{2}Thus, finally, t equals the following two values:

t = (-b + sqrt(b

^{2}-4*a*c))/(2*a)

t = (-b - sqrt(b

^{2}-4*a*c))/(2*a)

*Implementation Note: b*^{2}-4*a*c can occasionally have a negative value! This **needs** to be checked for at runtime or you will have arithmetic exceptions! Fortunately, this is a useful case: if b^{2}-4*a*c is negative, then these two Orbs will NEVER collide, past present or future, and you can tell the program that the time to collide is Positive Infinity (or some other suitably large number)Given these two values of t, we then make the following observations:

**If both values of t are positive**, the smaller value is the time until they collide.**If both values of t are negative**, the collision will never happen, because the orbs are travelling away from each other (and theoretically collided "in the past"). You can return Positive Infinity for this situation.**If one value is positive and the other is negative**, then the orbs are currently intersecting. Mathematically, this should never happen. Realistically, due to rounding errors and other weirdness, this can happen. Return 0 as the time until they collide.

One optimization I made is to calculate the derivative of their distance, which involves a much simpler calculation: dvy*dy+dvx*dx. If this value is positive, the orbs are moving away from each other (their distance is increasing) and we assume they never collide.

*Implementation Note: Yeah, that's not the complete calculation of the derivative of their distance. But the missing component is guaranteed to be a non-negative number, which means it will never affect the sign of the calculation (which is all we need).*So now you can determine the precise moment when two Orbs will collide. But what then? How do they collide?

In one dimension, given two bodies with mass and velocity, this is what their final velocities look like:

vf

_{1} = (vi

_{1}*(m

_{1}-m

_{2}) + 2*m

_{2}*vi

_{2}) / (m

_{1}+m

_{2})

vf

_{2} = (vi

_{2}*(m

_{2}-m

_{1}) + 2*m

_{1}*vi

_{1}) / (m

_{1}+m

_{2})

In two dimensions, you need to calculate the relative velocities of the two Orbs normal to the angle of impact. This can be difficult to visualize, but imagine rotating the system until one orb is located exactly PI/2 radians above the other. Then redraw the motion vectors relative to this rotated system. Now, the only thing you need to care about is the vertical (y) velocities of the Orbs, since in this system, the x Velocity has no effect on the collision.

So, the math: first find the point of Impact. This is helpfully given by the two equations:

px = (x

_{1}*r

_{2} + x

_{2}*r

_{1}) / (r

_{1}+r

_{2}),

py = (y

_{1}*r

_{2} + y

_{2}*r

_{1}) / (r

_{1}+r

_{2})

*Implementation Note: We are assuming that the orbs are perfectly touching each other. In my implementation, this is a reasonable assumption. However, due to rounding errors, it may not be true that the orbs are /exactly/ touching, either by virtue of intersecting slightly or being slightly separated. This is a floating point error, not a math error, and it doesn't really matter anyways.*So now the angle of impact: as I mentioned, we want one orb to sit just above the other Orb, with a perfectly vertical line that intersects the center of both Orbs. So once we calculate the angle of Impact (aI) we adjust it by adding PI/2 to it.

aI = atan2(py-y

_{1}, px-x

_{1}) + PI/2

*Implementation Note: java.lang.Math.atan2(double,double) returns angles in the full range of PI to -PI, unlike the mathematical arctangent function which only returns values within PI/2 to -PI/2 due to its inability to distinguish whether the y component or the x component was the negative value.*The "adjusted angles" of the Orbs, then, become

aj

_{1} = a

_{1} - aI

aj

_{2} = a

_{2} - aI

Now, given these adjusted angles, we calculate vi

_{1} to be equal to vy

_{1} (which is v

_{1}*sin(aj

_{1})) and vi

_{2} to be equal to vy

_{2} and plug them in to the one-dimensional formula I described above. Remember that at this point, because the x velocity is tangent to the point of collision, it has no effect on the collision (you could imagine this as a scenario where two orbs barely graze each other as they fly past each other; however fast they are moving, they don't affect the other sans friction and air resistance).

Once we have the new y velocities, we calculate the final velocity as vf

_{n} = sqrt(vfx

_{n}^{2}+vfy

_{n}^{2}), and the final angle a

_{n} = atan2(vfy

_{n}, vfx

_{n}) + aI. Remember that in those two formulas, vfx

_{n} is equal to vix

_{n}These are the two calculations you need to build a functional Physics Engine.

I leave it an an exercise to the reader to work out how collisions with Boundaries work. The collisions themselves are absurdly simple, merely an adjustment to the angle. Calculating the time until a collision is more complicated, though I'll give this hint: find the point on the orb that is nearest to the line, and calculate how long it takes for that point, given the Orb's velocity and angle, to collide with the line.

One more thing: to the best of my knowledge, the information I've presented here is correct. I know for certain that the algorithms I implemented in my program are correct, but I may have made a typo here or there in the process of transcribing my algorithms. Also, there are other considerations I haven't mentioned here: for example, if (through rounding errors) two orbs are intersecting and they collide, future calculations may still indicate that the orbs are intersecting, and therefore due to collide in 0 seconds, resulting in an infinite loop. There may be other issues that you need a solution for. All of these are things that will be implementation specific, and I can't offer general advice on fixing those issues.

Good luck, and let me know if you have any follow-up questions.