How to build a unit system where G = 1

How to build a unit system where G = 1#

Gravitational dynamics codes almost always work in units where the gravitational constant is exactly 1. It removes \(G\) from every force evaluation in the inner loop, and it keeps positions, masses and times at order unity instead of spanning forty orders of magnitude in SI.

unxt builds such a system for you: fix any two of length, mass and time, and it solves for the third so that \(G = 1\).

>>> import numpy as np
>>> import unxt as u
>>> from astropy import constants as const
>>> from unxt.unitsystems import DynamicalSimUSysFlag

Fix length and mass, solve for time#

Pass DynamicalSimUSysFlag as the first argument, then the two units you want to keep. For galactic dynamics, kiloparsecs and solar masses:

>>> usys = u.unitsystem(DynamicalSimUSysFlag, "kpc", "solMass")
>>> usys
LengthMassTimeUnitSystem(length=Unit("kpc"), mass=Unit("solMass"), time=Unit("...kpc(3/2) s kg(1/2) / (solMass(1/2) m(3/2))"))

The time unit is whatever it had to be. It is not a round number and it is not meant to be read — it is the time unit in which \(G\) comes out as 1:

>>> bool(np.isclose(const.G.decompose(usys).value, 1.0))
True

To see what that time unit actually is in seconds, decompose it:

>>> usys["time"].decompose()
Unit("...e+19 s")

or convert one of them to something familiar:

>>> u.Q(1.0, usys["time"]).uconvert("Myr")
Quantity(Array(471483., dtype=float32...), unit='Myr')

Fix length and time instead#

If your integrator’s timestep is the thing you care about, fix time and let the mass unit fall out:

>>> usys2 = u.unitsystem(DynamicalSimUSysFlag, "kpc", "Myr")
>>> usys2["mass"].decompose()
Unit("...e+41 kg")

>>> bool(np.isclose(const.G.decompose(usys2).value, 1.0))
True

Same guarantee, different free choice. Which two you fix is a question of which numbers you want to be readable.

Compare with a plain unit system#

The built-in galactic system uses the same kind of units but makes no promise about \(G\):

>>> float(const.G.decompose(u.unitsystem("galactic")).value)
4.4985021514695...e-12

That factor is exactly what a DynamicalSimUSysFlag system saves you from carrying through every force calculation.

Move values in and out#

Convert into the system with ustrip, passing the whole system so it picks the unit matching each quantity’s dimension:

>>> float(u.Q(220.0, "km/s").ustrip(usys))
106082.1...

A circular velocity of 220 km/s is about 1.06e5 in these units. Going the other way, build the quantity with the system’s own unit and convert out:

>>> u.Q(1.0, usys["time"]).uconvert("Myr")
Quantity(Array(471483., dtype=float32...), unit='Myr')

Derived units come out already reduced to astropy’s own base units, since the units they are built from are scaled composites with no name of their own:

>>> usys["velocity"]
Unit("2.07387 m / s")

>>> usys["acceleration"]
Unit("1.39383e-19 m / s2")

One velocity unit in this system is 2.07 m/s — which is the conversion ustrip(usys) just did for you.