The masses/inertias of each of those cartesian coordinates (m1 for x1
and y1, m2 for x2 and y2)

A potential energy function for your objects:

U = (m1 y1 + m2 y2) * g

And that’s it! Hamiltonian mechanics steps your generalized coordinates (θ1
and θ2) through time, without needing to do any simulation involving
x1/y1/x2/y2! And you don’t need to worry about tension or any other
stuff like that. All you need is a description of your coordinate system
itself, and the potential energy!

Thanks to ~~Alexander~~ William Rowan Hamilton, we can express our
system parameterized by arbitrary coordinates and get back equations of motions
as first-order differential equations. This library solves those first-order
differential equations for you using automatic differentiation and some matrix
manipulation.

Configurations are nice, but Hamiltonian dynamics is all about motion through
phase space, so let’s convert this configuration-space representation of the
state into a phase-space representation of the state:

phase0 :: Phase 2
phase0 = toPhase doublePendulum config0

And now we can ask for the state of our system at any amount of points in time!

Potential improvements

Time-dependent systems: Shouldn’t be an problem in theory/math; just
add a time parameter before all of the functions. This opens a lot of
doors, like deriving inertial forces for free (like the famous Coriolis
force and centrifugal force).

The only thing is that it makes the API pretty inconvenient, because it’d
require all of the functions to also take a time parameter. Of course, the
easy way out/ugly solution would be to just offer two versions of the same
function (one for time-independent systems and one for time-dependent
systems. But this is un-ideal.

Velocity-dependent potentials: Would give us the ability to model systems
with velocity-dependent Lagrangians like a charged particle in an
electromagnetic field, and also dissipative systems, like systems with
friction (dependent on signum v) and linear & quadratic wind resistance.

This issue is much harder, theoretically. It involves inverting arbitrary
functions forall a. RealFloat a => V.Vector n a -> V.Vector m a. It
might be possible with the help of some
bidirectionalization techniques, but I can’t get the bff
package to compile, and I’m not sure how to get bff-mono to work with
numeric functions.

If anyone is familiar with bidirectionalization techniques and is willing
to help out, please send me a message or open an issue! :)