Reputation: 61
I have just started trying to wrap my head around the basics of using quaternions in my rocket simulation to handle guidance maneuvers. In my sim, I keep an ECI-to-body quaternion transform updated at each time step. The body frame that I am using is a vehicle-carried NED frame. The positive x-axis points out the nose, the positive z-axis points down, and the positive y-axis completes the right hand rule. To make things simple, I launch the rocket vertically at a latitude of 0 deg and longitude of 0 deg. At launch, the ECI and body frames are exactly aligned. So at t=0, the rocket is pointed down the x-axis of the ECI frame, with the body frame's positive z-axis (which points out from the 'belly' of the rocket) pointing in the positive z-axis of the ECI frame. I use a quaternion of the form [w,x,y,z] that, when multiplied with a vector, rotates that vector and returns that rotated vector in the same frame. I believe this is called an active quaternion.
Lets say that a positive rotation about the body frame's y-axis of 30 degrees is performed. Updating the ECI-to-body quaternion to represent this rotation gives [0.965, 0, 0.258, 0]. After this rotation is applied, i then want to express the rockets's body axes in the ECI frame. How do i do this? I originally thought that i should compute the inverse of the ECI-to-body transform, which would give me a body-to-ECI transform. I thought that i could just multiply this body-to-ECI transform with the 3 body axis vectors ([1,0,0], [0,1,0], and [0,0,1]) to give me those three vectors expressed in the ECI frame. For example, if I compute the inverse of this ECI-to-body transform, i get [0.965, 0, -0.258, 0]. If i then multiply this body-to-ECI transform with the body vector [1,0,0], I get [0.866, 0, 0.5]. However, the z-component of this vector is positive, even though a positive, right-handed rotation about the body y-axis was performed. This means that the rocket should be pointed slightly in the negative z-direction in the ECI frame after the rotation. What am I doing wrong?
I notice that if I don't use the inverse (body-to-ECI), and instead use the ECI-to-body transform to multiply with the body frame basis vectors, I get the correct result. Is this the correct way of performing this computation? If I have thrust and aero body force vectors in the body frame, do I just multiply them with the ECI-to-body transform to express them in the ECI frame? That is very non-intuitive for me. Below is some Python code that details my problem:
import numpy as np
def normalize(q):
"""Normalize a quaternion to unit length."""
return q / np.linalg.norm(q)
def quaternion_conjugate(q):
"""Compute the conjugate of a quaternion."""
w, x, y, z = q
return np.array([w, -x, -y, -z])
def quaternion_multiply(q1, q2):
"""Multiply two quaternions."""
w1, x1, y1, z1 = q1
w2, x2, y2, z2 = q2
return np.array([
w1*w2 - x1*x2 - y1*y2 - z1*z2,
w1*x2 + x1*w2 + y1*z2 - z1*y2,
w1*y2 - x1*z2 + y1*w2 + z1*x2,
w1*z2 + x1*y2 - y1*x2 + z1*w2
])
def rotate_vector(q, v):
"""Rotate a 3D vector v using quaternion q."""
q = normalize(q)
q_conj = quaternion_conjugate(q)
q_v = np.hstack(([0], v))
q_rotated = quaternion_multiply(quaternion_multiply(q, q_v), q_conj)
return q_rotated[1:]
def rotate_about_y(q, theta):
"""Rotate quaternion q by theta (radians) about the y-axis."""
half_theta = theta / 2
q_rot = np.array([np.cos(half_theta), 0, np.sin(half_theta), 0])
return normalize(quaternion_multiply(q_rot, q))
# ECI and body frame initially align
eci_to_body = np.array([1,0,0,0])
# Right-handed rotation of 30 degrees about body frame +y-axis
theta = np.radians(30)
eci_to_body = rotate_about_y(eci_to_body, theta)
body_to_eci = quaternion_conjugate(eci_to_body)
print(f"eci_to_body = {eci_to_body}") # [0.965, 0, 0.258, 0]
print(f"body_to_eci = {body_to_eci}") # [0.965, 0, -0.258, 0]
# Body frame x-axis (out the nose)
pointing_vector_body = np.array([1,0,0])
# Express the body frame x-axis in the ECI frame using body_to_eci. This
# results in the pointing vector having the wrong sign for the z-component
pointing_vector_eci_wrong = rotate_vector(body_to_eci, pointing_vector_body)
print(f"pointing_vector_eci_wrong = {pointing_vector_eci_wrong}") # [0.866, 0, 0.5]
# Express the body frame x-axis in the ECI frame using eci_to_body. This
# gives the correct answer. Is this by accident, or is eci_to_body*body_vector
# the correct way of transforming body frame vectors to the ECI frame?
pointing_vector_eci_correct = rotate_vector(eci_to_body, pointing_vector_body)
print(f"pointing_vector_eci_correct = {pointing_vector_eci_correct}") # [0.866, 0, -0.5]
Upvotes: 0
Views: 66
Reputation: 2636
You make this statement:
"I use a quaternion of the form [w,x,y,z] that, when multiplied with a vector, rotates that vector and returns that rotated vector in the same frame. I believe this is called an active quaternion."
But then you describe using your quaternions to change vector coordinate frames rather than rotating a vector within the same coordinate frame. These are two different things. You can use active convention to represent coordinate frame transformations (but not vice-versa), however it is more typical to use passive quaternions for this, at least in the aerospace industry where I work. I.e., your coordinate transformation would look like this instead using Right-Chain convention:
v_BODY = q_ECI2BODY^-1 * v_ECI * q_ECI2BODY
and successive rotations would stack up on the right instead of the left. E.g., suppose you wanted a new body attitude then the stacking would be:
V_NEWBODY = q_BODY2NEWBODY^-1 * q_ECI2BODY^-1 * v_ECI * q_ECI2BODY * q_BODY2NEWBODY
You might want to try this change in your code to see if it matches your intuition.
Btw, I haven't checked your quaternion multiplication, but are you using Hamilton Right-Handed or JPL Left-Handed? I try to always use Hamilton Right-Handed myself, but I know that some software such as Unity uses JPL Left-Handed. (Why???)
EDIT: I just checked and yes, you are using Hamilton Right-Handed convention.
Upvotes: 0