# Minor Gripe

## Current state of the Elixir Rigid Physics library

After a few weeks of on-and-off work (maybe 2 weeks actual work punctuated by coder’s block and wonderful socializing here at RC), my timeline looks like this:

• Create stub physics engine server in Elixir as a GenServer(done)
• Create visualizer for engine state using Three.js and Phoenix Channels (done)
• Add first body type to physics engine: boxes! (done)
• Add framework for doing dynamics sans collisions in engine (in progress, more below)
• Add collision detection framework for finding contact manifolds
• Add more types of bodies beyond boxes
• Figure out how to spread the simulation across multiple processes
• Figure out how to spread the simulation across multiple nodes

I’m currently kinda blargh on this, since it’s a long road with a lot of work and not a lot of suprises. That said, it’s definetely a bit outside of my comfort zone and has been helping me improve my Elixir understanding (more intuition about processes, more knowledge about debugging, more familiarity with Phoenix Channels) and definitely brushing up on my physics.

## Assorted learnings on the project so far

Graphmath is incomplete in its quaternion support. The quaternion support in graphmath is solid-ish and has tests, but we’re still missing a lot of basic functionality that I think most people would expect. Things like creating a quaternion from an axis and an angle, making the identity quaternion, rotating an orientation quaternion by some amount, and so forth.

Writing eval functions in Elixir is surprisingly easy. It took maaaaaybe five lines of code to add server-side code evluation to the REPL, and that’s been my main workhorse for testing and interacting with the physics engine. Behold:

`````` def handle_in("server_eval", %{"body"=> body}, socket) do
try do
{ans, bindings} = Code.eval_string(body, [socket: socket ] ++ Map.get(socket.assigns, :bindings, []), __ENV__)
push(socket, "server_eval_result", %{body: inspect(ans)} )
rescue
error -> IO.inspect( error, label: "Error")
push(socket, "server_error_msg", %{body: inspect(error)} )
end
end``````

Erlang tools can do weird things to GenServers. To wit, starting `:debugger.start` or `:observer.start` seems to send a success message to the server…which, since it was missing a `handle_info`, would explode and cause me to lose my connection to the engine. This is obviously not great, so I added a catchall to prevent that.

Map comprehensions are a great way of dealing with the grumpiness of the Jason encoder. Like, sure, it’s obvious to you and me that Erlang tuples should be arrays in Javascript—but that’s apparently too great a conceptual leap for our tooling. So, we do something cute like this:

``````
defimpl Jason.Encoder, for: ElixirRigidPhysics.World do
def encode(%ElixirRigidPhysics.World{bodies: bodies} = world, opts) do
clean_bodies = for {ref, body} <- bodies, into: %{} do
clean_body = for {k, v} <- body, into: %{} do
case {k,v} do
{k, v} when is_reference(v) -> {k, "#{inspect v}"}
{k, v} when is_tuple(v)-> {k, Tuple.to_list(v)}
{k, v} -> {k,v}
end
end
{inspect(ref), clean_body}
end
Jason.Encode.map( %ElixirRigidPhysics.World{ world | bodies: clean_bodies}, opts)
end
end``````

And yes, one of my colleagues upon viewing that snippet starting whispering a reprisal of Let the bodies hit the floor. Naming variables is like a whole thing.

Records are really great if you don’t want maps. I haven’t profiled it yet (yes, yes, shame on me, Uncle Joe is turning over in his grave) but using records to store state for a GenServer seems to have worked out quite nicely at the expense of some minorly more annoying typing. Usage of that looks like:

``````  require Record
Record.defrecord(:sim, world: nil, subscribers: MapSet.new(), next_tick: nil)

### ...snip...

@impl true
def handle_info(:tick_simulation, s) do
new_world = ElixirRigidPhysics.Dynamics.step(sim(s, :world), @tick_rate / 1000)
thandle = Process.send_after(self(), :tick_simulation, @tick_rate)
update_subscribers(sim(s, :subscribers), new_world)
{:noreply, sim(s, world: new_world, next_tick: thandle)}
end``````

Implementing a GenServer should be done starting with casts, and upgrading to calls only as needed. Aforementioned colleague was unamused at the large amount of synchronous calls I’d been making, and pointed out that that was maybe not a great idea. Sure enough, while showing me how to use the debugger and attach to processes, a breakpoint caused the world process to become unresponsive.

That in and of itself isn’t a big deal, but the server process that minded it freaked out when its GenServer calls were timing out and promptly exploded itself in exasperation at being ignored. While I do empathize with this behavior, it makes debugging extremely difficult, and so we changed most of the methods to be casts to set the server’s expectations of attention to a (lower) more reasonable bar. GenServers are like online dating that way, I guess…if you send messages expecting a response, you’ll be disappointed and end up inconveniencing others.

HTML DOM text inputs don’t forward key press events for arrow keys. This was super annoying and took a lot of reading and experimentation to track down. I just wanted to add historical command navigation to my REPL, and for the life of me the blasted thing wouldn’t let me use arrow keys. I eventually figured out that I was listening to the wrong event.

