Skip to content

Fix x-axis during superrest rotation and minor changes #104

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
moble merged 3 commits into from
Jul 14, 2026
Merged

Fix x-axis during superrest rotation and minor changes #104

Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5ff39cc
Fix x-axis during rotation, remove unused 'fix_time_phase' option, ac…
keefemitman Apr 4, 2025
77a51e8
Update code to fix x vector to xz plane and update docstrings
keefemitman Apr 8, 2025
11f16c8
Add fix yz plane option
keefemitman Jul 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 11 additions & 7 deletions scri/asymptotic_bondi_data/map_to_superrest_frame.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -415,8 +415,13 @@ def rotation_from_spin_charge(chi, t):
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)
q = (1 - chi_f * quaternion.z).normalized()

x_rotated = quaternion.as_vector_part(q.inverse() * quaternion.x * q)
phase = np.angle(x_rotated[0] + 1j*x_rotated[1])
q = q * quaternion.from_rotation_vector(phase * np.array([0, 0, 1]))

@moble moble Apr 4, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand this.

The first definition of q is the (smallest) rotation that takes z directly onto chi_f. This should behave well, in the sense of making a very small change to x and y, when chi_f is close to z — that is, there is no coordinate singularity there. It should behave badly if chi_f is close to -z, and fail if chi_f=-z, but I think that's intrinsic to what you're trying to do. So I think that part is entirely natural.

But I don't see what the rest of it is doing. You do have the freedom to prepend any rotation about the z axis and still rotate z onto chi_ref. But why this phase? It's impossible to actually fix the x axis for general chi_ref. I can imagine trying to confine it to the x-y or x-z plane, for example, but this doesn't do that. What's the goal?

@keefemitman keefemitman Apr 4, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well this makes it so that chi_f is always point strictly along z, and the addition of lines 420 - 422 make it so that the original x vector, after the rotation, has no component in the y direction. So I thought this was confining the x-vector to be in the x-z- plane, and certainly from my tests this seems to be the case.

Do you disagree?

I found this to be useful when trying to preserve phase information while mapping to the remnant frame and extracting QNMs, as it makes the QNM phases more consistent as a function of progenitors.

@moble moble Apr 4, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here's a simple test case:

import numpy as np
import quaternion
chi_f = quaternion.from_vector_part([1,1,0]).normalized()

q = (1 - chi_f * quaternion.z).normalized()

x_rotated = quaternion.as_vector_part(q.inverse() * quaternion.x * q)
phase = np.angle(x_rotated[0] + 1j*x_rotated[1])
q = q * quaternion.from_rotation_vector(phase * np.array([0, 0, 1]))

(q * quaternion.x * q.inverse()).y

Which results in -0.7071067811865475, so it's definitely not in the x-z plane. I'm sure there's a simple fix to make this work, but I don't have the time or brainpower right now to think about this.

Also, I'm not convinced this is a good idea. When you impose artificial conditions like that on a rotation, you often get terrible numerical behavior. I guess it would be a small effect anyway when chi_f starts out close to z, which would be fine. But when you're getting chi_f closer to the diagonals between x and y, this starts introducing rotations of 90°, and then you'll start getting sign flips. So I'm very concerned about the robustness of this sort of thing.

@keefemitman keefemitman Apr 4, 2025 •
edited
Loading

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, don't you have your q and q.inverse() mixed up? If you replace quaternion.x in (q * quaternion.x * q.inverse()).y with chi_f you don't get z. If you do (q.inverse() * quaternion.x * q).y then it works.

@moble moble Apr 7, 2025 •
edited
Loading

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm. Now I'm confused about what this function is supposed to return. The return statement identifies q as a frame_rotation, which is defined here. The input system is decomposed in a frame (x, y, z), but its spin charge is off in some different direction chi_f; we want to transform to a new frame where z' = chi_f, so that definition says that

chi_f = z' = q * z * q.inverse()

which is what you get with this q (both before and after the new stuff).

Of course, the docstring does seem more consistent with the opposite interpretation. But in that case, I would have thought the return statement should use frame_rotation=q.inverse().components.

@keefemitman keefemitman Apr 7, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's suppose to return the rotation (as a BMSTransformation) such that abd.transform(BMSTransformation.frame_rotation.components) yields an abd object whose spin vector is pointing along the z axis (which is what the code certainly does as of now). Is this confusion just stemming from how the transform function does a passive rather than an active transformation?

@moble moble Apr 7, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's suppose to return the rotation (as a BMSTransformation) such that abd.transform(BMSTransformation.frame_rotation.components) yields an abd object whose spin vector is pointing along the z axis (which is what the code certainly does as of now).

Where that last z is really z_prime, right? Meaning that if you do

abd_prime = abd.transform(BMSTransformation.frame_rotation.components)
chi_prime = abd_prime.bondi_dimensionless_spin()[np.argmin(abs(t))]

then chi_prime[0] and chi_prime[1] are both approximately 0.0, right? If so, that's consistent with my interpretation.

So the transformation will take the field and decompose it in a new coordinate system, the basis vectors of which are related to the basis vectors of the coordinate system of the original abd (which are the same as those used in rotation_from_spin_charge) by

x_prime = q * quaternion.x * q.inverse()
y_prime = q * quaternion.y * q.inverse()
z_prime = q * quaternion.z * q.inverse()

and

z_prime = chi_f

But x_prime.y is not zero.

@keefemitman keefemitman Apr 7, 2025 •
edited
Loading

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah ok I understand now; I was using the inverted transformation formula when working this out, so in reality I was adjusting the quaternion such that y_prime.x was zero. If I just change the new code to be

y_rotated = quaternion.as_vector_part(q * quaternion.y * q.inverse())
phase = np.angle(1j*(y_rotated[0] + 1j*y_rotated[1]))
q = q * quaternion.from_rotation_vector(phase * np.array([0, 0, 1]))

then it will really make x_prime.y be zero. But it sounds like you don't want to do this anyways? So should I make this change or do you just want to just scrap this all together?

@moble moble Apr 8, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, that makes sense. If you and/or Aniket still like it, I have no objection to making this an option, with the default being the old behavior.

But could you also update the docstring with our understanding of what's supposed to be returned, including the statements that

chi(t) = q * z * q.inverse()

and if the option is True, then

(q * x * q.inverse()).y = 0

@keefemitman keefemitman Apr 8, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok sounds good! Just pushed an update.


return scri.bms_transformations.BMSTransformation(frame_rotation=q.components)


def rotation_from_vectors(vector, target_vector, t=None):
Expand Down Expand Up @@ -683,7 +688,6 @@ def map_to_superrest_frame(
order=["supertranslation", "rotation", "CoM_transformation"],
ell_max=None,
alpha_ell_max=None,
fix_time_phase_freedom=False,
modes=None,
print_conv=False,
):
Expand Down Expand Up @@ -739,15 +743,13 @@ def map_to_superrest_frame(
order : list, optional
Order in which to solve for the BMS transformations.
Default is ["rotation", "CoM_transformation", "supertranslation"].
If "time_phase" is included, then a time/phase optimization is performed.
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.
modes : list, optional
List of modes to include when performing the 2d alignment.
Default is every mode.
Expand All @@ -765,6 +767,9 @@ def map_to_superrest_frame(
"""
abd = self.copy()

if order == []:
return abd, scri.bms_transformations.BMSTransformation(), None

if target_strain_input is not None:
target_strain = target_strain_input.copy()
target_strain.t -= t_0
Expand Down Expand Up @@ -874,7 +879,6 @@ def map_to_superrest_frame(
)

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)
Expand Down