Calculating principle axes from inertia tensor.

[tomazos]

Further to recent discussions about inertia tensor, I've written a C++ function that given an inertia tensor (actually any symmetric matrix), returns the three eigenvalues in order of size. Feel free to use it in RP3D, I've included it below. I haven't yet figured out how to get the principle axes from the eigenvalues. I think they are a non-zero member of the nullspace of inertia_tensor - eigenvalue * identity_matrix but my linear algebra isn't quite good enough yet to figure out how to actually code that. The wikipedia page says:

Once again, the eigenvectors of A can be obtained by recourse to the Cayley–Hamilton theorem. If α1, α2, α3 are distinct eigenvalues of A, then (A - α1I)(A - α2I)(A - α3I) = 0. Thus the columns of the product of any two of these matrices will contain an eigenvector for the third eigenvalue. However, if α3 = α1, then (A - α1I)^2(A - α2I) = 0 and (A - α2I)(A - α1I)^2 = 0. Thus the generalized eigenspace of α1 is spanned by the columns of A - α2I while the ordinary eigenspace is spanned by the columns of (A - α1I)(A - α2I). The ordinary eigenspace of α2 is spanned by the columns of (A - α1I)^2.

But how that translates into code is still slightly beyond my grasp. Maybe you have a tip?

template <typename T, qualifier Q>
GLM_FUNC_QUALIFIER vec<3, T, Q> EigenvaluesSymmetric(mat<3, 3, T, Q> const& A) {
  T eig[3];
  auto sqr = [](T x) -> T { return x * x; };

  T p1 = sqr(A[0][1]) + sqr(A[0][2]) + sqr(A[1][2]);
  if (p1 == 0) {
    eig[0] = A[0][0];
    eig[1] = A[1][1];
    eig[2] = A[2][2];
    if (eig[0] < eig[1]) std::swap(eig[0], eig[1]);
    if (eig[1] < eig[2]) std::swap(eig[1], eig[2]);
    if (eig[0] < eig[1]) std::swap(eig[0], eig[1]);

    return vec<3, T, Q>(eig[0], eig[1], eig[2]);
  }

  T q = (A[0][0] + A[1][1] + A[2][2]) / T(3);

  T p2 = sqr(A[0][0] - q) + sqr(A[1][1] - q) + sqr(A[2][2] - q) + T(2) * p1;

  T p = std::sqrt(p2 / T(6));

  mat<3, 3, T, Q> B = (T(1) / p) * (A - q * mat<3, 3, T, Q>(1.0));

  T r = determinant(B) / T(2);

  T phi;
  if (r <= -T(1))
    phi = T(M_PI) / T(3);
  else if (r >= T(1))
    phi = 0;
  else
    phi = std::acos(r) / T(3);

  eig[0] = q + T(2) * p * std::cos(phi);
  eig[2] = q + T(2) * p * std::cos(phi + T(2) * M_PI / T(3));
  eig[1] = T(3) * q - eig[0] - eig[2];
  return vec<3, T, Q>(eig[0], eig[1], eig[2]);
}

[tomazos]

https://math.stackexchange.com/q/4024159/22933

[tomazos]

https://math.stackexchange.com/q/4024266/22933