## A quick ramble on dyanmics and rotation

I’ve been reading a lot of the docs and source code for other engines to crib ideas get inspiration. Cannon.js, which has a pretty nice and small codebase, has been very interesting and a good reference. Unfortunately, they did something I’m completely baffled by.

Say you’ve got your body you want to simulate. It’s a box, perhaps somewhat fittingly given my old line of work.

Anyays, the representation of it in your world (or, any world, really) is going to be a position and an orientation, better known in robotics literature as its pose. That’s a neat factoid, but not relevant to what we’re talking about now.

The position is where in the world frame (root coordinate system) it can be found. The orientation is how it is rotated compared to the default alignemnt of the frame—is it 45 degrees to the north, is it resting on its side, etc. Because reasons, I’m representing the orientation of a body as a quaternion.

Alright, neat. Our body has a pose. The dynamics bit, though, requires us to be able to move the body subject to some force. That in turn means applying some impulse to change the momentum (in our case, this really means velocity) of our body, and that in turn means integrating the velocity with the original pose to create our new pose.

So, for the position component, this isn’t too bad. We add a linear velocity field to our body, and then we can use something like forward Euler whereby we scale the velocity by the timestemp and add it to the position. Very fast, very simple, unstable over repeated iterations, but definetely Good Enough(tm) for our purposes.

The part that’s been annoying to me is that angular velocity bit.

Same setup: we have our orientation quaternion, and then we need some way of storing the current angular velocity. That’s our first hitch, though…how do we store that? Quaternion? Matrix? Rotor? What even is the derivative of a quaternion?

Well, I did what any good engineer would do and I looked at prior art. Cannon stores its angular velocity as a 3D vector. Okay, that’s weird, why not a quaternion, I thought to myself, since quaternions are for spinny things and vectors are for slidey things.

Well, a comment in their code suggests that they’re storing it as an axis of rotation, and then the magnitude of the vector is the rate of rotation about that axis. I’m personally unconvinced that this will work out in anything except tears, but I’ve seen the code run, so I have faith.

You might wonder why I’m skeptical. See, somewhere down the line in this implementation we’re going to have to account for things changeing that angular velocity on the way to changing the angular momentum of the body. Imagine a box getting hit on two adjacent sides at once, right? We expect that the sum of the changes in angular momentum produce a resultant angular velocity that is some scaled sum. So, in your head, you should be picturing a function that takes two axis-angle representations of impulse and the moment of intertia of our box and yields a new angular velocity vector. Something about that makes my brain hurt.

Anyways, back to the present, we need to—as my freshman geometry teacher would say—turn ducks into concrete.

We have an orientation quaternion, a timestep, and an angular velocity axis-angle thing and somehow need to moosh them to create a new orientation quaternion.

I assume that the way this would work is that we scale the angular velocity vector by the timestep to get how much change we expect in the current simulation frame (please God let this just be a scaling of the angle component but life is suffering and we must all eat at Arby’s). That done, we need to transform an orientation quaternion by the result.

Now, I think that that might be as simple as converting the scaled angular velocity vector to a quaternion itself and multiplying them together, but I’m unsure.

For reference, the code I’m trying to understand:

``````/**
* Rotate an absolute orientation quaternion given an angular velocity and a time step.
* @param  {Vec3} angularVelocity
* @param  {number} dt
* @param  {Vec3} angularFactor
* @param  {Quaternion} target
* @return {Quaternion} The "target" object
*/
Quaternion.prototype.integrate = function(angularVelocity, dt, angularFactor, target){
target = target || new Quaternion();

var ax = angularVelocity.x * angularFactor.x,
ay = angularVelocity.y * angularFactor.y,
az = angularVelocity.z * angularFactor.z,
bx = this.x,
by = this.y,
bz = this.z,
bw = this.w;

var half_dt = dt * 0.5;

target.x += half_dt * (ax * bw + ay * bz - az * by);
target.y += half_dt * (ay * bw + az * bx - ax * bz);
target.z += half_dt * (az * bw + ax * by - ay * bx);
target.w += half_dt * (- ax * bx - ay * by - az * bz);

return target;
};``````

Some things leap out at me:

• `half_dt` reminds me of the use of half-angles used in some descriptions of rotations. Unsure if that’s a red herring or not.
• The bits about rotating the target code look suspiciously like a quat-quat multiplication, but we’re missing a term!
• There’s no `aw` term, but maybe that’s because that element of the quaternion would be zero. In fact, further examination convinces me of this.
• Okay, why is the `aw` term zero? Closest I’ve got is that `<0,v>` is a valid way of converting an axis-angle vector `v` to a rotation quaternion. Notes say that the resulting quaternion needn’t be unit length, but why?

## Conclusion

Onwards ever onwards. I don’t really understand all the math at play, and that is very frustrating. That said, I think I’ll be able to get this to work.