From d38ab7fdd3174869380a0b5fd03cca938e4bf772 Mon Sep 17 00:00:00 2001
From: github-actions <github-actions@github.com>
Date: Wed, 18 Jan 2023 21:20:55 +0000
Subject: [PATCH 01/13] Bump version to v2022.8.4

---
 CITATION.cff   | 2 +-
 pyproject.toml | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/CITATION.cff b/CITATION.cff
index cd28be30..cddbf0c1 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -26,7 +26,7 @@ authors:
   given-names: "Keefe"
   orcid: "https://orcid.org/0000-0003-0276-3856"
 title: "scri"
-version: 2022.8.3
+version: 2022.8.4
 doi: 10.5281/zenodo.4041971
 date-released: 2022-10-18
 url: "https://github.com/moble/scri"
diff --git a/pyproject.toml b/pyproject.toml
index d5d68cad..535b2e69 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "scri"
-version = "2022.8.3"
+version = "2022.8.4"
 description = "Time-dependent functions of spin-weighted spherical harmonics"
 readme = "README.md"
 license = "MIT"

From 336576fabbbd183e44194ad8e8b607f4348587db Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Fri, 10 Feb 2023 12:59:42 -0800
Subject: [PATCH 02/13] Update superrest frame code to use BMSTransformations
 class and add option to map abd objects to the frame of another abd object
 via map_to_abd_frame.py

---
 scri/asymptotic_bondi_data/__init__.py        |   3 +
 scri/asymptotic_bondi_data/bms_charges.py     |  45 +-
 .../asymptotic_bondi_data/map_to_abd_frame.py | 431 ++++++++++++
 .../map_to_superrest_frame.py                 | 664 ++++++++++--------
 tests/test_abd_frame.py                       |  91 +++
 tests/test_superrest_frame.py                 |  79 ++-
 6 files changed, 996 insertions(+), 317 deletions(-)
 create mode 100644 scri/asymptotic_bondi_data/map_to_abd_frame.py
 create mode 100644 tests/test_abd_frame.py

diff --git a/scri/asymptotic_bondi_data/__init__.py b/scri/asymptotic_bondi_data/__init__.py
index c762511b..ac3aaa39 100644
--- a/scri/asymptotic_bondi_data/__init__.py
+++ b/scri/asymptotic_bondi_data/__init__.py
@@ -180,6 +180,7 @@ def interpolate(self, new_times):
         # interpolate frame data if necessary
         if self.frame.shape[0] == self.n_times:
             import quaternion
+
             new_abd.frame = quaternion.squad(self.frame, self.t, new_times)
         return new_abd
 
@@ -203,6 +204,7 @@ def interpolate(self, new_times):
         bondi_rest_mass,
         bondi_four_momentum,
         bondi_angular_momentum,
+        CWWY_angular_momentum,
         bondi_dimensionless_spin,
         bondi_boost_charge,
         bondi_CoM_charge,
@@ -210,3 +212,4 @@ def interpolate(self, new_times):
     )
 
     from .map_to_superrest_frame import map_to_superrest_frame
+    from .map_to_abd_frame import map_to_abd_frame
diff --git a/scri/asymptotic_bondi_data/bms_charges.py b/scri/asymptotic_bondi_data/bms_charges.py
index 65edbb67..b8c01c27 100644
--- a/scri/asymptotic_bondi_data/bms_charges.py
+++ b/scri/asymptotic_bondi_data/bms_charges.py
@@ -5,8 +5,10 @@
 ### class.  In particular, they assume that the first argument, `self` is an instance of
 ### AsymptoticBondData.  They should probably not be used outside of that class.
 
+import scri
 import numpy as np
 from math import sqrt, pi
+import spherical_functions as sf
 
 
 def mass_aspect(self, truncate_ell=max):
@@ -83,7 +85,7 @@ def bondi_four_momentum(self):
     return charge_vector_from_aspect(charge_aspect)
 
 
