Horizons¶
import matplotlib.pyplot as plt
import numpy as np
import sxs
horizons = sxs.load("SXS:BBH:0123/Lev/Horizons.h5")
Skipping download from 'https://data.black-holes.org/catalog.json' because local file is newer Found the following files to load from the SXS catalog: SXS:BBH:0123v4/Lev5/Horizons.h5
The horizons
object has three main attributes, which just store quantities related to the three horizons:
horizons.A, horizons.B, horizons.C
(<sxs.horizons.HorizonQuantities at 0x15d2044c0>, <sxs.horizons.HorizonQuantities at 0x162ee2490>, <sxs.horizons.HorizonQuantities at 0x15d204610>)
Each horizon object stores information about its masses (areal and Christodoulou), spins (dimensionless and dimensionful), and coordinate trajectory.
We can list all the (public) attributes of each horizon:
[attribute for attribute in dir(horizons.A) if not attribute.startswith("_")]
['areal_mass', 'chi_inertial', 'chi_inertial_mag', 'chi_mag_inertial', 'christodoulou_mass', 'coord_center_inertial', 'dimensionful_inertial_spin', 'dimensionful_inertial_spin_mag', 'time']
[Note that chi_mag_inertial
is just an alias for chi_inertial_mag
, kept for backwards compatibility. The _mag
attributes just return the Euclidean norm of the corresponding vectors.]
For example, we can plot both types of masses for all three horizons:
plt.plot(horizons.A.time, horizons.A.christodoulou_mass, label=r"Christodoulou mass $A$")
plt.plot(horizons.A.time, horizons.A.areal_mass, label=r"Areal mass $A$")
plt.plot(horizons.B.time, horizons.B.christodoulou_mass, label=r"Christodoulou mass $B$")
plt.plot(horizons.B.time, horizons.B.areal_mass, label=r"Areal mass $B$")
plt.plot(horizons.C.time, horizons.C.christodoulou_mass, label=r"Christodoulou mass $C$")
plt.plot(horizons.C.time, horizons.C.areal_mass, label=r"Areal mass $C$")
plt.title("Horizon masses")
plt.ylabel(r"Mass (simulation units)")
plt.xlabel(r"Time (simulation units)")
plt.ylim(bottom=0.0)
plt.legend();
Or we can plot the coordinate trajectory of horizon A:
plt.plot(horizons.A.time, horizons.A.coord_center_inertial)
plt.title("Horizon A coordinate center in the inertial frame")
plt.ylabel(r"Position vector components (simulation units)")
plt.xlabel(r"Time (simulation units)")
plt.legend([r"$x$", r"$y$", r"$z$"]);
We can measure the Newtonian center of mass motion as a function of time:
com = horizons.newtonian_com
We will see below that the motion is distinctly non-inertial, so we also fit the measured motion to an averaged inertial motion:
# x_i and v_i are the best-fit initial translation and velocity
x_i, v_i, t_i, t_f = horizons.average_com_motion()
# Here, we express these as functions of time for plotting purposes
com_average = sxs.TimeSeries(
x_i[np.newaxis] + v_i[np.newaxis] * horizons.A.time[:, np.newaxis],
horizons.A.time
)
Now, we can plot them:
# Plot the measured and averaged components (and keep the lines)
measured = plt.plot(horizons.A.time, com)
averaged = plt.plot(horizons.A.time, com_average, ls='dashed')
# Copy colors from solid to dashed lines
for l1, l2 in zip(measured, averaged):
l2.set_color(l1.get_color())
plt.title("Newtonian center of mass of A-B horizon system")
plt.ylabel(r"Position vector components (simulation units)")
plt.xlabel(r"Time (simulation units)")
plt.legend([r"$x$", r"$y$", r"$z$"]);
There are some remarkable features of this plot:
- There are large oscillations that are not consistent with any motions we expect on physical grounds.
- The oscillations change suddenly just before merger — right around the time a gauge change is applied in the simulation. This suggests, of course, that the oscillations can be at least partially attributed to gauge weirdness.
- Beyond the oscillations there is also significant drift, ending up near merger displaced by around $0.1M$ from the origin of coordinates.
The averaged motion is, of course, just a Poincaré transformation, so we can remove it from the data. When we remove this transformation from the waveforms we find that certain unexpected features in the waveform disappear. This was the subject of this paper as well as a paper by CJ Woodford et al., and we will see more about this in the next notebook.
Backwards compatibility¶
Previously, the recommended way to load data from a file was to open it with h5py
, and then directly extract datasets. The user might do something like this:
with h5py.File("Horizons.h5", "r") as horizons:
time = horizons["AhA.dir/ArealMass.dat"][:, 0]
areal_mass = horizons["AhA.dir/ArealMass.dat"][:, 1]
With newer data, this code will break because the format of the files has changed. However, we can use the horizons
object to provide an interface that is as close as possible to the old interface so that code that works with the old files can continue to work with new files, with a minimal set of changes.
In particular, it is still possible to extract that data from the horizons
object:
horizons["AhA.dir/ArealMass.dat"]
array([[0.00000000e+00, 5.06401589e-01], [1.00000000e+00, 5.06444626e-01], [2.00000000e+00, 5.06550360e-01], ..., [4.36400000e+03, 5.08236829e-01], [4.36450000e+03, 5.08720136e-01], [4.36500000e+03, 5.09431577e-01]])
This actually reconstructs the arrays that the HDF5 file would contain on the fly. The same is true for all of the datasets found in the old-format Horizons.h5
files. (And we will see that it is also true for waveform objects.)
In fact, we even have a simple little context manager that takes the place of h5py.File
:
with sxs.loadcontext("Horizons.h5") as horizons:
time = horizons["AhA.dir/ArealMass.dat"][:, 0]
areal_mass = horizons["AhA.dir/ArealMass.dat"][:, 1]
Thus, only one function call would need to change to use newer files. (And in fact, the loadcontext
function is even a bit nicer because it can handle the downloading and caching for you, just like sxs.load
.)
However, be aware that this may not an be efficient use of memory, and is almost certainly slower than the newer interfaces. Wherever possible, you should update your code to use newer interfaces. Failing to do so will leave you open to ridicule from your peers and loved ones.
Continue with the introduction to waveforms.