Fix x-axis during superrest rotation and minor changes

[keefemitman]

In addition, remove unused 'fix_time_phase' option and return the input when order=[].

Needed for PR #103.

[moble]

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]

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]

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]

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]

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]

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]

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]

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]

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]

Ok sounds good! Just pushed an update.

[duetosymmetry]

I don't see why this is needed for #103?

[keefemitman]

I thought Aniket wanted this? (specifically the removing fix_time_phase)

[akhairna]

Sorry for the late response. I needed the time_phase freedom fixing to perform some testing in my work. This isn't actually needed for the previous PR. But I wanted to rebase my #101 and #103 on top of #104 before getting mine merged. @moble, @duetosymmetry - I think once #104 is merged, then #101 and #103 can be merged sequentially. If there is a better suggestion for me, please let me know.

[moble]

Something's gone wrong here: if chi = z, I would expect to get q=1 at the end, but this gives me q=z.

I think it should be phase = np.atan2(-y_rotated[0], y_rotated[1]).

[keefemitman]

Ah yeah good catch; np.angle(-1j*(y_rotated[0] + 1j*y_rotated[1])) works, but using atan2 is fine by me!