-def bondi_angular_momentum(self, output_dimensionless=False):
+def bondi_angular_momentum(self):
     """Compute the total Bondi angular momentum vector
 
         i (ψ₁ + σ ðσ̄)
@@ -99,6 +101,37 @@ def bondi_angular_momentum(self, output_dimensionless=False):
     return charge_vector_from_aspect(charge_aspect)[:, 1:]
 
 
+def CWWY_angular_momentum(self):
+    """Compute the Chen/Wang/Wang/Yau angular momentum vector.
+
+    See Eq. (5) of https://arxiv.org/abs/2102.03235
+
+    """
+    ell_max = 1  # Compute only the parts we need, ell<=1
+    supertranslation_potential_data = scri.asymptotic_bondi_data.map_to_bms_frame.𝔇inverse(
+        np.array(self.sigma.ethbar_GHP.ethbar_GHP + self.sigma.bar.eth_GHP.eth_GHP), self.ell_max
+    )
+    supertranslation_potential = scri.ModesTimeSeries(
+        sf.SWSH_modes.Modes(
+            supertranslation_potential_data,
+            spin_weight=0,
+            ell_min=0,
+            ell_max=self.ell_max,
+            multiplication_truncator=max,
+        ),
+        time=self.t,
+    )
+    charge_aspect = (
+        1j
+        * (
+            self.psi1.truncate_ell(ell_max)
+            + self.sigma.multiply(self.sigma.bar.eth_GHP, truncator=lambda tup: ell_max)
+            + supertranslation_potential.multiply(self.mass_aspect().eth_GHP, truncator=lambda tup: ell_max)
+        )
+    ).view(np.ndarray)
+    return charge_vector_from_aspect(charge_aspect)[:, 1:]
+
+
 def bondi_dimensionless_spin(self):
     """Compute the dimensionless Bondi spin vector"""
     N = self.bondi_boost_charge()
@@ -112,7 +145,7 @@ def bondi_dimensionless_spin(self):
     vhat = v.copy()
     t_idx = v_norm != 0  # Get the indices for timesteps with non-zero velocity
     vhat[t_idx] = v[t_idx] / v_norm[t_idx, np.newaxis]
-    gamma = (1 / np.sqrt(1 - v_norm ** 2))[:, np.newaxis]
+    gamma = (1 / np.sqrt(1 - v_norm**2))[:, np.newaxis]
     J_dot_vhat = np.einsum("ij,ij->i", J, vhat)[:, np.newaxis]
     spin_charge = (gamma * (J + np.cross(v, N)) - (gamma - 1) * J_dot_vhat * vhat) / M_sqr
     return spin_charge
@@ -209,7 +242,9 @@ def supermomentum(self, supermomentum_def, **kwargs):
     if supermomentum_def.lower() in ["bondi-sachs", "bs"]:
         supermomentum = self.psi2 + self.sigma.grid_multiply(self.sigma.bar.dot, **kwargs)
     elif supermomentum_def.lower() in ["moreschi", "m"]:
-        supermomentum = self.psi2 + self.sigma.grid_multiply(self.sigma.bar.dot, **kwargs) + self.sigma.bar.eth_GHP.eth_GHP
+        supermomentum = (
+            self.psi2 + self.sigma.grid_multiply(self.sigma.bar.dot, **kwargs) + self.sigma.bar.eth_GHP.eth_GHP
+        )
     elif supermomentum_def.lower() in ["geroch", "g"]:
         supermomentum = (
             self.psi2
@@ -217,7 +252,9 @@ def supermomentum(self, supermomentum_def, **kwargs):
             + 0.5 * (self.sigma.bar.eth_GHP.eth_GHP - self.sigma.ethbar_GHP.ethbar_GHP)
         )
     elif supermomentum_def.lower() in ["geroch-winicour", "gw"]:
-        supermomentum = self.psi2 + self.sigma.grid_multiply(self.sigma.bar.dot, **kwargs) - self.sigma.ethbar_GHP.ethbar_GHP
+        supermomentum = (
+            self.psi2 + self.sigma.grid_multiply(self.sigma.bar.dot, **kwargs) - self.sigma.ethbar_GHP.ethbar_GHP
+        )
     else:
         raise ValueError(
             f"Supermomentum defintion '{supermomentum_def}' not recognized. Please choose one of "
diff --git a/scri/asymptotic_bondi_data/map_to_abd_frame.py b/scri/asymptotic_bondi_data/map_to_abd_frame.py
new file mode 100644
index 00000000..bec484f9
--- /dev/null
+++ b/scri/asymptotic_bondi_data/map_to_abd_frame.py
@@ -0,0 +1,431 @@
+import numpy as np
+
+import sxs
+import scri
+import spinsfast
+import spherical_functions as sf
+from spherical_functions import LM_index
+
+import quaternion
+from quaternion.calculus import derivative
+from quaternion.calculus import indefinite_integral as integrate
+
+from scipy.interpolate import CubicSpline
+
+from functools import partial
+
+from scri.asymptotic_bondi_data.map_to_superrest_frame import time_translation, rotation
+
+
+def cost(δt_δϕ, args):
+    modes_A, modes_B, t_reference, δϕ_factor, δΨ_factor, normalization = args
+
+    # Take the sqrt because least_squares squares the inputs...
+    diff = integrate(
+        np.sum(
+            abs(modes_A(t_reference + δt_δϕ[0]) * np.exp(1j * δt_δϕ[1]) ** δϕ_factor * δΨ_factor - modes_B) ** 2, axis=1
+        ),
+        t_reference,
+    )[-1]
+    return np.sqrt(diff / normalization)
+
+
+def align2d(h_A, h_B, t1, t2, n_brute_force_δt=None, n_brute_force_δϕ=None, include_modes=None, nprocs=4):
+    """Align waveforms by shifting in time and phase
+
+    This function determines the optimal time and phase offset to apply to `h_A` by minimizing
+    the averaged (over time) L² norm (over the sphere) of the difference of the h_Aveforms.
+
+    The integral is taken from time `t1` to `t2`.
+
+    Note that the input waveforms are assumed to be initially aligned at least well
+    enough that:
+
+      1) the time span from `t1` to `t2` in the two waveforms will overlap at
+         least slightly after the second waveform is shifted in time; and
+      2) waveform `h_B` contains all the times corresponding to `t1` to `t2`
+         in waveform `h_A`.
+
+    The first of these can usually be assured by simply aligning the peaks prior to
+    calling this function:
+
+        h_A.t -= h_A.max_norm_time() - h_B.max_norm_time()
+
+    The second assumption will be satisfied as long as `t1` is not too close to the
+    beginning of `h_B` and `t2` is not too close to the end.
+
+    Parameters
+    ----------
+    h_A : WaveformModes
+    h_B : WaveformModes
+        Waveforms to be aligned
+    t1 : float
+    t2 : float
+        Beginning and end of integration interval
+    n_brute_force_δt : int, optional
+        Number of evenly spaced δt values between (t1-t2) and (t2-t1) to sample
+        for the initial guess.  By default, this is just the maximum number of
+        time steps in the range (t1, t2) in the input waveforms.  If this is
+        too small, an incorrect local minimum may be found.
+    n_brute_force_δϕ : int, optional
+        Number of evenly spaced δϕ values between 0 and 2π to sample
+        for the initial guess.  By default, this is 2 * ell_max + 1.
+    include_modes: list, optional
+        A list containing the (ell, m) modes to be included in the L² norm.
+    nprocs: int, optional
+        Number of cpus to use.
+        Default is 4. 'None' corresponds to the maximum number. '-1' corresponds to no parallelization.
+
+    Returns
+    -------
+    optimum: OptimizeResult
+        Result of scipy.optimize.least_squares
+    h_A_prime: WaveformModes
+        Resulting waveform after transforming `h_A` using `optimum`
+
+    Notes
+    -----
+    Choosing the time interval is usually the most difficult choice to make when
+    aligning waveforms.  Assuming you want to align during inspiral, the times
+    must span sufficiently long that the waveforms' norm (equivalently, orbital
+    frequency changes) significantly from `t1` to `t2`.  This means that you
+    cannot always rely on a specific number of orbits, for example.  Also note
+    that neither number should be too close to the beginning or end of either
+    waveform, to provide some "wiggle room".
+
+    Precession generally causes no problems for this function.  In principle,
+    eccentricity, center-of-mass offsets, boosts, or other supertranslations could
+    cause problems, but this function begins with a brute-force method of finding
+    the optimal time offset that will avoid local minima in all but truly
+    outrageous situations.  In particular, as long as `t1` and `t2` are separated
+    by enough, there should never be a problem.
+
+    """
+    from scipy.optimize import least_squares
+
+    import multiprocessing as mp
+
+    h_A_copy = h_A.copy()
+    h_B_copy = h_B.copy()
+
+    # Check that (t1, t2) makes sense and is actually contained in both waveforms
+    if t2 <= t1:
+        raise ValueError(f"(t1,t2)=({t1}, {t2}) is out of order")
+    if h_A_copy.t[0] > t1 or h_A_copy.t[-1] < t2:
+        raise ValueError(
+            f"(t1,t2)=({t1}, {t2}) not contained in h_A_copy.t, which spans ({h_A_copy.t[0]}, {h_A_copy.t[-1]})"
+        )
+    if h_B_copy.t[0] > t1 or h_B_copy.t[-1] < t2:
+        raise ValueError(
+            f"(t1,t2)=({t1}, {t2}) not contained in h_B_copy.t, which spans ({h_B_copy.t[0]}, {h_B_copy.t[-1]})"
+        )
+
+    # Figure out time offsets to try
+    δt_lower = max(t1 - t2, h_A_copy.t[0] - t1)
+    δt_upper = min(t2 - t1, h_A_copy.t[-1] - t2)
+
+    # We'll start by brute forcing, sampling time offsets evenly at as many
+    # points as there are time steps in (t1,t2) in the input waveforms
+    if n_brute_force_δt is None:
+        n_brute_force_δt = max(
+            sum((h_A_copy.t >= t1) & (h_A_copy.t <= t2)), sum((h_B_copy.t >= t1) & (h_B_copy.t <= t2))
+        )
+    δt_brute_force = np.linspace(δt_lower, δt_upper, num=n_brute_force_δt)
+
+    if n_brute_force_δϕ is None:
+        n_brute_force_δϕ = 2 * h_A_copy.ell_max + 1
+    δϕ_brute_force = np.linspace(0, 2 * np.pi, n_brute_force_δϕ, endpoint=False)
+
+    δt_δϕ_brute_force = np.array(np.meshgrid(δt_brute_force, δϕ_brute_force)).T.reshape(-1, 2)
+
+    t_reference = h_B_copy.t[np.argmin(abs(h_B_copy.t - t1)) : np.argmin(abs(h_B_copy.t - t2)) + 1]
+
+    # Remove certain modes, if requested
+    ell_max = min(h_A_copy.ell_max, h_B_copy.ell_max)
+    if include_modes != None:
+        for L in range(2, ell_max + 1):
+            for M in range(-L, L + 1):
+                if not (L, M) in include_modes:
+                    h_A_copy.data[:, LM_index(L, M, h_A_copy.ell_min)] *= 0
+                    h_B_copy.data[:, LM_index(L, M, h_B_copy.ell_min)] *= 0
+
+    # Define the cost function
+    modes_A = CubicSpline(h_A_copy.t, h_A_copy[:, 2 : ell_max + 1].data)
+    modes_B = CubicSpline(h_B_copy.t, h_B_copy[:, 2 : ell_max + 1].data)(t_reference)
+
+    normalization = integrate(CubicSpline(h_B_copy.t, h_B_copy[:, 2 : ell_max + 1].norm())(t_reference), t_reference)[
+        -1
+    ]
+    if normalization == 0:
+        normalization = t_reference[-1] - t_reference[0]
+
+    δϕ_factor = np.array([M for L in range(h_A_copy.ell_min, ell_max + 1) for M in range(-L, L + 1)])
+
+    optimums = []
+    h_A_primes = []
+    for δΨ_factor in [-1, +1]:
+        # Optimize by brute force with multiprocessing
+        cost_wrapper = partial(cost, args=[modes_A, modes_B, t_reference, δϕ_factor, δΨ_factor, normalization])
+
+        if nprocs == -1:
+            cost_brute_force = [cost_wrapper(δt_δϕ) for δt_δϕ in δt_δϕ_brute_force]
+        else:
+            if nprocs is None:
+                nprocs = mp.cpu_count()
+
+            pool = mp.Pool(processes=nprocs)
+            cost_brute_force = pool.map(cost_wrapper, δt_δϕ_brute_force)
+            pool.close()
+            pool.join()
+
+        δt_δϕ = δt_δϕ_brute_force[np.argmin(cost_brute_force)]
+
+        # Optimize explicitly
+        optimum = least_squares(cost_wrapper, δt_δϕ, bounds=[(δt_lower, 0), (δt_upper, 2 * np.pi)], max_nfev=50000)
+        optimums.append(optimum)
+
+        h_A_prime = h_A.copy()
+        h_A_prime.t = h_A.t - optimum.x[0]
+        h_A_prime.data = h_A[:, 2 : ell_max + 1].data * np.exp(1j * optimum.x[1]) ** δϕ_factor * δΨ_factor
+        h_A_prime.ell_min = 2
+        h_A_prime.ell_max = ell_max
+        h_A_primes.append(h_A_prime)
+
+    idx = np.argmin(abs(np.array([optimum.cost for optimum in optimums])))
+
+    return h_A_primes[idx], optimums[idx].fun, optimums[idx]
+
+
+def rel_err_between_abds(abd1, abd2, t1, t2):
+    t_array = abd1.t[np.argmin(abs(abd1.t - t1)) : np.argmin(abs(abd1.t - t2)) + 1]
+
+    abd1_interp = abd1.interpolate(t_array)
+    abd2_interp = abd2.interpolate(t_array)
+
+    abd_diff = abd1_interp.copy()
+    abd_diff.sigma -= abd2_interp.sigma
+    abd_diff.psi0 -= abd2_interp.psi0
+    abd_diff.psi1 -= abd2_interp.psi1
+    abd_diff.psi2 -= abd2_interp.psi2
+    abd_diff.psi3 -= abd2_interp.psi3
+    abd_diff.psi4 -= abd2_interp.psi4
+
+    rel_err = 0
+    rel_err += integrate(abd_diff.sigma.norm(), abd_diff.t)[-1] / (
+        integrate(abd1_interp.sigma.norm(), abd1_interp.t)[-1] or (abd1_interp.t[-1] - abd1_interp.t[0])
+    )
+    rel_err += integrate(abd_diff.psi0.norm(), abd_diff.t)[-1] / (
+        integrate(abd1_interp.psi0.norm(), abd1_interp.t)[-1] or (abd1_interp.t[-1] - abd1_interp.t[0])
+    )
+    rel_err += integrate(abd_diff.psi1.norm(), abd_diff.t)[-1] / (
+        integrate(abd1_interp.psi1.norm(), abd1_interp.t)[-1] or (abd1_interp.t[-1] - abd1_interp.t[0])
+    )
+    rel_err += integrate(abd_diff.psi2.norm(), abd_diff.t)[-1] / (
+        integrate(abd1_interp.psi2.norm(), abd1_interp.t)[-1] or (abd1_interp.t[-1] - abd1_interp.t[0])
+    )
+    rel_err += integrate(abd_diff.psi3.norm(), abd_diff.t)[-1] / (
+        integrate(abd1_interp.psi3.norm(), abd1_interp.t)[-1] or (abd1_interp.t[-1] - abd1_interp.t[0])
+    )
+    rel_err += integrate(abd_diff.psi4.norm(), abd_diff.t)[-1] / (
+        integrate(abd1_interp.psi4.norm(), abd1_interp.t)[-1] or (abd1_interp.t[-1] - abd1_interp.t[0])
+    )
+    rel_err /= 6
+
+    return rel_err
+
+
+def map_to_abd_frame(
+    self,
+    target_abd,
+    t_0=0,
+    padding_time=250,
+    N_itr_maxes={
+        "abd": 2,
+        "superrest": 2,
+        "CoM_transformation": 10,
+        "rotation": 10,
+        "supertranslation": 10,
+    },
+    rel_err_tols={"CoM_transformation": 1e-12, "rotation": 1e-12, "supertranslation": 1e-12},
+    ell_max=None,
+    alpha_ell_max=None,
+    fix_time_phase_freedom=True,
+    nprocs=4,
+    print_conv=False,
+):
+    """Transform an abd object to a target abd object using data at t=0.
+
+    This computes the transformations necessary to map an abd object to a target abd object.
+    It uses the function map_to_bms_frame with the target charges computed from the target abd object.
+
+    Parameters
+    ----------
+    target_abd : AsymptoticBondiData
+        Target AsymptoticBondiData to map to.
+    t_0: float, optional
+        When to map to the target BMS frame.
+        Default is 0.
+    padding_time : float, optional
+        Amount by which to pad around t=0 to speed up computations, i.e.,
+        distance from t=0 in each direction to be included in self.interpolate(...).
+        This also determines the range over which certain BMS charges will be computed.
+        Default is 250.
+    N_itr_maxes : dict, optional
+        Maximum number of iterations to perform for each transformation.
+        For 'abd' , this is the number of iterations to use for the abd procedure.
+        For 'superrest', this is the number of iterations to use for the superrest procedure.
+        For the other options, these are the number of iterations to use for finding each individual transformation.
+        Default is
+        N_itr_maxes = {
+            'abd':                2,
+            'superrest':          2,
+            'CoM_transformation': 10,
+            'rotation':           10,
+            'supertranslation':   10,
+        }.
+    rel_err_tols : dict, optional
+        Relative error tolerances for each transformation.
+        Default is
+        rel_err_tols = {
+            'CoM_transformation': 1e-12,
+            'rotation':           1e-12,
+            'supertranslation':   1e-12
+        }.
+    ell_max : int, optional
+        Maximum ell to use for SWSH/Grid transformations.
+        Default is self.ell_max.
+    alpha_ell_max : int, optional
+        Maximum ell of the supertranslation to use.
+        Default is self.ell_max.
+    fix_time_phase_freedom : bool, optional
+        Whether or not to fix the time and phase freedom using a 2d minimization scheme.
+        Default is True.
+    nprocs : int, optional
+        Number of cpus to use during parallelization for fixing the time and phase freedom.
+        Default is 4. 'None' corresponds to the maximum number. '-1' corresponds to no parallelization.
+    print_conv: bool, defaults to False
+        Whether or not to print the termination criterion. Default is False.
+
+    Returns
+    -------
+    abd_prime : AsymptoticBondiData
+        Result of self.transform(...) where the input transformations are
+        the transformations found in the BMSTransformations object.
+    transformations : BMSTransformation
+        BMS transformation to map to the target BMS frame.
+
+    """
+    abd = self.copy()
+
+    time_translation = scri.bms_transformations.BMSTransformation()
+    if fix_time_phase_freedom:
+        # ensure that they are reasonable close
+        time_translation = scri.bms_transformations.BMSTransformation(
+            supertranslation=[
+                sf.constant_as_ell_0_mode(
+                    abd.t[np.argmax(abd.bondi_four_momentum()[:, 0])]
+                    - target_abd.t[np.argmax(target_abd.bondi_four_momentum()[:, 0])]
+                )
+            ]
+        )
+
+    BMS_transformation = (time_translation * scri.bms_transformations.BMSTransformation()).reorder(
+        ["supertranslation", "frame_rotation", "boost_velocity"]
+    )
+
+    abd_interp = abd.interpolate(
+        abd.t[
+            np.argmin(abs(abd.t - (t_0 - 1.5 * padding_time))) : np.argmin(abs(abd.t - (t_0 + 1.5 * padding_time))) + 1
+        ]
+    )
+
+    target_abd_superrest, transformation2, rel_err2 = target_abd.map_to_superrest_frame(
+        t_0=t_0, padding_time=padding_time
+    )
+
+    itr = 0
+    rel_err = np.inf
+    rel_errs = [np.inf]
+    while itr < N_itr_maxes["abd"]:
+        if itr == 0:
+            abd_interp_prime = abd_interp.transform(
+                supertranslation=BMS_transformation.supertranslation,
+                frame_rotation=BMS_transformation.frame_rotation.components,
+                boost_velocity=BMS_transformation.boost_velocity,
+            )
+
+        # find the transformations that map to the superrest frame
+        abd_interp_superrest, transformation1, rel_err1 = abd_interp_prime.map_to_superrest_frame(
+            t_0=t_0, padding_time=padding_time
+        )
+
+        # compose these transformations in the right order
+        BMS_transformation = (transformation2.inverse() * transformation1 * BMS_transformation).reorder(
+            ["supertranslation", "frame_rotation", "boost_velocity"]
+        )
+
+        # obtain the transformed abd object
+        abd_interp_prime = abd_interp.transform(
+            supertranslation=BMS_transformation.supertranslation,
+            frame_rotation=BMS_transformation.frame_rotation.components,
+            boost_velocity=BMS_transformation.boost_velocity,
+        )
+
+        if fix_time_phase_freedom:
+            # find the time/phase transformations
+            news_interp_prime = scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM(
+                2.0 * abd_interp_prime.sigma.bar, dataType=scri.h
+            )
+            target_news = scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM(
+                2.0 * target_abd.sigma.bar, dataType=scri.h
+            )
+
+            _, rel_err, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
+                news_interp_prime,
+                target_news,
+                t_0 - padding_time,
+                t_0 + padding_time,
+                n_brute_force_δt=None,
+                n_brute_force_δϕ=None,
+                include_modes=None,
+                nprocs=nprocs,
+            )
+
+            time_phase_transformation = scri.bms_transformations.BMSTransformation(
+                supertranslation=[sf.constant_as_ell_0_mode(res.x[0])],
+                frame_rotation=quaternion.from_rotation_vector(res.x[1] * np.array([0, 0, 1])).components,
+            )
+
+            BMS_transformation = (time_phase_transformation * BMS_transformation).reorder(
+                ["supertranslation", "frame_rotation", "boost_velocity"]
+            )
+
+            # obtain the transformed abd object
+            abd_interp_prime = abd_interp.transform(
+                supertranslation=BMS_transformation.supertranslation,
+                frame_rotation=BMS_transformation.frame_rotation.components,
+                boost_velocity=BMS_transformation.boost_velocity,
+            )
+        else:
+            rel_err = rel_err_between_abds(target_abd, abd_interp_prime, t_0 - padding_time, t_0 + padding_time)
+
+        if rel_err < min(rel_errs):
+            best_BMS_transformation = BMS_transformation.copy()
+            best_rel_err = rel_err
+        rel_errs.append(rel_err)
+
+        itr += 1
+
+    if print_conv:
+        if not itr < N_itr_max:
+            print(f"BMS: maximum number of iterations reached; the min error was {best_rel_err}.")
+        else:
+            print(f"BMS: tolerance achieved in {itr} iterations!")
+
+    abd_prime = abd.transform(
+        supertranslation=best_BMS_transformation.supertranslation,
+        frame_rotation=best_BMS_transformation.frame_rotation.components,
+        boost_velocity=best_BMS_transformation.boost_velocity,
+    )
+
+    return abd_prime, best_BMS_transformation, best_rel_err
diff --git a/scri/asymptotic_bondi_data/map_to_superrest_frame.py b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
index 16d5aa54..1d326cc5 100644
--- a/scri/asymptotic_bondi_data/map_to_superrest_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
@@ -1,5 +1,3 @@
-import copy
-import warnings
 import numpy as np
 
 import sxs
@@ -10,9 +8,8 @@
 from scri.asymptotic_bondi_data.bms_charges import charge_vector_from_aspect
 
 import quaternion
-from quaternion.calculus import derivative
+from quaternion.calculus import indefinite_integral as integrate
 
-from scipy.integrate import simpson
 from scipy.interpolate import CubicSpline
 
 
@@ -31,7 +28,7 @@ def MT_to_WM(h_mts, sxs_version=False, dataType=scri.h):
     if not sxs_version:
         h = scri.WaveformModes(
             t=h_mts.t,
-            data=np.array(h_mts)[:, sf.LM_index(abs(h_mts.s), -abs(h_mts.s), 0) :],
+            data=np.array(h_mts)[:, LM_index(abs(h_mts.s), -abs(h_mts.s), 0) :],
             ell_min=abs(h_mts.s),
             ell_max=h_mts.ell_max,
             frameType=scri.Inertial,
@@ -153,7 +150,7 @@ def compute_bondi_rest_mass_and_conformal_factor(PsiM, ell_max):
     return M_Grid, K_Grid
 
 
-def compute_Moreschi_supermomentum(abd, alpha, ell_max):
+def compute_Moreschi_supermomentum(PsiM, alpha, ell_max):
     """Compute the Moreschi supermomentum in a supertranslated frame.
 
     The transformation of the Moreschi supermomentum can be found in
@@ -169,7 +166,6 @@ def compute_Moreschi_supermomentum(abd, alpha, ell_max):
     ell_max: int
         Maximum ell to use when converting data from SWSHs to data on the spherical grid.
     """
-    PsiM = abd.supermomentum("Moreschi")
     M_Grid, K_Grid = compute_bondi_rest_mass_and_conformal_factor(np.array(PsiM), ell_max)
 
     PsiM_Grid = PsiM.grid().real
@@ -252,18 +248,39 @@ def supertranslation_to_map_to_superrest_frame(
         Whether or not to print the termination criterion. Default is False.
     """
     alpha_Grid = np.zeros((2 * ell_max + 1, 2 * ell_max + 1), dtype=float)
+    best_alpha_Grid = np.zeros((2 * ell_max + 1, 2 * ell_max + 1), dtype=float)
+
+    PsiM = abd.supermomentum("Moreschi")
 
     itr = 0
     rel_err = np.inf
-    rel_errs = []
+    rel_errs = [np.inf]
     while itr < N_itr_max and not rel_err < rel_err_tol:
         prev_alpha_Grid = sf.SWSH_grids.Grid(alpha_Grid.copy(), spin_weight=0)
 
-        PsiM = compute_Moreschi_supermomentum(abd, alpha_Grid, ell_max)
+        if itr == 0:
+            PsiM_interp = compute_Moreschi_supermomentum(PsiM, alpha_Grid, ell_max)
+
+            M_Grid, K_Grid = compute_bondi_rest_mass_and_conformal_factor(np.array(PsiM_interp), ell_max)
+
+            if target_PsiM is not None:
+                target_PsiM_Grid = WM_to_MT(target_PsiM).grid().real
+                target_PsiM_Grid_interp = sf.SWSH_grids.Grid(
+                    np.zeros((2 * target_PsiM.ell_max + 1, 2 * target_PsiM.ell_max + 1), dtype=float), spin_weight=0
+                )
+                for i in range(2 * target_PsiM.ell_max + 1):
+                    for j in range(2 * target_PsiM.ell_max + 1):
+                        alpha_i_j = prev_alpha_Grid[i, j]
+                        target_PsiM_Grid_interp[i, j] = CubicSpline(target_PsiM.t, target_PsiM_Grid[:, i, j])(alpha_i_j)
+                M_Grid = -target_PsiM_Grid_interp.view(np.ndarray)
+
+        alpha_Grid += compute_alpha_perturbation(PsiM_interp, M_Grid, K_Grid, ell_max)
+
+        PsiM_interp = compute_Moreschi_supermomentum(PsiM, alpha_Grid, ell_max)
 
-        M_Grid, K_Grid = compute_bondi_rest_mass_and_conformal_factor(np.array(PsiM), ell_max)
+        M_Grid, K_Grid = compute_bondi_rest_mass_and_conformal_factor(np.array(PsiM_interp), ell_max)
 
-        if target_PsiM != None:
+        if target_PsiM is not None:
             target_PsiM_Grid = WM_to_MT(target_PsiM).grid().real
             target_PsiM_Grid_interp = sf.SWSH_grids.Grid(
                 np.zeros((2 * target_PsiM.ell_max + 1, 2 * target_PsiM.ell_max + 1), dtype=float), spin_weight=0
@@ -274,20 +291,15 @@ def supertranslation_to_map_to_superrest_frame(
                     target_PsiM_Grid_interp[i, j] = CubicSpline(target_PsiM.t, target_PsiM_Grid[:, i, j])(alpha_i_j)
             M_Grid = -target_PsiM_Grid_interp.view(np.ndarray)
 
-        alpha_Grid += compute_alpha_perturbation(PsiM, M_Grid, K_Grid, ell_max)
-
-        rel_err = (
-            spinsfast.map2salm(
-                (
-                    abs(sf.SWSH_grids.Grid(alpha_Grid.copy(), spin_weight=0) - prev_alpha_Grid)
-                    / abs(sf.SWSH_grids.Grid(alpha_Grid.copy(), spin_weight=0))
-                ).view(np.ndarray),
-                0,
-                ell_max,
-            )[LM_index(0, 0, 0)]
-            / np.sqrt(4 * np.pi)
-        ).real
-        rel_errs.append(rel_err)
+            target_PsiM_interp = compute_Moreschi_supermomentum(target_PsiM, alpha_Grid, ell_max)
+
+            rel_err = np.linalg.norm(PsiM_interp[4:] - target_PsiM_interp[4:]) / np.linalg.norm(target_PsiM_interp[4:])
+        else:
+            rel_err = np.linalg.norm(PsiM_interp[4:])
+
+        if rel_err < min(rel_errs):
+            best_alpha_Grid = alpha_Grid.copy()
+            rel_errs.append(rel_err)
 
         itr += 1
 
@@ -299,10 +311,10 @@ def supertranslation_to_map_to_superrest_frame(
         else:
             print(f"supertranslation: tolerance achieved in {itr} iterations!")
 
-    supertranslation = spinsfast.map2salm(alpha_Grid.view(np.ndarray), 0, ell_max)
+    supertranslation = spinsfast.map2salm(best_alpha_Grid.view(np.ndarray), 0, ell_max)
     supertranslation[0:4] = 0
 
-    return supertranslation, rel_errs
+    return scri.bms_transformations.BMSTransformation(supertranslation=supertranslation), rel_errs
 
 
 def transformation_from_CoM_charge(G, t):
@@ -319,17 +331,16 @@ def transformation_from_CoM_charge(G, t):
     """
     polynomial_fit = np.polyfit(t, G, deg=1)
 
-    CoM_transformation = {
-        "space_translation": polynomial_fit[1],
-        "boost_velocity": polynomial_fit[0]
-    }
+    CoM_transformation = scri.bms_transformations.BMSTransformation(
+        supertranslation=-np.insert(sf.vector_as_ell_1_modes(polynomial_fit[1]), 0, 0),
+        boost_velocity=polynomial_fit[0],
+        order=["supertranslation", "boost_velocity", "frame_rotation"],
+    )
 
     return CoM_transformation
 
 
-def com_transformation_to_map_to_superrest_frame(
-    abd, N_itr_max=10, rel_err_tol=1e-12, ell_max=None, space_translation=True, boost_velocity=True, print_conv=False
-):
+def com_transformation_to_map_to_superrest_frame(abd, N_itr_max=10, rel_err_tol=1e-12, print_conv=False):
     """Determine the space translation and boost needed to map an abd object to the superrest frame.
 
     These are found through an iterative solve; e.g., compute the transformations needed to minimize
@@ -342,134 +353,75 @@ def com_transformation_to_map_to_superrest_frame(
     abd: AsymptoticBondiData
         AsymptoticBondiData object from which the Moreschi supermomentum will be computed.
     N_itr_max: int, defaults to 10
-        Maximum numebr of iterations to perform. Default is 10.
+        Maximum number of iterations to perform. Default is 10.
     rel_err_tol: float, defaults to 1e-12
         Minimum relativie error tolerance between transformation iterations. Default is 1e-12.
-    ell_max: int, defaults to 12
-        Maximum ell to use when converting data from SWSHs to data on the spherical grid. Default is 12.
-    space_translation: bool, defaults to True
-        Whether or not to return the space translation.
-    boost_velocity: bool, defaults to True
-        Whether or not to return the boost velocity.
     print_conv: bool, defaults to False
         Whether or not to print the termination criterion. Default is False.
     """
-    if not space_translation and not boost_velocity:
-        raise ValueError("space_translation and boost_velocity cannot both be False.")
-
-    CoM_transformation = {
-        "space_translation": np.zeros(3),
-        "boost_velocity": np.zeros(3)
-    }
+    CoM_transformation = scri.bms_transformations.BMSTransformation()
+    best_CoM_transformation = scri.bms_transformations.BMSTransformation()
 
     itr = 0
-    rel_err = np.array(2 * [np.inf])
-    rel_errs = []
-    while itr < N_itr_max and not (rel_err < rel_err_tol).sum() == int(space_translation) + int(boost_velocity):
-        prev_CoM_transformation = copy.deepcopy(CoM_transformation)
-
+    rel_err = np.inf
+    rel_errs = [np.inf]
+    while itr < N_itr_max and not rel_err < rel_err_tol:
         if itr == 0:
             abd_prime = abd.copy()
-        else:
-            abd_prime = abd.transform(
-                space_translation=CoM_transformation["space_translation"],
-                boost_velocity=CoM_transformation["boost_velocity"],
-            )
-
-        G_prime = abd_prime.bondi_CoM_charge() / abd_prime.bondi_four_momentum()[:, 0, None]
+            G_prime = abd_prime.bondi_CoM_charge() / abd_prime.bondi_four_momentum()[:, 0, None]
 
         new_CoM_transformation = transformation_from_CoM_charge(G_prime, abd_prime.t)
-        for transformation in CoM_transformation:
-            if not space_translation:
-                if transformation == "space_translation":
-                    continue
-            if not boost_velocity:
-                if transformation == "boost_velocity":
-                    continue
-            CoM_transformation[transformation] += new_CoM_transformation[transformation]
-
-        for i, transformation in enumerate(CoM_transformation):
-            with warnings.catch_warnings():
-                warnings.filterwarnings("ignore", message="invalid value encountered in double_scalars")
-                err = np.sqrt(
-                    abs(
-                        np.dot(
-                            CoM_transformation[transformation] - prev_CoM_transformation[transformation],
-                            CoM_transformation[transformation] - prev_CoM_transformation[transformation],
-                        )
-                        / np.dot(CoM_transformation[transformation], CoM_transformation[transformation])
-                    )
-                )
+        CoM_transformation = (new_CoM_transformation * CoM_transformation).reorder(
+            ["supertranslation", "frame_rotation", "boost_velocity"]
+        )
+        # remove proper supertranslation and frame rotation components
+        CoM_transformation.supertranslation[4:] *= 0
+        CoM_transformation.frame_rotation = quaternion.quaternion(1, 0, 0, 0)
+
+        abd_prime = abd.transform(
+            supertranslation=CoM_transformation.supertranslation,
+            frame_rotation=CoM_transformation.frame_rotation.components,
+            boost_velocity=CoM_transformation.boost_velocity,
+        )
 
-                if err == np.inf:
-                    err = 1
-                rel_err[i] = err
+        G_prime = abd_prime.bondi_CoM_charge() / abd_prime.bondi_four_momentum()[:, 0, None]
 
-        rel_errs.append(copy.deepcopy(rel_err))
+        rel_err = integrate(np.linalg.norm(G_prime, axis=-1), abd_prime.t)[-1] / (abd_prime.t[-1] - abd_prime.t[0])
+        if rel_err < min(rel_errs):
+            best_CoM_transformation = CoM_transformation.copy()
+        rel_errs.append(rel_err)
 
         itr += 1
 
-    transformation_name = "CoM"
-    if not space_translation:
-        transformation_name = "boost"
-        transformation = CoM_transformation["boost_velocity"]
-        rel_errs = np.array(rel_errs)[:, 1:]
-    elif not boost_velocity:
-        transformation_name = "space_translation"
-        transformation = CoM_transformation["space_translation"]
-        rel_errs = np.array(rel_errs)[:, :1]
-    else:
-        transformation = CoM_transformation
-        rel_errs = np.array(rel_errs)
-
     if print_conv:
         if not itr < N_itr_max:
-            print(
-                f"{transformation_name}: maximum number of iterations reached; the min error was {min(rel_errs.flatten())}."
-            )
+            print(f"CoM: maximum number of iterations reached; the min error was {min(rel_errs)}.")
         else:
-            print(f"{transformation_name}: tolerance achieved in {itr} iterations!")
-
-    return transformation, rel_errs
-
-
-def rotation_from_spin_charge(chi, t):
-    """Obtain the rotation from the remnant BH's spin vector.
-
-    This finds the rotation that aligns the z-component of the spin vector with the z-axis.
+            print(f"CoM: tolerance achieved in {itr} iterations!")
 
-    Parameters
-    ----------
-    chi: ndarray, real, shape (..., 3)
-        Remnant BH's spin vector.
-    t: ndarray, real
-        Time array corresponding to the size of the spin vector.
-    """
-    chi_f = quaternion.quaternion(*chi[np.argmin(abs(t))]).normalized()
-    return (1 - chi_f * quaternion.z).normalized()
+    return best_CoM_transformation, rel_errs
 
 
-def rotation_from_omega_vectors(omega, target_omega, t=None):
-    """Obtain the rotation from the two angular velocity vectors.
+def rotation_from_vectors(vector, target_vector, t=None):
+    """Obtain the rotation from the two vectors.
 
-    This finds the rotation that best aligns two angular velocity vectors over the time interval t.
+    This finds the rotation that best aligns two vectors over the time interval t.
 
     Parameters
     ----------
-    omega: ndarray, real, shape (..., 3)
-        Angular velocity vector.
-    target_omega: ndarray, real, shape (..., 3)
-        Target angular velocity vector.
+    vector: ndarray, real, shape (..., 3)
+        Angular vector.
+    target_vector: ndarray, real, shape (..., 3)
+        Target vector.
     t: ndarray, real
         Time array corresponding to the size of the spin vector. Default is None.
     """
-    q = quaternion.optimal_alignment_in_Euclidean_metric(omega, target_omega, t=None)
-    return q
+    q = quaternion.optimal_alignment_in_Euclidean_metric(vector, target_vector, t=t).inverse()
 
+    return scri.bms_transformations.BMSTransformation(frame_rotation=q.components)
 
-def rotation_to_map_to_superrest_frame(
-    abd, target_h=None, N_itr_max=10, rel_err_tol=[1e-12, 1e-5], ell_max=None, print_conv=False
-):
+
+def rotation_to_map_to_superrest_frame(abd, target_strain=None, N_itr_max=10, rel_err_tol=1e-12, print_conv=False):
     """Determine the rotation needed to map an abd object to the superrest frame.
 
     This is found through an iterative solve; e.g., compute the transformation needed to align
@@ -487,86 +439,101 @@ def rotation_to_map_to_superrest_frame(
     ----------
     abd: AsymptoticBondiData
         AsymptoticBondiData object from which the Moreschi supermomentum will be computed.
-    target_h: WaveformModes
+    target_strain: WaveformModes, optional
         WaveformModes object from which the target angular velocity will be computed. Default is None.
     N_itr_max: int, defaults to 10
         Maximum number of iterations to perform. Default is 10.
-    rel_err_tol: array_like, defaults to [1e-12, 1e-5]
+    rel_err_tol: array_like, defaults to 1e-12
         First value is minimum relative error tolerance between transformation iterations; second value is the
         minimum relative error tolerance between the NR angular velocity and the target angular velocity.
-        Default is [1e-12, 1e-5].
-    ell_max: int, defaults to 12
-        Maximum ell to use when converting data from SWSHs to data on the spherical grid. Default is 12.
+        Default is 1e-12.
     print_conv: bool, defaults to False
         Whether or not to print the termination criterion. Default is False.
     """
-    rotation = quaternion.quaternion(1, 0, 0, 0)
+    rotation_transformation = scri.bms_transformations.BMSTransformation()
+    best_rotation_transformation = scri.bms_transformations.BMSTransformation()
 
-    # if there is no target_h, just map the spin charge to the z-axis;
+    # if there is no target_strain, just map the spin charge to the z-axis;
     # otherwise, align the angular momentum fluxes over the full window.
-    if target_h == None:
+    if target_strain is not None:
+        target_news = MT_to_WM(
+            scri.ModesTimeSeries(
+                sf.SWSH_modes.Modes(
+                    target_strain.data_dot,
+                    spin_weight=-2,
+                    ell_min=2,
+                    ell_max=target_strain.ell_max,
+                    multiplication_truncator=max,
+                ),
+                time=target_strain.t,
+            ).dot,
+            dataType=scri.hdot,
+        )
+
+        target_omega = target_news.angular_velocity()
+        target_omega_spline = CubicSpline(
+            target_omega / np.linalg.norm(target_omega, axis=-1)[:, None], target_strain.t
+        )
+
         itr = 0
         rel_err = np.inf
-        rel_errs = []
-        rel_errs_omega = []
-        while itr < N_itr_max and not rel_err < rel_err_tol[0]:
-            prev_rotation = copy.deepcopy(rotation)
-
+        rel_errs = [np.inf]
+        while itr < N_itr_max and not rel_err < rel_err_tol:
             if itr == 0:
                 abd_prime = abd.copy()
-            else:
-                abd_prime = abd.transform(frame_rotation=rotation.components)
+                news = MT_to_WM(2.0 * abd.sigma.bar.dot, dataType=scri.hdot)
+                omega = news.angular_velocity()
+                omega = omega / np.linalg.norm(omega, axis=-1)[:, None]
 
-            chi_prime = abd_prime.bondi_dimensionless_spin()
+            rotation_transformation = (
+                rotation_from_vectors(omega, target_omega_spline(news.t), news.t) * rotation_transformation
+            ).reorder(["supertranslation", "frame_rotation", "boost_velocity"])
+            # remove supertranslation and CoM components
+            rotation_transformation.supertranslation *= 0
+            rotation_transformation.boost_velocity *= 0
 
-            rotation = rotation_from_spin_charge(chi_prime, abd_prime.t) * rotation
+            abd_prime = abd.transform(frame_rotation=rotation_transformation.frame_rotation.components)
 
-            rel_err = quaternion.rotation_intrinsic_distance(rotation, prev_rotation)
+            news = MT_to_WM(2.0 * abd_prime.sigma.bar.dot, dataType=scri.hdot)
+            omega = news.angular_velocity()
+            omega = omega / np.linalg.norm(omega, axis=-1)[:, None]
 
+            rel_err = integrate(np.linalg.norm(omega - target_omega_spline(news.t), axis=-1), news.t)[-1] / (
+                abd_prime.t[-1] - abd_prime.t[0]
+            )
+            if rel_err < min(rel_errs):
+                best_rotation_transformation = rotation_transformation.copy()
             rel_errs.append(rel_err)
 
             itr += 1
-
-        if print_conv:
-            if not itr < N_itr_max:
-                print(
-                    f"rotation: maximum number of iterations reached; the min error was {min(np.array(rel_errs).flatten())}."
-                )
-            else:
-                print(f"rotation: tolerance achieved in {itr} iterations!")
-
-        return rotation.components, rel_errs
     else:
-        target_news = MT_to_WM(WM_to_MT(target_h).dot, dataType=scri.hdot)
-
-        target_omega = target_news.angular_velocity()
-        target_omega = target_omega / np.linalg.norm(target_omega, axis=1)[:, None]
-
         itr = 0
         rel_err = np.inf
-        rel_err_omega = np.inf
-        rel_errs = []
-        rel_errs_omega = []
-        while itr < N_itr_max and not (rel_err < rel_err_tol[0] or rel_err_omega < rel_err_tol[1]):
-            prev_rotation = copy.deepcopy(rotation)
-
+        rel_errs = [np.inf]
+        while itr < N_itr_max and not rel_err < rel_err_tol:
             if itr == 0:
                 abd_prime = abd.copy()
-            else:
-                abd_prime = abd.transform(frame_rotation=rotation.components)
-            news = MT_to_WM(2.0 * abd_prime.sigma.bar.dot, dataType=scri.hdot)
+                chi_prime = abd_prime.bondi_dimensionless_spin()
+                chi_prime = chi_prime / np.linalg.norm(chi_prime, axis=-1)[:, None]
 
-            omega = news.angular_velocity()
-            omega = omega / np.linalg.norm(omega, axis=1)[:, None]
+            rotation_transformation = (
+                rotation_from_vectors(chi_prime, [[0, 0, 1]] * abd_prime.t.size, abd_prime.t) * rotation_transformation
+            ).reorder(["supertranslation", "frame_rotation", "boost_velocity"])
+            # remove supertranslation and CoM components
+            rotation_transformation.supertranslation *= 0
+            rotation_transformation.boost_velocity *= 0
 
-            rotation = rotation_from_omega_vectors(omega, target_omega, news.t) * rotation
+            abd_prime = abd.transform(frame_rotation=rotation_transformation.frame_rotation.components)
 
-            rel_err = quaternion.rotor_intrinsic_distance(rotation, prev_rotation)
-            rel_err_omega = simpson(np.linalg.norm(omega - target_omega, axis=1), news.t) / simpson(
-                np.linalg.norm(target_omega, axis=1), news.t
-            )
+            chi_prime = abd_prime.bondi_dimensionless_spin()
+            chi_prime = chi_prime / np.linalg.norm(chi_prime, axis=-1)[:, None]
 
-            rel_errs.append([rel_err, rel_err_omega])
+            rel_err = integrate(np.linalg.norm(chi_prime - [[0, 0, 1]] * abd_prime.t.size, axis=-1), abd_prime.t)[
+                -1
+            ] / (abd_prime.t[-1] - abd_prime.t[0])
+            if rel_err < min(rel_errs):
+                best_rotation_transformation = rotation_transformation.copy()
+            rel_errs.append(rel_err)
 
             itr += 1
 
@@ -578,7 +545,7 @@ def rotation_to_map_to_superrest_frame(
         else:
             print(f"rotation: tolerance achieved in {itr} iterations!")
 
-    return rotation.components, rel_errs
+    return best_rotation_transformation, rel_errs
 
 
 def time_translation(abd, t_0=0):
@@ -602,24 +569,102 @@ def time_translation(abd, t_0=0):
     return abd_prime
 
 
+def rotation(abd, phi=0):
+    """Rotate an abd object.
+
+    This is faster than using abd.transform().
+    """
+    q = quaternion.from_rotation_vector(-phi * np.array([0, 0, 1]))
+
+    h = MT_to_WM(2.0 * abd.sigma.bar, False, scri.h)
+    Psi4 = MT_to_WM(0.5 * (-np.sqrt(2)) ** 4 * abd.psi4, False, scri.psi4)
+    Psi3 = MT_to_WM(0.5 * (-np.sqrt(2)) ** 3 * abd.psi3, False, scri.psi3)
+    Psi2 = MT_to_WM(0.5 * (-np.sqrt(2)) ** 2 * abd.psi2, False, scri.psi2)
+    Psi1 = MT_to_WM(0.5 * (-np.sqrt(2)) ** 1 * abd.psi1, False, scri.psi1)
+    Psi0 = MT_to_WM(0.5 * (-np.sqrt(2)) ** 0 * abd.psi0, False, scri.psi0)
+
+    h.rotate_physical_system(q)
+    Psi4.rotate_physical_system(q)
+    Psi3.rotate_physical_system(q)
+    Psi2.rotate_physical_system(q)
+    Psi1.rotate_physical_system(q)
+    Psi0.rotate_physical_system(q)
+
+    abd_rot = abd.copy()
+    abd_rot.sigma = 0.5 * WM_to_MT(h).bar
+    abd_rot.psi4 = 2 * (-1.0 / np.sqrt(2)) ** 4 * WM_to_MT(Psi4)
+    abd_rot.psi3 = 2 * (-1.0 / np.sqrt(2)) ** 3 * WM_to_MT(Psi3)
+    abd_rot.psi2 = 2 * (-1.0 / np.sqrt(2)) ** 2 * WM_to_MT(Psi2)
+    abd_rot.psi1 = 2 * (-1.0 / np.sqrt(2)) ** 1 * WM_to_MT(Psi1)
+    abd_rot.psi0 = 2 * (-1.0 / np.sqrt(2)) ** 0 * WM_to_MT(Psi0)
+
+    return abd_rot
+
+
+def rel_err_for_abd_in_superrest(abd, target_PsiM, target_strain):
+    G = abd.bondi_CoM_charge() / abd.bondi_four_momentum()[:, 0, None]
+    rel_err_CoM_transformation = integrate(np.linalg.norm(G, axis=-1), abd.t)[-1] / (abd.t[-1] - abd.t[0])
+
+    if target_strain is not None:
+        target_news = MT_to_WM(
+            scri.ModesTimeSeries(
+                sf.SWSH_modes.Modes(
+                    target_strain.data_dot,
+                    spin_weight=-2,
+                    ell_min=2,
+                    ell_max=target_strain.ell_max,
+                    multiplication_truncator=max,
+                ),
+                time=target_strain.t,
+            ).dot,
+            dataType=scri.hdot,
+        )
+
+        target_omega = target_news.angular_velocity()
+        target_omega_spline = CubicSpline(
+            target_omega / np.linalg.norm(target_omega, axis=-1)[:, None], target_strain.t
+        )
+
+        news = MT_to_WM(2.0 * abd.sigma.bar.dot, dataType=scri.hdot)
+        omega = news.angular_velocity()
+        omega = omega / np.linalg.norm(omega, axis=-1)[:, None]
+
+        rel_err_rotation = integrate(np.linalg.norm(omega - target_omega_spline(news.t), axis=-1), news.t)[-1] / (
+            abd.t[-1] - abd.t[0]
+        )
+    else:
+        spin = abd.bondi_dimensionless_spin()
+        spin = spin / np.linalg.norm(spin, axis=-1)[:, None]
+        rel_err_rotation = integrate(np.linalg.norm(spin - [[0, 0, 1]] * abd.t.size, axis=-1), abd.t)[-1] / (
+            abd.t[-1] - abd.t[0]
+        )
+
+    if target_PsiM is not None:
+        rel_err_PsiM = np.linalg.norm(
+            np.array(abd.supermomentum("Moreschi"))[np.argmin(abs(abd.t - 0)), 4:]
+            - np.array(target_PsiM)[np.argmin(abs(abd.t - 0)), 4:]
+        )
+    else:
+        rel_err_PsiM = np.linalg.norm(np.array(abd.supermomentum("Moreschi"))[np.argmin(abs(abd.t - 0)), 4:])
+
+    return rel_err_CoM_transformation, rel_err_rotation, rel_err_PsiM
+
+
 def map_to_superrest_frame(
     self,
     t_0=0,
-    target_h_input=None,
     target_PsiM_input=None,
+    target_strain_input=None,
+    padding_time=250,
     N_itr_maxes={
-        "supertranslation": 10,
-        "com_transformation": 10,
+        "superrest": 2,
+        "CoM_transformation": 10,
         "rotation": 10,
+        "supertranslation": 10,
     },
-    rel_err_tols={
-        "supertranslation": 1e-12,
-        "com_transformation": 1e-12,
-        "rotation": [1e-12, 1e-5],
-    },
+    rel_err_tols={"CoM_transformation": 1e-12, "rotation": 1e-12, "supertranslation": 1e-12},
     ell_max=None,
     alpha_ell_max=None,
-    padding_time=250,
     print_conv=False,
 ):
     """Transform an abd object to the superrest frame.
@@ -640,28 +685,36 @@ def map_to_superrest_frame(
     t_0 : float, optional
         When to map to the superrest frame.
         Default is 0.
-    target_h_input : WaveformModes, optional
-        Target strain used to constrain the rotation
-        freedom via the angular momentum flux.
-        Default is aligning to the z-axis.
     target_PsiM_input : WaveformModes, optional
         Target Moreschi supermomentum to map to.
         Default is 0.
+    target_strain_input : WaveformModes, optional
+        Target strain used to constrain the rotation
+        freedom via the angular momentum flux.
+        Default is aligning to the z-axis.
+    padding_time : float, optional
+        Amount by which to pad around t_0 to speed up computations, i.e.,
+        distance from t_0 in each direction to be included in self.interpolate(...).
+        This also determines the range over which certain BMS charges will be computed.
+        Default is 250.
     N_itr_maxes : dict, optional
         Maximum number of iterations to perform for each transformation.
+        For 'superrest', this is the number of iterations to use for the superrest procedure.
+        For the other options, these are the number of iterations to use for finding each individual transformation.
         Default is
         N_itr_maxes = {
+            'superrest':          2,
+            'CoM_transformation': 10,
+            'rotation':           10,
             'supertranslation':   10,
-            'com_transformation': 10,
-            'rotation':           10
         }.
     rel_err_tols : dict, optional
         Relative error tolerances for each transformation.
         Default is
-        N_itr_maxes = {
-            'supertranslation':   1e-12,
-            'com_transformation': 1e-12,
-            'rotation':           1e-12
+        rel_err_tols = {
+            'CoM_transformation': 1e-12,
+            'rotation':           1e-12,
+            'supertranslation':   1e-12
         }.
     ell_max : int, optional
         Maximum ell to use for SWSH/Grid transformations.
@@ -669,137 +722,146 @@ def map_to_superrest_frame(
     alpha_ell_max : int, optional
         Maximum ell of the supertranslation to use.
         Default is self.ell_max.
-    padding_time : float, optional
-        Amount by which to pad around t_0 to speed up computations, i.e.,
-        distance from t_0 in each direction to be included in self.interpolate(...).
-        This also determines the range over which certain BMS charges will be computed.
-        Default is 100.
+    print_conv: bool, defaults to False
+        Whether or not to print the termination criterion. Default is False.
 
     Returns
     -------
     abd_prime : AsymptoticBondiData
         Result of self.transform(...) where the input transformations are
-        the transformations found in the transformations dictionary.
-    transformations : dict
-        Dictionary of transformations and their relative errors whose keys are
-            * 'transformations'
-                * 'space_translation'
-                * 'supertranslation'
-                * 'frame_rotation'
-                * 'boost_velocity'
-            * 'rel_errs'
-                (same as above)
+        the transformations found in the BMSTransformations object.
+    transformations : BMSTransformation
+        BMS transformation to map to the target BMS frame.
 
     """
-    # apply a time translation so that we're mapping
-    # to the superrest frame at u = 0
-    abd = time_translation(self, t_0)
+    abd = self.copy()
 
-    if target_h_input != None:
-        target_h = target_h_input.copy()
-        target_h.t -= t_0
+    if target_strain_input is not None:
+        target_strain = target_strain_input.copy()
+        target_strain.t -= t_0
     else:
-        target_h = None
+        target_strain = None
 
-    if target_PsiM_input != None:
+    if target_PsiM_input is not None:
         target_PsiM = target_PsiM_input.copy()
         target_PsiM.t -= t_0
     else:
         target_PsiM = None
 
-    if ell_max == None:
+    if ell_max is None:
         ell_max = abd.ell_max
 
-    if alpha_ell_max == None:
+    if alpha_ell_max is None:
         alpha_ell_max = ell_max
 
     abd_interp = abd.interpolate(
-        abd.t[np.argmin(abs(abd.t - (-padding_time))) : np.argmin(abs(abd.t - (+padding_time))) + 1]
+        abd.t[np.argmin(abs(abd.t - (t_0 - padding_time))) : np.argmin(abs(abd.t - (t_0 + padding_time))) + 1]
     )
 
-    # space_translation
-    space_translation, space_rel_errs = com_transformation_to_map_to_superrest_frame(
-        abd_interp,
-        N_itr_max=N_itr_maxes["com_transformation"],
-        rel_err_tol=rel_err_tols["com_transformation"],
-        ell_max=ell_max,
-        space_translation=True,
-        boost_velocity=False,
-        print_conv=print_conv,
+    # apply a time translation so that we're mapping
+    # to the superrest frame at u = 0
+    time_translation = scri.bms_transformations.BMSTransformation(supertranslation=[sf.constant_as_ell_0_mode(t_0)])
+    BMS_transformation = time_translation * scri.bms_transformations.BMSTransformation().reorder(
+        ["supertranslation", "frame_rotation", "boost_velocity"]
     )
 
-    # supertranslation
-    abd_prime = abd_interp.transform(space_translation=space_translation)
+    itr = 0
+    rel_err = [np.inf, np.inf, np.inf]
+    rel_errs = [[np.inf, np.inf, np.inf]]
+    best_rel_err = [np.inf, np.inf, np.inf]
+    while itr < N_itr_maxes["superrest"] and not (
+        rel_err[0] < rel_err_tols["CoM_transformation"]
+        and rel_err[1] < rel_err_tols["rotation"]
+        and rel_err[2] < rel_err_tols["supertranslation"]
+    ):
+        if itr == 0:
+            abd_interp_prime = abd_interp.transform(
+                supertranslation=BMS_transformation.supertranslation,
+                frame_rotation=BMS_transformation.frame_rotation.components,
+                boost_velocity=BMS_transformation.boost_velocity,
+            )
+
+        # CoM transformation
+        CoM_transformation, CoM_rel_errs = com_transformation_to_map_to_superrest_frame(
+            abd_interp_prime,
+            N_itr_max=N_itr_maxes["CoM_transformation"],
+            rel_err_tol=rel_err_tols["CoM_transformation"],
+            print_conv=print_conv,
+        )
 
-    alpha, alpha_rel_errs = supertranslation_to_map_to_superrest_frame(
-        abd_prime,
-        target_PsiM,
-        N_itr_max=N_itr_maxes["supertranslation"],
-        rel_err_tol=rel_err_tols["supertranslation"],
-        ell_max=ell_max,
-        print_conv=print_conv,
-    )
+        BMS_transformation = (CoM_transformation * BMS_transformation).reorder(
+            ["supertranslation", "frame_rotation", "boost_velocity"]
+        )
 
-    alpha[0] = 0
-    alpha[1:4] = sf.vector_as_ell_1_modes(space_translation)
-
-    # rotation
-    abd_prime = abd_interp.transform(supertranslation=alpha[: LM_index(alpha_ell_max, alpha_ell_max, 0) + 1])
-
-    # do interpolation here because .transform() changes the size of the time array
-    target_h_interp = None
-    if target_h != None:
-        target_h_interp = target_h.interpolate(abd_prime.t)
-    rotation, rot_rel_errs = rotation_to_map_to_superrest_frame(
-        abd_prime,
-        target_h_interp,
-        N_itr_max=N_itr_maxes["rotation"],
-        rel_err_tol=rel_err_tols["rotation"],
-        ell_max=ell_max,
-        print_conv=print_conv,
-    )
+        abd_interp_prime = abd_interp.transform(
+            supertranslation=BMS_transformation.supertranslation,
+            frame_rotation=BMS_transformation.frame_rotation.components,
+            boost_velocity=BMS_transformation.boost_velocity,
+        )
 
-    # com_transformation
-    abd_prime = abd_interp.transform(
-        supertranslation=alpha[: LM_index(alpha_ell_max, alpha_ell_max, 0) + 1], frame_rotation=rotation
-    )
+        # rotation
+        rotation, rot_rel_errs = rotation_to_map_to_superrest_frame(
+            abd_interp_prime,
+            target_strain=target_strain,
+            N_itr_max=N_itr_maxes["rotation"],
+            rel_err_tol=rel_err_tols["rotation"],
+            print_conv=print_conv,
+        )
 
-    CoM_transformation, CoM_rel_errs = com_transformation_to_map_to_superrest_frame(
-        abd_prime,
-        N_itr_max=N_itr_maxes["com_transformation"],
-        rel_err_tol=rel_err_tols["com_transformation"],
-        ell_max=ell_max,
-        space_translation=True,
-        boost_velocity=True,
-        print_conv=print_conv,
+        BMS_transformation = (rotation * BMS_transformation).reorder(
+            ["supertranslation", "frame_rotation", "boost_velocity"]
+        )
+
+        abd_interp_prime = abd_interp.transform(
+            supertranslation=BMS_transformation.supertranslation,
+            frame_rotation=BMS_transformation.frame_rotation.components,
+            boost_velocity=BMS_transformation.boost_velocity,
+        )
+
+        # supertranslation
+        supertranslation, supertranslation_rel_errs = supertranslation_to_map_to_superrest_frame(
+            abd_interp_prime,
+            target_PsiM,
+            N_itr_max=N_itr_maxes["supertranslation"],
+            rel_err_tol=rel_err_tols["supertranslation"],
+            ell_max=ell_max,
+            print_conv=print_conv,
+        )
+
+        BMS_transformation = (supertranslation * BMS_transformation).reorder(
+            ["supertranslation", "frame_rotation", "boost_velocity"]
+        )
+
+        abd_interp_prime = abd_interp.transform(
+            supertranslation=BMS_transformation.supertranslation,
+            frame_rotation=BMS_transformation.frame_rotation.components,
+            boost_velocity=BMS_transformation.boost_velocity,
+        )
+
+        rel_err = rel_err_for_abd_in_superrest(abd_interp_prime, target_PsiM, target_strain)
+        if np.mean(rel_err) < min([np.mean(r) for r in rel_errs]):
+            best_BMS_transformation = BMS_transformation.copy()
+            best_rel_err = rel_err
+        rel_errs.append(rel_err)
+
+        itr += 1
+
+    if print_conv:
+        if not itr < N_itr_maxes["superrest"]:
+            print(f"superrest: maximum number of iterations reached; the min error was {best_rel_err}.")
+        else:
+            print(f"superrest: tolerance achieved in {itr} iterations!")
+
+    # undo the time translation
+    best_BMS_transformation = (time_translation.inverse() * best_BMS_transformation).reorder(
+        ["supertranslation", "frame_rotation", "boost_velocity"]
     )
 
     # transform abd
     abd_prime = abd.transform(
-        supertranslation=alpha[: LM_index(alpha_ell_max, alpha_ell_max, 0) + 1], frame_rotation=rotation
-    )
-    abd_prime = abd_prime.transform(
-        space_translation=CoM_transformation["space_translation"], boost_velocity=CoM_transformation["boost_velocity"]
+        supertranslation=best_BMS_transformation.supertranslation,
+        frame_rotation=best_BMS_transformation.frame_rotation.components,
+        boost_velocity=best_BMS_transformation.boost_velocity,
     )
 
-    # undo the initial time translation
-    abd_prime = time_translation(abd_prime, -t_0)
-
-    alpha[0:4] = 0
-
-    transformations = {
-        "transformations": {
-            "space_translation": space_translation,
-            "supertranslation": alpha,
-            "frame_rotation": rotation,
-            "CoM_transformation": CoM_transformation,
-        },
-        "rel_errs": {
-            "space_translation": space_rel_errs,
-            "supertranslation": alpha_rel_errs,
-            "frame_rotation": rot_rel_errs,
-            "CoM_transformation": CoM_rel_errs,
-        },
-    }
-
-    return abd_prime, transformations
+    return abd_prime, best_BMS_transformation, best_rel_err
diff --git a/tests/test_abd_frame.py b/tests/test_abd_frame.py
new file mode 100644
index 00000000..0b7e1b5d
--- /dev/null
+++ b/tests/test_abd_frame.py
@@ -0,0 +1,91 @@
+import scri
+import pytest
+import numpy as np
+import spinsfast
+import spherical_functions as sf
+import quaternion
+from quaternion.calculus import derivative
+
+ABD = scri.AsymptoticBondiData
+
+
+def kerr_schild(mass, spin, ell_max=8):
+    psi2 = np.zeros(sf.LM_total_size(0, ell_max), dtype=complex)
+    psi1 = np.zeros(sf.LM_total_size(0, ell_max), dtype=complex)
+    psi0 = np.zeros(sf.LM_total_size(0, ell_max), dtype=complex)
+
+    # In the Moreschi-Boyle convention
+    psi2[0] = -sf.constant_as_ell_0_mode(mass)
+    psi1[2] = -np.sqrt(2) * (3j * spin / 2) * (np.sqrt((8 / 3) * np.pi))
+    psi0[6] = 2 * (3 * spin**2 / mass / 2) * (np.sqrt((32 / 15) * np.pi))
+
+    return psi2, psi1, psi0
+
+
+def test_abd_to_abd():
+    mass = 2.0
+    spin = 0.456
+    ell_max = 8
+    u = np.linspace(-100, 100, num=5000)
+    psi2, psi1, psi0 = kerr_schild(mass, spin, ell_max)
+    abd = ABD.from_initial_values(u, ell_max=ell_max, psi2=psi2, psi1=psi1)
+
+    supertranslation = np.array(
+        [0.0, 3e-2 - 1j * 5e-3, 1e-3, -3e-2 - 1j * 5e-3, 2e-4 + 1j * 1e-4, 1j * 3e-3, 1e-2, 1j * 3e-3, 2e-4 - 1j * 1e-4]
+    )
+    frame_rotation = quaternion.quaternion(1, 2, 3, 4).normalized().components
+    boost_velocity = np.array([2e-4, -3e-4, -5e-4])
+
+    abd_target = abd.transform(
+        supertranslation=supertranslation, frame_rotation=frame_rotation, boost_velocity=boost_velocity
+    )
+
+    # We don't preform time/phase fixing because this data is radiation-free so it doesn't make sense
+    abd_recovered, transformation, rel_err = abd.map_to_abd_frame(
+        abd_target,
+        t_0=0,
+        padding_time=20,
+        N_itr_maxes={"abd": 1, "superrest": 1, "CoM_transformation": 10, "rotation": 10, "supertranslation": 10},
+        fix_time_phase_freedom=False,
+        nprocs=-1,
+    )
+
+    abd_target_interp = abd_target.interpolate(
+        abd_target.t[
+            np.argmin(abs(abd_target.t - max(abd_target.t[0], abd_recovered.t[0]))) : np.argmin(
+                abs(abd_target.t - min(abd_target.t[-1], abd_recovered.t[-1]))
+            )
+            + 1
+        ]
+    )
+    abd_recovered_interp = abd_recovered.interpolate(
+        abd_recovered.t[
+            np.argmin(abs(abd_recovered.t - max(abd_target.t[0], abd_recovered.t[0]))) : np.argmin(
+                abs(abd_recovered.t - min(abd_target.t[-1], abd_recovered.t[-1]))
+            )
+            + 1
+        ]
+    )
+
+    assert np.allclose(
+        np.array(
+            [
+                abd_target_interp.sigma,
+                abd_target_interp.psi4,
+                abd_target_interp.psi3,
+                abd_target_interp.psi2,
+                abd_target_interp.psi1,
+                abd_target_interp.psi0,
+            ]
+        ),
+        np.array(
+            [
+                abd_recovered_interp.sigma,
+                abd_recovered_interp.psi4,
+                abd_recovered_interp.psi3,
+                abd_recovered_interp.psi2,
+                abd_recovered_interp.psi1,
+                abd_recovered_interp.psi0,
+            ]
+        ),
+    )
diff --git a/tests/test_superrest_frame.py b/tests/test_superrest_frame.py
index 03257ff4..fbeebd1c 100644
--- a/tests/test_superrest_frame.py
+++ b/tests/test_superrest_frame.py
@@ -3,6 +3,7 @@
 import numpy as np
 import spinsfast
 import spherical_functions as sf
+import quaternion
 from quaternion.calculus import derivative
 
 ABD = scri.AsymptoticBondiData
@@ -28,30 +29,84 @@ def test_abd_kerr_superrest_frame():
     u = np.linspace(-100, 100, num=5000)
     psi2, psi1, psi0 = kerr_schild(mass, spin, ell_max)
     abd = ABD.from_initial_values(u, ell_max=ell_max, psi2=psi2, psi1=psi1)
+
     tolerance = 1e-12
 
     supertranslation = np.array(
         [0.0, 3e-2 - 1j * 5e-3, 1e-3, -3e-2 - 1j * 5e-3, 2e-4 + 1j * 1e-4, 1j * 3e-3, 1e-2, 1j * 3e-3, 2e-4 - 1j * 1e-4]
     )
-    frame_rotation = np.array([np.sqrt(1 - (0.02**2 + 0.1**2 + 0.04**2)), 0.02, 0.1, 0.04])
+    frame_rotation = quaternion.quaternion(1, 2, 3, 4).normalized().components
     boost_velocity = np.array([2e-4, -3e-5, 2e-4])
     abd_prime = abd.transform(
         supertranslation=supertranslation, frame_rotation=frame_rotation, boost_velocity=boost_velocity
     )
 
-    abd_recovered, _ = abd_prime.map_to_superrest_frame(padding_time=20)
+    abd_recovered, transformation, rel_errs = abd_prime.map_to_superrest_frame(
+        t_0=0,
+        padding_time=20,
+        N_itr_maxes={"superrest": 1, "CoM_transformation": 10, "rotation": 10, "supertranslation": 10},
+    )
+
+    PsiM_recovered = abd_recovered.supermomentum("Moreschi")[np.argmin(abs(abd_recovered.t))]
+    PsiM_recovered[0:4] = 0
+    PsiM_recovered_S2 = spinsfast.salm2map(
+        PsiM_recovered.view(np.ndarray), 0, ell_max, 2 * ell_max + 1, 2 * ell_max + 1
+    )
+    PsiM_recovered_norm = spinsfast.map2salm((abs(PsiM_recovered_S2) ** 2).view(np.ndarray), 0, ell_max)[0] / np.sqrt(
+        4 * np.pi
+    )
+    assert np.allclose(PsiM_recovered_norm, 0.0, atol=tolerance, rtol=tolerance)
+
+    chi_recovered = abd_recovered.bondi_dimensionless_spin()
+    chi_recovered = chi_recovered / np.linalg.norm(chi_recovered, axis=-1)[:, None]
+    assert np.allclose(chi_recovered, [[0, 0, 1]] * abd_recovered.t.size, atol=tolerance, rtol=tolerance)
+
+    G_recovered = abd_recovered.bondi_CoM_charge() / abd_recovered.bondi_four_momentum()[:, 0, None]
+    assert np.allclose(G_recovered[np.argmin(abs(abd_recovered.t))], 0.0, atol=tolerance, rtol=tolerance)
+    assert np.allclose(
+        derivative(G_recovered, abd_recovered.t)[np.argmin(abs(abd_recovered.t))], 0.0, atol=tolerance, rtol=tolerance
+    )
+
+
+def test_abd_kerr_target_superrest_frame():
+    mass = 2.0
+    spin = 0.456
+    ell_max = 8
+    u = np.linspace(-100, 100, num=5000)
+    psi2, psi1, psi0 = kerr_schild(mass, spin, ell_max)
+    abd = ABD.from_initial_values(u, ell_max=ell_max, psi2=psi2, psi1=psi1)
+
+    tolerance = 1e-12
+
+    supertranslation = np.array(
+        [0.0, 3e-2 - 1j * 5e-3, 1e-3, -3e-2 - 1j * 5e-3, 2e-4 + 1j * 1e-4, 1j * 3e-3, 1e-2, 1j * 3e-3, 2e-4 - 1j * 1e-4]
+    )
+    abd_target = abd.transform(
+        supertranslation=supertranslation,
+    )
+
+    abd_recovered, transformation, rel_errs = abd.map_to_superrest_frame(
+        t_0=0,
+        padding_time=20,
+        target_PsiM_input=abd_target.supermomentum("Moreschi"),
+        N_itr_maxes={
+            "superrest": 1,
+            "CoM_transformation": 10,
+            "rotation": 10,
+            "supertranslation": 10,
+        },
+    )
 
-    PsiM = abd_recovered.supermomentum("Moreschi")[np.argmin(abs(abd_recovered.t))]
-    PsiM[0:4] = 0
-    PsiM_S2 = spinsfast.salm2map(PsiM.view(np.ndarray), 0, ell_max, 2 * ell_max + 1, 2 * ell_max + 1)
-    PsiM_norm = spinsfast.map2salm((abs(PsiM_S2) ** 2).view(np.ndarray), 0, ell_max)[0] / np.sqrt(4 * np.pi)
-    assert np.allclose(PsiM_norm, 0.0, atol=tolerance, rtol=tolerance)
+    PsiM_diff = abd_recovered.supermomentum("Moreschi")[np.argmin(abs(abd_recovered.t))]
+    PsiM_diff -= abd_target.supermomentum("Moreschi")[np.argmin(abs(abd_recovered.t))]
+    PsiM_diff[0:4] = 0
+    PsiM_diff_S2 = spinsfast.salm2map(PsiM_diff.view(np.ndarray), 0, ell_max, 2 * ell_max + 1, 2 * ell_max + 1)
+    PsiM_diff_norm = spinsfast.map2salm((abs(PsiM_diff_S2) ** 2).view(np.ndarray), 0, ell_max)[0] / np.sqrt(4 * np.pi)
+    assert np.allclose(PsiM_diff_norm, 0.0, atol=tolerance, rtol=tolerance)
 
-    chi = abd_recovered.bondi_dimensionless_spin()[np.argmin(abs(abd_recovered.t))]
-    theta = np.arccos(chi[2] / np.linalg.norm(chi))
-    if theta > np.pi / 2.0:
-        theta -= np.pi
-    assert np.allclose(theta, 0.0, atol=tolerance, rtol=tolerance)
+    chi_recovered = abd_recovered.bondi_dimensionless_spin()
+    chi_recovered = chi_recovered / np.linalg.norm(chi_recovered, axis=-1)[:, None]
+    assert np.allclose(chi_recovered, [[0, 0, 1]] * abd_recovered.t.size, atol=tolerance, rtol=tolerance)
 
     G = abd_recovered.bondi_CoM_charge() / abd_recovered.bondi_four_momentum()[:, 0, None]
     assert np.allclose(G[np.argmin(abs(abd_recovered.t))], 0.0, atol=tolerance, rtol=tolerance)

From 6c8f09e5e51953f54b9801899521c6d10f0c83f2 Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Tue, 14 Feb 2023 17:31:05 -0800
Subject: [PATCH 03/13] Fix code (for real data) by fixing time/phase freedom
 in superrest frame

---
 scri/asymptotic_bondi_data/bms_charges.py     |   2 +-
 .../asymptotic_bondi_data/map_to_abd_frame.py |  75 ++++++++--
 .../map_to_superrest_frame.py                 | 139 ++++++++++--------
 scri/asymptotic_bondi_data/transformations.py |   2 +-
 scri/waveform_base.py                         |   8 +-
 tests/test_abd_frame.py                       |   2 +-
 6 files changed, 149 insertions(+), 79 deletions(-)

diff --git a/scri/asymptotic_bondi_data/bms_charges.py b/scri/asymptotic_bondi_data/bms_charges.py
index b8c01c27..d7b46956 100644
--- a/scri/asymptotic_bondi_data/bms_charges.py
+++ b/scri/asymptotic_bondi_data/bms_charges.py
@@ -108,7 +108,7 @@ def CWWY_angular_momentum(self):
 
     """
     ell_max = 1  # Compute only the parts we need, ell<=1
-    supertranslation_potential_data = scri.asymptotic_bondi_data.map_to_bms_frame.𝔇inverse(
+    supertranslation_potential_data = scri.asymptotic_bondi_data.map_to_superrest_frame.𝔇inverse(
         np.array(self.sigma.ethbar_GHP.ethbar_GHP + self.sigma.bar.eth_GHP.eth_GHP), self.ell_max
     )
     supertranslation_potential = scri.ModesTimeSeries(
diff --git a/scri/asymptotic_bondi_data/map_to_abd_frame.py b/scri/asymptotic_bondi_data/map_to_abd_frame.py
index bec484f9..cce70439 100644
--- a/scri/asymptotic_bondi_data/map_to_abd_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_abd_frame.py
@@ -247,6 +247,7 @@ def map_to_abd_frame(
         "supertranslation": 10,
     },
     rel_err_tols={"CoM_transformation": 1e-12, "rotation": 1e-12, "supertranslation": 1e-12},
+    order=["rotation", "CoM_transformation", "supertranslation"],
     ell_max=None,
     alpha_ell_max=None,
     fix_time_phase_freedom=True,
@@ -291,6 +292,9 @@ def map_to_abd_frame(
             'rotation':           1e-12,
             'supertranslation':   1e-12
         }.
+    order : list, optional
+        Order in which to solve for the BMS transformations.
+        Default is ["rotation", "CoM_transformation", "supertranslation"].
     ell_max : int, optional
         Maximum ell to use for SWSH/Grid transformations.
         Default is self.ell_max.
@@ -317,14 +321,18 @@ def map_to_abd_frame(
     """
     abd = self.copy()
 
+    target_strain = scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM(
+        2.0 * target_abd.sigma.bar, dataType=scri.h
+    )
+
     time_translation = scri.bms_transformations.BMSTransformation()
     if fix_time_phase_freedom:
         # ensure that they are reasonable close
         time_translation = scri.bms_transformations.BMSTransformation(
             supertranslation=[
                 sf.constant_as_ell_0_mode(
-                    abd.t[np.argmax(abd.bondi_four_momentum()[:, 0])]
-                    - target_abd.t[np.argmax(target_abd.bondi_four_momentum()[:, 0])]
+                    abd.t[np.argmax(abd.bondi_four_momentum()[np.argmin(abs(abd.t - t_0)), 0])]
+                    - target_abd.t[np.argmax(target_abd.bondi_four_momentum()[np.argmin(abs(target_abd.t - t_0)), 0])]
                 )
             ]
         )
@@ -340,7 +348,18 @@ def map_to_abd_frame(
     )
 
     target_abd_superrest, transformation2, rel_err2 = target_abd.map_to_superrest_frame(
-        t_0=t_0, padding_time=padding_time
+        t_0=t_0,
+        padding_time=padding_time,
+        N_itr_maxes=N_itr_maxes,
+        rel_err_tols=rel_err_tols,
+        ell_max=ell_max,
+        alpha_ell_max=alpha_ell_max,
+        print_conv=print_conv,
+        order=order,
+    )
+
+    target_strain_superrest = scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM(
+        2.0 * target_abd_superrest.sigma.bar, dataType=scri.h
     )
 
     itr = 0
@@ -356,13 +375,44 @@ def map_to_abd_frame(
 
         # find the transformations that map to the superrest frame
         abd_interp_superrest, transformation1, rel_err1 = abd_interp_prime.map_to_superrest_frame(
-            t_0=t_0, padding_time=padding_time
+            t_0=t_0,
+            padding_time=padding_time,
+            N_itr_maxes=N_itr_maxes,
+            rel_err_tols=rel_err_tols,
+            ell_max=ell_max,
+            alpha_ell_max=alpha_ell_max,
+            print_conv=print_conv,
+            order=order,
         )
 
+        if fix_time_phase_freedom:
+            # 2d align now that we're in the superrest frame
+            strain_interp_superrest = scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM(
+                2.0 * abd_interp_superrest.sigma.bar, dataType=scri.h
+            )
+
+            _, rel_err, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
+                strain_interp_superrest,
+                target_strain_superrest,
+                t_0 - padding_time,
+                t_0 + padding_time,
+                n_brute_force_δt=None,
+                n_brute_force_δϕ=None,
+                include_modes=None,
+                nprocs=nprocs,
+            )
+
+            time_phase_transformation = scri.bms_transformations.BMSTransformation(
+                supertranslation=[sf.constant_as_ell_0_mode(res.x[0])],
+                frame_rotation=quaternion.from_rotation_vector(res.x[1] * np.array([0, 0, 1])).components,
+            )
+        else:
+            time_phase_transformation = scri.bms_transformations.BMSTransformation()
+
         # compose these transformations in the right order
-        BMS_transformation = (transformation2.inverse() * transformation1 * BMS_transformation).reorder(
-            ["supertranslation", "frame_rotation", "boost_velocity"]
-        )
+        BMS_transformation = (
+            transformation2.inverse() * (time_phase_transformation * (transformation1 * BMS_transformation))
+        ).reorder(["supertranslation", "frame_rotation", "boost_velocity"])
 
         # obtain the transformed abd object
         abd_interp_prime = abd_interp.transform(
@@ -373,16 +423,13 @@ def map_to_abd_frame(
 
         if fix_time_phase_freedom:
             # find the time/phase transformations
-            news_interp_prime = scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM(
+            strain_interp_prime = scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM(
                 2.0 * abd_interp_prime.sigma.bar, dataType=scri.h
             )
-            target_news = scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM(
-                2.0 * target_abd.sigma.bar, dataType=scri.h
-            )
 
             _, rel_err, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
-                news_interp_prime,
-                target_news,
+                strain_interp_prime,
+                target_strain,
                 t_0 - padding_time,
                 t_0 + padding_time,
                 n_brute_force_δt=None,
@@ -417,7 +464,7 @@ def map_to_abd_frame(
         itr += 1
 
     if print_conv:
-        if not itr < N_itr_max:
+        if not itr < N_itr_maxes["abd"]:
             print(f"BMS: maximum number of iterations reached; the min error was {best_rel_err}.")
         else:
             print(f"BMS: tolerance achieved in {itr} iterations!")
diff --git a/scri/asymptotic_bondi_data/map_to_superrest_frame.py b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
index 1d326cc5..bc066e54 100644
--- a/scri/asymptotic_bondi_data/map_to_superrest_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
@@ -85,7 +85,7 @@ def 𝔇(h, ell_max):
             value = 0
         else:
             value = ((ell + 2) * (ell + 1) * (ell) * (ell - 1)) / 4.0
-        h_with_operator[LM_index(ell, -ell, 0) : LM_index(ell, ell, 0) + 1] *= value
+        h_with_operator[..., LM_index(ell, -ell, 0) : LM_index(ell, ell, 0) + 1] *= value
 
     return h_with_operator
 
@@ -299,7 +299,7 @@ def supertranslation_to_map_to_superrest_frame(
 
         if rel_err < min(rel_errs):
             best_alpha_Grid = alpha_Grid.copy()
-            rel_errs.append(rel_err)
+        rel_errs.append(rel_err)
 
         itr += 1
 
@@ -402,6 +402,21 @@ def com_transformation_to_map_to_superrest_frame(abd, N_itr_max=10, rel_err_tol=
     return best_CoM_transformation, rel_errs
 
 
+def rotation_from_spin_charge(chi, t):
+    """Obtain the rotation from the remnant BH's spin vector.
+    This finds the rotation that aligns the z-component of the spin vector with the z-axis.
+    Parameters
+    ----------
+    chi: ndarray, real, shape (..., 3)
+        Remnant BH's spin vector.
+    t: ndarray, real
+        Time array corresponding to the size of the spin vector.
+    """
+    chi_f = quaternion.quaternion(*chi[np.argmin(abs(t))]).normalized()
+    q = (1 - chi_f * quaternion.z).normalized().components
+    return scri.bms_transformations.BMSTransformation(frame_rotation=q)
+
+
 def rotation_from_vectors(vector, target_vector, t=None):
     """Obtain the rotation from the two vectors.
 
@@ -517,7 +532,7 @@ def rotation_to_map_to_superrest_frame(abd, target_strain=None, N_itr_max=10, re
                 chi_prime = chi_prime / np.linalg.norm(chi_prime, axis=-1)[:, None]
 
             rotation_transformation = (
-                rotation_from_vectors(chi_prime, [[0, 0, 1]] * abd_prime.t.size, abd_prime.t) * rotation_transformation
+                rotation_from_spin_charge(chi_prime, abd_prime.t) * rotation_transformation
             ).reorder(["supertranslation", "frame_rotation", "boost_velocity"])
             # remove supertranslation and CoM components
             rotation_transformation.supertranslation *= 0
@@ -663,8 +678,10 @@ def map_to_superrest_frame(
         "supertranslation": 10,
     },
     rel_err_tols={"CoM_transformation": 1e-12, "rotation": 1e-12, "supertranslation": 1e-12},
+    order=["rotation", "CoM_transformation", "supertranslation"],
     ell_max=None,
     alpha_ell_max=None,
+    fix_time_phase_freedom=False,
     print_conv=False,
 ):
     """Transform an abd object to the superrest frame.
@@ -716,12 +733,18 @@ def map_to_superrest_frame(
             'rotation':           1e-12,
             'supertranslation':   1e-12
         }.
+    order : list, optional
+        Order in which to solve for the BMS transformations.
+        Default is ["rotation", "CoM_transformation", "supertranslation"].
     ell_max : int, optional
         Maximum ell to use for SWSH/Grid transformations.
         Default is self.ell_max.
     alpha_ell_max : int, optional
         Maximum ell of the supertranslation to use.
         Default is self.ell_max.
+    fix_time_phase_freedom : bool, optional
+        Whether or not to fix the time and phase freedom using a 2d minimization scheme.
+        Default is True.
     print_conv: bool, defaults to False
         Whether or not to print the termination criterion. Default is False.
 
@@ -781,62 +804,62 @@ def map_to_superrest_frame(
                 boost_velocity=BMS_transformation.boost_velocity,
             )
 
-        # CoM transformation
-        CoM_transformation, CoM_rel_errs = com_transformation_to_map_to_superrest_frame(
-            abd_interp_prime,
-            N_itr_max=N_itr_maxes["CoM_transformation"],
-            rel_err_tol=rel_err_tols["CoM_transformation"],
-            print_conv=print_conv,
-        )
-
-        BMS_transformation = (CoM_transformation * BMS_transformation).reorder(
-            ["supertranslation", "frame_rotation", "boost_velocity"]
-        )
-
-        abd_interp_prime = abd_interp.transform(
-            supertranslation=BMS_transformation.supertranslation,
-            frame_rotation=BMS_transformation.frame_rotation.components,
-            boost_velocity=BMS_transformation.boost_velocity,
-        )
-
-        # rotation
-        rotation, rot_rel_errs = rotation_to_map_to_superrest_frame(
-            abd_interp_prime,
-            target_strain=target_strain,
-            N_itr_max=N_itr_maxes["rotation"],
-            rel_err_tol=rel_err_tols["rotation"],
-            print_conv=print_conv,
-        )
-
-        BMS_transformation = (rotation * BMS_transformation).reorder(
-            ["supertranslation", "frame_rotation", "boost_velocity"]
-        )
-
-        abd_interp_prime = abd_interp.transform(
-            supertranslation=BMS_transformation.supertranslation,
-            frame_rotation=BMS_transformation.frame_rotation.components,
-            boost_velocity=BMS_transformation.boost_velocity,
-        )
-
-        # supertranslation
-        supertranslation, supertranslation_rel_errs = supertranslation_to_map_to_superrest_frame(
-            abd_interp_prime,
-            target_PsiM,
-            N_itr_max=N_itr_maxes["supertranslation"],
-            rel_err_tol=rel_err_tols["supertranslation"],
-            ell_max=ell_max,
-            print_conv=print_conv,
-        )
-
-        BMS_transformation = (supertranslation * BMS_transformation).reorder(
-            ["supertranslation", "frame_rotation", "boost_velocity"]
-        )
+        for transformation in order:
+            if transformation == "supertranslation":
+                new_transformation, supertranslation_rel_errs = supertranslation_to_map_to_superrest_frame(
+                    abd_interp_prime,
+                    target_PsiM,
+                    N_itr_max=N_itr_maxes["supertranslation"],
+                    rel_err_tol=rel_err_tols["supertranslation"],
+                    ell_max=ell_max,
+                    print_conv=print_conv,
+                )
+            elif transformation == "rotation":
+                new_transformation, rot_rel_errs = rotation_to_map_to_superrest_frame(
+                    abd_interp_prime,
+                    target_strain=target_strain,
+                    N_itr_max=N_itr_maxes["rotation"],
+                    rel_err_tol=rel_err_tols["rotation"],
+                    print_conv=print_conv,
+                )
+            elif transformation == "CoM_transformation":
+                new_transformation, CoM_rel_errs = com_transformation_to_map_to_superrest_frame(
+                    abd_interp_prime,
+                    N_itr_max=N_itr_maxes["CoM_transformation"],
+                    rel_err_tol=rel_err_tols["CoM_transformation"],
+                    print_conv=print_conv,
+                )
+            elif transformation == "time_phase":
+                if target_strain is not None:
+                    strain_interp_prime = scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM(
+                        2.0 * abd_interp_prime.sigma.bar, dataType=scri.h
+                    )
+
+                    _, rel_err, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
+                        strain_interp_prime,
+                        target_strain,
+                        t_0 - padding_time,
+                        t_0 + padding_time,
+                        n_brute_force_δt=None,
+                        n_brute_force_δϕ=None,
+                        include_modes=None,
+                        nprocs=nprocs,
+                    )
+
+                    new_transformation = scri.bms_transformations.BMSTransformation(
+                        supertranslation=[sf.constant_as_ell_0_mode(res.x[0])],
+                        frame_rotation=quaternion.from_rotation_vector(res.x[1] * np.array([0, 0, 1])).components,
+                    )
+
+            BMS_transformation = (new_transformation * BMS_transformation).reorder(
+                ["supertranslation", "frame_rotation", "boost_velocity"]
+            )
 
-        abd_interp_prime = abd_interp.transform(
-            supertranslation=BMS_transformation.supertranslation,
-            frame_rotation=BMS_transformation.frame_rotation.components,
-            boost_velocity=BMS_transformation.boost_velocity,
-        )
+            abd_interp_prime = abd_interp.transform(
+                supertranslation=BMS_transformation.supertranslation,
+                frame_rotation=BMS_transformation.frame_rotation.components,
+                boost_velocity=BMS_transformation.boost_velocity,
+            )
 
         rel_err = rel_err_for_abd_in_superrest(abd_interp_prime, target_PsiM, target_strain)
         if np.mean(rel_err) < min([np.mean(r) for r in rel_errs]):
diff --git a/scri/asymptotic_bondi_data/transformations.py b/scri/asymptotic_bondi_data/transformations.py
index d45b4e48..0eda6349 100644
--- a/scri/asymptotic_bondi_data/transformations.py
+++ b/scri/asymptotic_bondi_data/transformations.py
@@ -40,7 +40,7 @@ def _process_transformation_kwargs(input_ell_max, **kwargs):
                 i_neg = sf.LM_index(ell, -m, 0)
                 a = supertranslation[i_pos]
                 b = supertranslation[i_neg]
-                if abs(a - (-1.0) ** m * b.conjugate()) > 3e-16 + 1e-15 * abs(b):
+                if abs(a - (-1.0) ** m * b.conjugate()) > 1e-12 + 1e-12 * abs(b):
                     raise ValueError(
                         f"\nsupertranslation[{i_pos}]={a}  # (ell,m)=({ell},{m})\n"
                         + "supertranslation[{}]={}  # (ell,m)=({},{})\n".format(i_neg, b, ell, -m)
diff --git a/scri/waveform_base.py b/scri/waveform_base.py
index 32bfd3bf..7a48486e 100644
--- a/scri/waveform_base.py
+++ b/scri/waveform_base.py
@@ -961,7 +961,7 @@ def interpolate(self, tprime):
         W.t = np.copy(tprime)
         W.frame = quaternion.squad(self.frame, self.t, W.t)
         W.data = np.empty((W.n_times,) + self.data.shape[1:], dtype=self.data.dtype)
-        W.data_2d[:] = CubicSpline(self.t, self.data_2d.view(float))(W.t).view(complex)
+        W.data_2d[:] = CubicSpline(self.t, self.data_2d.view(complex))(W.t).view(complex)
         W.__history_depth__ -= 1
         W._append_history(f"{W} = {self}.interpolate({tprime})")
         return W
@@ -1011,12 +1011,12 @@ def SI_units(self, current_unit_mass_in_solar_masses, distance_from_source_in_me
         # quantities are curvature quantities, so they have dimensions 1/m^2, and thus have `m_scaling=2`.
         if self.r_is_scaled_out:
             if self.m_is_scaled_out:
-                amplitude_scaling = (R_over_M ** -self.r_scaling) * (M_in_meters ** -self.m_scaling)
+                amplitude_scaling = (R_over_M**-self.r_scaling) * (M_in_meters**-self.m_scaling)
             else:
-                amplitude_scaling = R_over_M ** -self.r_scaling
+                amplitude_scaling = R_over_M**-self.r_scaling
         else:
             if self.m_is_scaled_out:
-                amplitude_scaling = M_in_meters ** -self.m_scaling
+                amplitude_scaling = M_in_meters**-self.m_scaling
             else:
                 amplitude_scaling = 1.0
 
diff --git a/tests/test_abd_frame.py b/tests/test_abd_frame.py
index 0b7e1b5d..776f67e5 100644
--- a/tests/test_abd_frame.py
+++ b/tests/test_abd_frame.py
@@ -40,7 +40,7 @@ def test_abd_to_abd():
         supertranslation=supertranslation, frame_rotation=frame_rotation, boost_velocity=boost_velocity
     )
 
-    # We don't preform time/phase fixing because this data is radiation-free so it doesn't make sense
+    # We don't perform time/phase fixing because this data is radiation-free so it doesn't make sense
     abd_recovered, transformation, rel_err = abd.map_to_abd_frame(
         abd_target,
         t_0=0,

From b86a17d10dd5147bd78e5f20af107190a382229f Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Mon, 6 Mar 2023 11:40:45 -0800
Subject: [PATCH 04/13] Fixing DInverse bug

---
 scri/asymptotic_bondi_data/map_to_superrest_frame.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/scri/asymptotic_bondi_data/map_to_superrest_frame.py b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
index bc066e54..bb5f23a2 100644
--- a/scri/asymptotic_bondi_data/map_to_superrest_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
@@ -98,7 +98,7 @@ def 𝔇inverse(h, ell_max):
             value = 0
         else:
             value = 4.0 / ((ell + 2) * (ell + 1) * (ell) * (ell - 1))
-        h_with_operator[LM_index(ell, -ell, 0) : LM_index(ell, ell, 0) + 1] *= value
+        h_with_operator[..., LM_index(ell, -ell, 0) : LM_index(ell, ell, 0) + 1] *= value
 
     return h_with_operator
 

From 56fc309a476cc2e274292fe48739ed23273047ec Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Wed, 15 Mar 2023 10:42:25 -0700
Subject: [PATCH 05/13] Changing waveform_base.py to use new version of
 CubicSpline; imposing reality condition on supertranslations

---
 scri/asymptotic_bondi_data/transformations.py | 11 ++++-------
 scri/waveform_base.py                         |  2 +-
 2 files changed, 5 insertions(+), 8 deletions(-)

diff --git a/scri/asymptotic_bondi_data/transformations.py b/scri/asymptotic_bondi_data/transformations.py
index 0eda6349..20c9d97a 100644
--- a/scri/asymptotic_bondi_data/transformations.py
+++ b/scri/asymptotic_bondi_data/transformations.py
@@ -33,19 +33,16 @@ def _process_transformation_kwargs(input_ell_max, **kwargs):
                 "documentation for details).\n           Thus, it must be an array with length given by a "
                 "perfect square; its length is {len(supertranslation)}"
             )
-        # Check that the resulting supertranslation will be real
+        # Impose reality 
         for ell in range(ell_max_supertranslation + 1):
             for m in range(ell + 1):
                 i_pos = sf.LM_index(ell, m, 0)
                 i_neg = sf.LM_index(ell, -m, 0)
                 a = supertranslation[i_pos]
                 b = supertranslation[i_neg]
-                if abs(a - (-1.0) ** m * b.conjugate()) > 1e-12 + 1e-12 * abs(b):
-                    raise ValueError(
-                        f"\nsupertranslation[{i_pos}]={a}  # (ell,m)=({ell},{m})\n"
-                        + "supertranslation[{}]={}  # (ell,m)=({},{})\n".format(i_neg, b, ell, -m)
-                        + "Will result in an imaginary supertranslation."
-                    )
+                supertranslation[i_pos] = (a + (-1.0) ** m * b.conjugate()) / 2.0
+                supertranslation[i_neg] = (-1.0) ** m * supertranslation[i_pos].conjugate()
+                
     spacetime_translation = np.zeros((4,), dtype=float)
     spacetime_translation[0] = sf.constant_from_ell_0_mode(supertranslation[0]).real
     spacetime_translation[1:4] = -sf.vector_from_ell_1_modes(supertranslation[1:4]).real
diff --git a/scri/waveform_base.py b/scri/waveform_base.py
index 7a48486e..ec345c4d 100644
--- a/scri/waveform_base.py
+++ b/scri/waveform_base.py
@@ -961,7 +961,7 @@ def interpolate(self, tprime):
         W.t = np.copy(tprime)
         W.frame = quaternion.squad(self.frame, self.t, W.t)
         W.data = np.empty((W.n_times,) + self.data.shape[1:], dtype=self.data.dtype)
-        W.data_2d[:] = CubicSpline(self.t, self.data_2d.view(complex))(W.t).view(complex)
+        W.data_2d[:] = CubicSpline(self.t, self.data_2d)(W.t)
         W.__history_depth__ -= 1
         W._append_history(f"{W} = {self}.interpolate({tprime})")
         return W

From fc592a50ce500fc29fcf974ba7086f8c220a2f78 Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Mon, 27 Mar 2023 13:56:27 -0700
Subject: [PATCH 06/13] Change align2d code to use sxs

---
 .../asymptotic_bondi_data/map_to_abd_frame.py | 186 +-----------------
 1 file changed, 4 insertions(+), 182 deletions(-)

diff --git a/scri/asymptotic_bondi_data/map_to_abd_frame.py b/scri/asymptotic_bondi_data/map_to_abd_frame.py
index cce70439..26c73b59 100644
--- a/scri/asymptotic_bondi_data/map_to_abd_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_abd_frame.py
@@ -12,190 +12,12 @@
 
 from scipy.interpolate import CubicSpline
 
+from sxs.waveforms.alignment import align2d
+
 from functools import partial
 
 from scri.asymptotic_bondi_data.map_to_superrest_frame import time_translation, rotation
 
-
-def cost(δt_δϕ, args):
-    modes_A, modes_B, t_reference, δϕ_factor, δΨ_factor, normalization = args
-
-    # Take the sqrt because least_squares squares the inputs...
-    diff = integrate(
-        np.sum(
-            abs(modes_A(t_reference + δt_δϕ[0]) * np.exp(1j * δt_δϕ[1]) ** δϕ_factor * δΨ_factor - modes_B) ** 2, axis=1
-        ),
-        t_reference,
-    )[-1]
-    return np.sqrt(diff / normalization)
-
-
-def align2d(h_A, h_B, t1, t2, n_brute_force_δt=None, n_brute_force_δϕ=None, include_modes=None, nprocs=4):
-    """Align waveforms by shifting in time and phase
-
-    This function determines the optimal time and phase offset to apply to `h_A` by minimizing
-    the averaged (over time) L² norm (over the sphere) of the difference of the h_Aveforms.
-
-    The integral is taken from time `t1` to `t2`.
-
-    Note that the input waveforms are assumed to be initially aligned at least well
-    enough that:
-
-      1) the time span from `t1` to `t2` in the two waveforms will overlap at
-         least slightly after the second waveform is shifted in time; and
-      2) waveform `h_B` contains all the times corresponding to `t1` to `t2`
-         in waveform `h_A`.
-
-    The first of these can usually be assured by simply aligning the peaks prior to
-    calling this function:
-
-        h_A.t -= h_A.max_norm_time() - h_B.max_norm_time()
-
-    The second assumption will be satisfied as long as `t1` is not too close to the
-    beginning of `h_B` and `t2` is not too close to the end.
-
-    Parameters
-    ----------
-    h_A : WaveformModes
-    h_B : WaveformModes
-        Waveforms to be aligned
-    t1 : float
-    t2 : float
-        Beginning and end of integration interval
-    n_brute_force_δt : int, optional
-        Number of evenly spaced δt values between (t1-t2) and (t2-t1) to sample
-        for the initial guess.  By default, this is just the maximum number of
-        time steps in the range (t1, t2) in the input waveforms.  If this is
-        too small, an incorrect local minimum may be found.
-    n_brute_force_δϕ : int, optional
-        Number of evenly spaced δϕ values between 0 and 2π to sample
-        for the initial guess.  By default, this is 2 * ell_max + 1.
-    include_modes: list, optional
-        A list containing the (ell, m) modes to be included in the L² norm.
-    nprocs: int, optional
-        Number of cpus to use.
-        Default is 4. 'None' corresponds to the maximum number. '-1' corresponds to no parallelization.
-
-    Returns
-    -------
-    optimum: OptimizeResult
-        Result of scipy.optimize.least_squares
-    h_A_prime: WaveformModes
-        Resulting waveform after transforming `h_A` using `optimum`
-
-    Notes
-    -----
-    Choosing the time interval is usually the most difficult choice to make when
-    aligning waveforms.  Assuming you want to align during inspiral, the times
-    must span sufficiently long that the waveforms' norm (equivalently, orbital
-    frequency changes) significantly from `t1` to `t2`.  This means that you
-    cannot always rely on a specific number of orbits, for example.  Also note
-    that neither number should be too close to the beginning or end of either
-    waveform, to provide some "wiggle room".
-
-    Precession generally causes no problems for this function.  In principle,
-    eccentricity, center-of-mass offsets, boosts, or other supertranslations could
-    cause problems, but this function begins with a brute-force method of finding
-    the optimal time offset that will avoid local minima in all but truly
-    outrageous situations.  In particular, as long as `t1` and `t2` are separated
-    by enough, there should never be a problem.
-
-    """
-    from scipy.optimize import least_squares
-
-    import multiprocessing as mp
-
-    h_A_copy = h_A.copy()
-    h_B_copy = h_B.copy()
-
-    # Check that (t1, t2) makes sense and is actually contained in both waveforms
-    if t2 <= t1:
-        raise ValueError(f"(t1,t2)=({t1}, {t2}) is out of order")
-    if h_A_copy.t[0] > t1 or h_A_copy.t[-1] < t2:
-        raise ValueError(
-            f"(t1,t2)=({t1}, {t2}) not contained in h_A_copy.t, which spans ({h_A_copy.t[0]}, {h_A_copy.t[-1]})"
-        )
-    if h_B_copy.t[0] > t1 or h_B_copy.t[-1] < t2:
-        raise ValueError(
-            f"(t1,t2)=({t1}, {t2}) not contained in h_B_copy.t, which spans ({h_B_copy.t[0]}, {h_B_copy.t[-1]})"
-        )
-
-    # Figure out time offsets to try
-    δt_lower = max(t1 - t2, h_A_copy.t[0] - t1)
-    δt_upper = min(t2 - t1, h_A_copy.t[-1] - t2)
-
-    # We'll start by brute forcing, sampling time offsets evenly at as many
-    # points as there are time steps in (t1,t2) in the input waveforms
-    if n_brute_force_δt is None:
-        n_brute_force_δt = max(
-            sum((h_A_copy.t >= t1) & (h_A_copy.t <= t2)), sum((h_B_copy.t >= t1) & (h_B_copy.t <= t2))
-        )
-    δt_brute_force = np.linspace(δt_lower, δt_upper, num=n_brute_force_δt)
-
-    if n_brute_force_δϕ is None:
-        n_brute_force_δϕ = 2 * h_A_copy.ell_max + 1
-    δϕ_brute_force = np.linspace(0, 2 * np.pi, n_brute_force_δϕ, endpoint=False)
-
-    δt_δϕ_brute_force = np.array(np.meshgrid(δt_brute_force, δϕ_brute_force)).T.reshape(-1, 2)
-
-    t_reference = h_B_copy.t[np.argmin(abs(h_B_copy.t - t1)) : np.argmin(abs(h_B_copy.t - t2)) + 1]
-
-    # Remove certain modes, if requested
-    ell_max = min(h_A_copy.ell_max, h_B_copy.ell_max)
-    if include_modes != None:
-        for L in range(2, ell_max + 1):
-            for M in range(-L, L + 1):
-                if not (L, M) in include_modes:
-                    h_A_copy.data[:, LM_index(L, M, h_A_copy.ell_min)] *= 0
-                    h_B_copy.data[:, LM_index(L, M, h_B_copy.ell_min)] *= 0
-
-    # Define the cost function
-    modes_A = CubicSpline(h_A_copy.t, h_A_copy[:, 2 : ell_max + 1].data)
-    modes_B = CubicSpline(h_B_copy.t, h_B_copy[:, 2 : ell_max + 1].data)(t_reference)
-
-    normalization = integrate(CubicSpline(h_B_copy.t, h_B_copy[:, 2 : ell_max + 1].norm())(t_reference), t_reference)[
-        -1
-    ]
-    if normalization == 0:
-        normalization = t_reference[-1] - t_reference[0]
-
-    δϕ_factor = np.array([M for L in range(h_A_copy.ell_min, ell_max + 1) for M in range(-L, L + 1)])
-
-    optimums = []
-    h_A_primes = []
-    for δΨ_factor in [-1, +1]:
-        # Optimize by brute force with multiprocessing
-        cost_wrapper = partial(cost, args=[modes_A, modes_B, t_reference, δϕ_factor, δΨ_factor, normalization])
-
-        if nprocs == -1:
-            cost_brute_force = [cost_wrapper(δt_δϕ) for δt_δϕ in δt_δϕ_brute_force]
-        else:
-            if nprocs is None:
-                nprocs = mp.cpu_count()
-
-            pool = mp.Pool(processes=nprocs)
-            cost_brute_force = pool.map(cost_wrapper, δt_δϕ_brute_force)
-            pool.close()
-            pool.join()
-
-        δt_δϕ = δt_δϕ_brute_force[np.argmin(cost_brute_force)]
-
-        # Optimize explicitly
-        optimum = least_squares(cost_wrapper, δt_δϕ, bounds=[(δt_lower, 0), (δt_upper, 2 * np.pi)], max_nfev=50000)
-        optimums.append(optimum)
-
-        h_A_prime = h_A.copy()
-        h_A_prime.t = h_A.t - optimum.x[0]
-        h_A_prime.data = h_A[:, 2 : ell_max + 1].data * np.exp(1j * optimum.x[1]) ** δϕ_factor * δΨ_factor
-        h_A_prime.ell_min = 2
-        h_A_prime.ell_max = ell_max
-        h_A_primes.append(h_A_prime)
-
-    idx = np.argmin(abs(np.array([optimum.cost for optimum in optimums])))
-
-    return h_A_primes[idx], optimums[idx].fun, optimums[idx]
-
-
 def rel_err_between_abds(abd1, abd2, t1, t2):
     t_array = abd1.t[np.argmin(abs(abd1.t - t1)) : np.argmin(abs(abd1.t - t2)) + 1]
 
@@ -391,7 +213,7 @@ def map_to_abd_frame(
                 2.0 * abd_interp_superrest.sigma.bar, dataType=scri.h
             )
 
-            _, rel_err, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
+            rel_err, _, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
                 strain_interp_superrest,
                 target_strain_superrest,
                 t_0 - padding_time,
@@ -427,7 +249,7 @@ def map_to_abd_frame(
                 2.0 * abd_interp_prime.sigma.bar, dataType=scri.h
             )
 
-            _, rel_err, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
+            rel_err, _, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
                 strain_interp_prime,
                 target_strain,
                 t_0 - padding_time,

From 63ec1530ecb64ac986b017dbe67b2ec6ad3e0741 Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Tue, 4 Apr 2023 14:35:41 -0700
Subject: [PATCH 07/13] Enable writing/reading BMS transformtions to/from file

---
 scri/bms_transformations.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/scri/bms_transformations.py b/scri/bms_transformations.py
index 6bfbad11..167caa76 100644
--- a/scri/bms_transformations.py
+++ b/scri/bms_transformations.py
@@ -1,4 +1,5 @@
 import copy
+import h5py
 import scipy
 import spinsfast
 from scri.asymptotic_bondi_data.transformations import _process_transformation_kwargs, boosted_grid, conformal_factors
@@ -589,3 +590,34 @@ def __mul__(self, other):
         )
 
         return bms_composed
+
+    def to_file(self, filename):
+        dt = h5py.special_dtype(vlen=str)
+        with h5py.File(filename, "w") as hf:
+            hf.create_dataset("supertranslation", data=self.supertranslation)
+            hf.create_dataset("frame_rotation", data=self.frame_rotation.components)
+            hf.create_dataset("boost_velocity", data=self.boost_velocity)
+            hf.create_dataset("order", data=self.order)
+            hf.create_dataset("ell_max", data=self.ell_max)
+
+        return
+
+    def from_file(self, filename):
+        with h5py.File(filename, "r") as hf:
+            supertranslation = np.array(hf.get("supertranslation"))
+            frame_rotation = np.array(hf.get("frame_rotation"))
+            boost_velocity = np.array(hf.get("boost_velocity"))
+            order = [x.decode("utf-8") for x in np.array(hf.get("order"))]
+            ell_max = int(np.array(hf.get("ell_max")))
+
+        BMS = BMSTransformation(
+            frame_rotation=frame_rotation,
+            boost_velocity=boost_velocity,
+            supertranslation=supertranslation,
+            order=order,
+            ell_max=ell_max,
+        )
+
+        self.__dict__.update(BMS.__dict__)
+
+        return

From 1bed408192d23e8bccf1664728ec7571f8f0ab60 Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Fri, 14 Apr 2023 16:55:45 -0700
Subject: [PATCH 08/13] Fix angular velocity CubicSpline bug

---
 scri/asymptotic_bondi_data/map_to_superrest_frame.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/scri/asymptotic_bondi_data/map_to_superrest_frame.py b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
index bb5f23a2..69805d47 100644
--- a/scri/asymptotic_bondi_data/map_to_superrest_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
@@ -487,7 +487,7 @@ def rotation_to_map_to_superrest_frame(abd, target_strain=None, N_itr_max=10, re
 
         target_omega = target_news.angular_velocity()
         target_omega_spline = CubicSpline(
-            target_omega / np.linalg.norm(target_omega, axis=-1)[:, None], target_strain.t
+            target_strain.t, target_omega / np.linalg.norm(target_omega, axis=-1)[:, None]
         )
 
         itr = 0
@@ -637,7 +637,7 @@ def rel_err_for_abd_in_superrest(abd, target_PsiM, target_strain):
 
         target_omega = target_news.angular_velocity()
         target_omega_spline = CubicSpline(
-            target_omega / np.linalg.norm(target_omega, axis=-1)[:, None], target_strain.t
+            target_strain.t, target_omega / np.linalg.norm(target_omega, axis=-1)[:, None]
         )
 
         news = MT_to_WM(2.0 * abd.sigma.bar.dot, dataType=scri.hdot)

From fddcca30a77817c407d5185279c4e3891e13d448 Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Wed, 26 Apr 2023 15:19:06 -0700
Subject: [PATCH 09/13] fixes sxs compatibility issues

---
 .../asymptotic_bondi_data/map_to_abd_frame.py | 14 +++---
 .../map_to_superrest_frame.py                 | 43 ++++++++++++-------
 2 files changed, 35 insertions(+), 22 deletions(-)

diff --git a/scri/asymptotic_bondi_data/map_to_abd_frame.py b/scri/asymptotic_bondi_data/map_to_abd_frame.py
index 26c73b59..561bc45c 100644
--- a/scri/asymptotic_bondi_data/map_to_abd_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_abd_frame.py
@@ -16,7 +16,7 @@
 
 from functools import partial
 
-from scri.asymptotic_bondi_data.map_to_superrest_frame import time_translation, rotation
+from scri.asymptotic_bondi_data.map_to_superrest_frame import time_translation, rotation, MT_to_WM, WM_to_MT
 
 def rel_err_between_abds(abd1, abd2, t1, t2):
     t_array = abd1.t[np.argmin(abs(abd1.t - t1)) : np.argmin(abs(abd1.t - t2)) + 1]
@@ -213,9 +213,9 @@ def map_to_abd_frame(
                 2.0 * abd_interp_superrest.sigma.bar, dataType=scri.h
             )
 
-            rel_err, _, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
-                strain_interp_superrest,
-                target_strain_superrest,
+            rel_err, _, res = align2d(
+                MT_to_WM(WM_to_MT(strain_interp_superrest), sxs_version=True),
+                MT_to_WM(WM_to_MT(target_strain_superrest), sxs_version=True),
                 t_0 - padding_time,
                 t_0 + padding_time,
                 n_brute_force_δt=None,
@@ -249,9 +249,9 @@ def map_to_abd_frame(
                 2.0 * abd_interp_prime.sigma.bar, dataType=scri.h
             )
 
-            rel_err, _, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
-                strain_interp_prime,
-                target_strain,
+            rel_err, _, res = align2d(
+                MT_to_WM(WM_to_MT(strain_interp_prime), sxs_version=True),
+                MT_to_WM(WM_to_MT(target_strain), sxs_version=True),
                 t_0 - padding_time,
                 t_0 + padding_time,
                 n_brute_force_δt=None,
diff --git a/scri/asymptotic_bondi_data/map_to_superrest_frame.py b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
index 69805d47..1227021e 100644
--- a/scri/asymptotic_bondi_data/map_to_superrest_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
@@ -5,6 +5,8 @@
 import spinsfast
 import spherical_functions as sf
 from spherical_functions import LM_index
+
+from sxs.waveforms.alignment import align2d
 from scri.asymptotic_bondi_data.bms_charges import charge_vector_from_aspect
 
 import quaternion
@@ -291,7 +293,7 @@ def supertranslation_to_map_to_superrest_frame(
                     target_PsiM_Grid_interp[i, j] = CubicSpline(target_PsiM.t, target_PsiM_Grid[:, i, j])(alpha_i_j)
             M_Grid = -target_PsiM_Grid_interp.view(np.ndarray)
 
-            target_PsiM_interp = compute_Moreschi_supermomentum(target_PsiM, alpha_Grid, ell_max)
+            target_PsiM_interp = spinsfast.map2salm(target_PsiM_Grid_interp.view(np.ndarray), 0, ell_max)
 
             rel_err = np.linalg.norm(PsiM_interp[4:] - target_PsiM_interp[4:]) / np.linalg.norm(target_PsiM_interp[4:])
         else:
@@ -657,7 +659,7 @@ def rel_err_for_abd_in_superrest(abd, target_PsiM, target_strain):
     if target_PsiM is not None:
         rel_err_PsiM = np.linalg.norm(
             np.array(abd.supermomentum("Moreschi"))[np.argmin(abs(abd.t - 0)), 4:]
-            - np.array(target_PsiM)[np.argmin(abs(abd.t - 0)), 4:]
+            - np.array(target_PsiM.data)[np.argmin(abs(target_PsiM.t - 0)), 4:]
         )
     else:
         rel_err_PsiM = np.linalg.norm(np.array(abd.supermomentum("Moreschi"))[np.argmin(abs(abd.t - 0)), 4:])
@@ -778,7 +780,7 @@ def map_to_superrest_frame(
         alpha_ell_max = ell_max
 
     abd_interp = abd.interpolate(
-        abd.t[np.argmin(abs(abd.t - (t_0 - padding_time))) : np.argmin(abs(abd.t - (t_0 + padding_time))) + 1]
+        abd.t[np.argmin(abs(abd.t - (t_0 - (padding_time + 200)))) : np.argmin(abs(abd.t - (t_0 + (padding_time + 200)))) + 1]
     )
 
     # apply a time translation so that we're mapping
@@ -792,11 +794,17 @@ def map_to_superrest_frame(
     rel_err = [np.inf, np.inf, np.inf]
     rel_errs = [[np.inf, np.inf, np.inf]]
     best_rel_err = [np.inf, np.inf, np.inf]
-    while itr < N_itr_maxes["superrest"] and not (
-        rel_err[0] < rel_err_tols["CoM_transformation"]
-        and rel_err[1] < rel_err_tols["rotation"]
-        and rel_err[2] < rel_err_tols["supertranslation"]
-    ):
+    while itr < N_itr_maxes["superrest"]:
+        if type(rel_err) == tuple:
+            if (
+                    rel_err[0] < rel_err_tols["CoM_transformation"]
+                    and rel_err[1] < rel_err_tols["rotation"]
+                    and rel_err[2] < rel_err_tols["supertranslation"]
+            ):
+                break
+        else:
+            pass
+            
         if itr == 0:
             abd_interp_prime = abd_interp.transform(
                 supertranslation=BMS_transformation.supertranslation,
@@ -835,15 +843,15 @@ def map_to_superrest_frame(
                         2.0 * abd_interp_prime.sigma.bar, dataType=scri.h
                     )
 
-                    _, rel_err, res = scri.asymptotic_bondi_data.map_to_abd_frame.align2d(
-                        strain_interp_prime,
-                        target_strain,
-                        t_0 - padding_time,
-                        t_0 + padding_time,
+                    rel_err, _, res = align2d(
+                        MT_to_WM(WM_to_MT(strain_interp_prime), sxs_version=True),
+                        MT_to_WM(WM_to_MT(target_strain), sxs_version=True),
+                        0 - padding_time,
+                        0 + padding_time,
                         n_brute_force_δt=None,
                         n_brute_force_δϕ=None,
                         include_modes=None,
-                        nprocs=nprocs,
+                        nprocs=4,
                     )
 
                     new_transformation = scri.bms_transformations.BMSTransformation(
@@ -861,7 +869,12 @@ def map_to_superrest_frame(
                 boost_velocity=BMS_transformation.boost_velocity,
             )
 
-        rel_err = rel_err_for_abd_in_superrest(abd_interp_prime, target_PsiM, target_strain)
+        if target_strain is not None and order[-1] == "time_phase":
+            # rel_err is obtained from align2d, so do nothing
+            pass
+        else:
+            rel_err = rel_err_for_abd_in_superrest(abd_interp_prime, target_PsiM, target_strain)
+            
         if np.mean(rel_err) < min([np.mean(r) for r in rel_errs]):
             best_BMS_transformation = BMS_transformation.copy()
             best_rel_err = rel_err

From 00c1950b7dc1d48237b3c16bf6b33723226de8e0 Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Fri, 28 Apr 2023 12:11:18 -0700
Subject: [PATCH 10/13] Changes default order of superrest frame
 transformations

---
 scri/asymptotic_bondi_data/map_to_superrest_frame.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/scri/asymptotic_bondi_data/map_to_superrest_frame.py b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
index 1227021e..6d14a179 100644
--- a/scri/asymptotic_bondi_data/map_to_superrest_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
@@ -680,7 +680,7 @@ def map_to_superrest_frame(
         "supertranslation": 10,
     },
     rel_err_tols={"CoM_transformation": 1e-12, "rotation": 1e-12, "supertranslation": 1e-12},
-    order=["rotation", "CoM_transformation", "supertranslation"],
+    order=["supertranslation", "rotation", "CoM_transformation"],
     ell_max=None,
     alpha_ell_max=None,
     fix_time_phase_freedom=False,

From 0946610a764d18b817a88f5f2eefcac4dfdba45a Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Mon, 8 May 2023 12:05:06 -0700
Subject: [PATCH 11/13] Fix initial time translation guess

---
 scri/asymptotic_bondi_data/map_to_abd_frame.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/scri/asymptotic_bondi_data/map_to_abd_frame.py b/scri/asymptotic_bondi_data/map_to_abd_frame.py
index 561bc45c..eec80fa0 100644
--- a/scri/asymptotic_bondi_data/map_to_abd_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_abd_frame.py
@@ -153,8 +153,8 @@ def map_to_abd_frame(
         time_translation = scri.bms_transformations.BMSTransformation(
             supertranslation=[
                 sf.constant_as_ell_0_mode(
-                    abd.t[np.argmax(abd.bondi_four_momentum()[np.argmin(abs(abd.t - t_0)), 0])]
-                    - target_abd.t[np.argmax(target_abd.bondi_four_momentum()[np.argmin(abs(target_abd.t - t_0)), 0])]
+                    abd.t[np.argmax(abd.bondi_four_momentum()[:, 0])]
+                    - target_abd.t[np.argmax(target_abd.bondi_four_momentum()[:, 0])]
                 )
             ]
         )

From 2ae2bae7ba17e3ed8351b34e422894c3107fd3d6 Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Mon, 8 May 2023 16:30:07 -0700
Subject: [PATCH 12/13] Add modes option; fix time translation initial guess

---
 scri/asymptotic_bondi_data/map_to_abd_frame.py       | 9 +++++----
 scri/asymptotic_bondi_data/map_to_superrest_frame.py | 6 +++++-
 2 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/scri/asymptotic_bondi_data/map_to_abd_frame.py b/scri/asymptotic_bondi_data/map_to_abd_frame.py
index eec80fa0..e7951a68 100644
--- a/scri/asymptotic_bondi_data/map_to_abd_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_abd_frame.py
@@ -69,7 +69,7 @@ def map_to_abd_frame(
         "supertranslation": 10,
     },
     rel_err_tols={"CoM_transformation": 1e-12, "rotation": 1e-12, "supertranslation": 1e-12},
-    order=["rotation", "CoM_transformation", "supertranslation"],
+    order=["supertranslation", "rotation", "CoM_transformation"],
     ell_max=None,
     alpha_ell_max=None,
     fix_time_phase_freedom=True,
@@ -116,7 +116,7 @@ def map_to_abd_frame(
         }.
     order : list, optional
         Order in which to solve for the BMS transformations.
-        Default is ["rotation", "CoM_transformation", "supertranslation"].
+        Default is ["supertranslation", "rotation", "CoM_transformation"].
     ell_max : int, optional
         Maximum ell to use for SWSH/Grid transformations.
         Default is self.ell_max.
@@ -150,11 +150,12 @@ def map_to_abd_frame(
     time_translation = scri.bms_transformations.BMSTransformation()
     if fix_time_phase_freedom:
         # ensure that they are reasonable close
+        energy = abd.bondi_four_momentum()[:, 0]
+        target_energy = target_abd.bondi_four_momentum()[:, 0]
         time_translation = scri.bms_transformations.BMSTransformation(
             supertranslation=[
                 sf.constant_as_ell_0_mode(
-                    abd.t[np.argmax(abd.bondi_four_momentum()[:, 0])]
-                    - target_abd.t[np.argmax(target_abd.bondi_four_momentum()[:, 0])]
+                    abd.t[np.argmin(abs(energy - target_energy[np.argmin(abs(target_abd.t - t_0))]))] - t_0
                 )
             ]
         )
diff --git a/scri/asymptotic_bondi_data/map_to_superrest_frame.py b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
index 6d14a179..ada76e2e 100644
--- a/scri/asymptotic_bondi_data/map_to_superrest_frame.py
+++ b/scri/asymptotic_bondi_data/map_to_superrest_frame.py
@@ -684,6 +684,7 @@ def map_to_superrest_frame(
     ell_max=None,
     alpha_ell_max=None,
     fix_time_phase_freedom=False,
+    modes=None,
     print_conv=False,
 ):
     """Transform an abd object to the superrest frame.
@@ -747,6 +748,9 @@ def map_to_superrest_frame(
     fix_time_phase_freedom : bool, optional
         Whether or not to fix the time and phase freedom using a 2d minimization scheme.
         Default is True.
+    modes : list, optional
+        List of modes to include when performing the 2d alignment.
+        Default is every mode.
     print_conv: bool, defaults to False
         Whether or not to print the termination criterion. Default is False.
 
@@ -850,7 +854,7 @@ def map_to_superrest_frame(
                         0 + padding_time,
                         n_brute_force_δt=None,
                         n_brute_force_δϕ=None,
-                        include_modes=None,
+                        include_modes=modes,
                         nprocs=4,
                     )
 

From 5541dd9669b62172525dbbc6891301ec79c4dcbb Mon Sep 17 00:00:00 2001
From: Keefe Mitman <kmitman@caltech.edu>
Date: Mon, 8 May 2023 17:43:55 -0700
Subject: [PATCH 13/13] Fix abd to abd test

---
 tests/test_abd_frame.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_abd_frame.py b/tests/test_abd_frame.py
index 776f67e5..6353e017 100644
--- a/tests/test_abd_frame.py
+++ b/tests/test_abd_frame.py
@@ -45,7 +45,7 @@ def test_abd_to_abd():
         abd_target,
         t_0=0,
         padding_time=20,
-        N_itr_maxes={"abd": 1, "superrest": 1, "CoM_transformation": 10, "rotation": 10, "supertranslation": 10},
+        N_itr_maxes={"abd": 2, "superrest": 1, "CoM_transformation": 10, "rotation": 10, "supertranslation": 10},
         fix_time_phase_freedom=False,
         nprocs=-1,
     )