From db3e89c2ecbf1b3946d5586be0d994a343cec8b8 Mon Sep 17 00:00:00 2001 From: Andrew Tribick Date: Wed, 25 Oct 2023 00:26:35 +0200 Subject: [PATCH] Add quaternion constants in celmath, use sincos in rotation matrix functions --- src/celengine/axisarrow.cpp | 9 +++--- src/celengine/observer.cpp | 2 +- src/celengine/planetgrid.cpp | 4 +-- src/celengine/skygrid.cpp | 18 ++++++----- src/celephem/customorbit.cpp | 4 +-- src/celephem/customrotation.cpp | 6 ++-- src/celephem/samporient.cpp | 8 ++--- src/celephem/spicerotation.cpp | 7 +++-- src/celestia/gtk/dialog-eclipse.cpp | 2 +- src/celestia/win32/wineclipses.cpp | 2 +- src/celmath/geomutil.h | 47 +++++++++++++++++++++++------ 11 files changed, 70 insertions(+), 39 deletions(-) diff --git a/src/celengine/axisarrow.cpp b/src/celengine/axisarrow.cpp index 27111d54622..fc82a414b71 100644 --- a/src/celengine/axisarrow.cpp +++ b/src/celengine/axisarrow.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -308,19 +309,19 @@ AxesReferenceMark::render(Renderer* renderer, Eigen::Matrix4f labelTransformMatrix = labelTransform.matrix(); // x-axis - Eigen::Matrix4f xModelView = modelView * celmath::rotate(Eigen::AngleAxisf(90.0_deg, Eigen::Vector3f::UnitY())); + Eigen::Matrix4f xModelView = modelView * celmath::rotate(celmath::YRot90); glVertexAttrib4f(CelestiaGLProgram::ColorAttributeIndex, 1.0f, 0.0f, 0.0f, opacity); prog->setMVPMatrices(projection, xModelView); GetArrowVAO().draw(); // y-axis - Eigen::Matrix4f yModelView = modelView * celmath::rotate(Eigen::AngleAxisf(180.0_deg, Eigen::Vector3f::UnitY())); + Eigen::Matrix4f yModelView = modelView * celmath::rotate(celmath::YRot180); glVertexAttrib4f(CelestiaGLProgram::ColorAttributeIndex, 0.0f, 1.0f, 0.0f, opacity); prog->setMVPMatrices(projection, yModelView); GetArrowVAO().draw(); // z-axis - Eigen::Matrix4f zModelView = modelView * celmath::rotate(Eigen::AngleAxisf(-90.0_deg, Eigen::Vector3f::UnitX())); + Eigen::Matrix4f zModelView = modelView * celmath::rotate(celmath::XRot90Conjugate); glVertexAttrib4f(CelestiaGLProgram::ColorAttributeIndex, 0.0f, 0.0f, 1.0f, opacity); prog->setMVPMatrices(projection, zModelView); GetArrowVAO().draw(); @@ -455,7 +456,7 @@ BodyAxisArrows::BodyAxisArrows(const Body& _body) : Eigen::Quaterniond BodyAxisArrows::getOrientation(double tdb) const { - return (Eigen::Quaterniond(Eigen::AngleAxis(celestia::numbers::pi, Eigen::Vector3d::UnitY())) * body.getEclipticToBodyFixed(tdb)).conjugate(); + return (celmath::YRot180 * body.getEclipticToBodyFixed(tdb)).conjugate(); } diff --git a/src/celengine/observer.cpp b/src/celengine/observer.cpp index 6d9b3ecffe6..d8a045ce0e1 100644 --- a/src/celengine/observer.cpp +++ b/src/celengine/observer.cpp @@ -538,7 +538,7 @@ void Observer::setLocationFilter(uint64_t _locationFilter) void Observer::reverseOrientation() { - setOrientation(getOrientation() * Quaterniond(AngleAxisd(celestia::numbers::pi, Vector3d::UnitY()))); + setOrientation(getOrientation() * celmath::YRot180); reverseFlag = !reverseFlag; } diff --git a/src/celengine/planetgrid.cpp b/src/celengine/planetgrid.cpp index b84f6565d18..0d2ffb2ae72 100644 --- a/src/celengine/planetgrid.cpp +++ b/src/celengine/planetgrid.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -99,8 +100,7 @@ PlanetographicGrid::render(Renderer* renderer, InitializeGeometry(*renderer); // Compatibility - Eigen::Quaterniond q(Eigen::AngleAxis(celestia::numbers::pi, Eigen::Vector3d::UnitY())); - q *= body.getEclipticToBodyFixed(tdb); + Eigen::Quaterniond q = celmath::YRot180 * body.getEclipticToBodyFixed(tdb); Eigen::Quaternionf qf = q.cast(); // The grid can't be rendered exactly on the planet sphere, or diff --git a/src/celengine/skygrid.cpp b/src/celengine/skygrid.cpp index 85166602b15..0e0a507163a 100644 --- a/src/celengine/skygrid.cpp +++ b/src/celengine/skygrid.cpp @@ -335,10 +335,6 @@ SkyGrid::render(Renderer& renderer, int windowWidth, int windowHeight) { - // 90 degree rotation about the x-axis used to transform coordinates - // to Celestia's system. - Quaterniond xrot90 = XRotation(-celestia::numbers::pi / 2.0); - auto vfov = static_cast(renderer.getProjectionMode()->getFOV(observer.getZoom())); double viewAspectRatio = static_cast(windowWidth) / static_cast(windowHeight); @@ -369,7 +365,13 @@ SkyGrid::render(Renderer& renderer, Vector3d c3( w, h, -1.0); Quaterniond cameraOrientation = renderer.getCameraOrientation(); - Matrix3d r = (cameraOrientation * xrot90 * m_orientation.conjugate() * xrot90.conjugate()).toRotationMatrix().transpose(); + + // 90 degree rotation about the x-axis used to transform coordinates + // to Celestia's system. + Matrix3d r = (cameraOrientation * + celmath::XRot90Conjugate * + m_orientation.conjugate() * + celmath::XRot90).toRotationMatrix().transpose(); // Transform the frustum corners by the camera and grid // rotations. @@ -481,13 +483,15 @@ SkyGrid::render(Renderer& renderer, int endDec = (int) std::floor(DEG_MIN_SEC_TOTAL * (maxDec / celestia::numbers::pi) / (double) decIncrement) * decIncrement; // Get the orientation at single precision - Quaterniond q = xrot90 * m_orientation * xrot90.conjugate(); + Quaterniond q = celmath::XRot90Conjugate * m_orientation * celmath::XRot90; Quaternionf orientationf = q.cast(); // Radius of sphere is arbitrary, with the constraint that it shouldn't // intersect the near or far plane of the view frustum. Matrix4f m = renderer.getModelViewMatrix() * - celmath::rotate((xrot90 * m_orientation.conjugate() * xrot90.conjugate()).cast()) * + celmath::rotate((celmath::XRot90Conjugate * + m_orientation.conjugate() * + celmath::XRot90).cast()) * celmath::scale(1000.0f); Matrices matrices = {&renderer.getProjectionMatrix(), &m}; diff --git a/src/celephem/customorbit.cpp b/src/celephem/customorbit.cpp index f5ae9622088..239896be6f4 100644 --- a/src/celephem/customorbit.cpp +++ b/src/celephem/customorbit.cpp @@ -1366,8 +1366,8 @@ class PhobosOrbit : public CachingOrbit // Orientation of the orbital plane with respect to the Laplacian plane Eigen::Matrix3d Rorbit = (celmath::YRotation(node) * - celmath::XRotation(celmath::degToRad(i)) * - celmath::YRotation(w)).toRotationMatrix(); + celmath::XRotation(celmath::degToRad(i)) * + celmath::YRotation(w)).toRotationMatrix(); // Rotate to the Earth's equatorial plane constexpr double N = celmath::degToRad(refplane_RA); diff --git a/src/celephem/customrotation.cpp b/src/celephem/customrotation.cpp index 07132a5c4fe..aca1c0684e0 100644 --- a/src/celephem/customrotation.cpp +++ b/src/celephem/customrotation.cpp @@ -90,7 +90,7 @@ class IAURotationModel : public CachingRotationModel double inclination = 90.0 - poleDec; if (flipped) - return celmath::XRotation(celestia::numbers::pi) * + return celmath::XRot180 * celmath::XRotation(celmath::degToRad(-inclination)) * celmath::YRotation(celmath::degToRad(-node)); else @@ -176,9 +176,7 @@ class EarthRotationModel : public CachingRotationModel Eigen::Quaterniond q = celmath::XRotation(obliquity) * celmath::ZRotation(-precession) * eclRotation.conjugate(); // convert to Celestia's coordinate system - return celmath::XRotation(celestia::numbers::pi / 2.0) * - q * - celmath::XRotation(-celestia::numbers::pi / 2.0); + return celmath::XRot90 * q * celmath::XRot90Conjugate; } double getPeriod() const override diff --git a/src/celephem/samporient.cpp b/src/celephem/samporient.cpp index 9a6c10d4494..e343f7c714c 100644 --- a/src/celephem/samporient.cpp +++ b/src/celephem/samporient.cpp @@ -111,14 +111,12 @@ class SampledOrientation : public RotationModel void SampledOrientation::addSample(double t, const Eigen::Quaternionf& q) { - // 90 degree rotation about x-axis to convert orientation to Celestia's - // coordinate system. - static const Eigen::Quaternionf coordSysCorrection = celmath::XRotation((float) (celestia::numbers::pi / 2.0)); - // TODO: add a check for out of sequence samples OrientationSample& samp = samples.emplace_back(); samp.t = t; - samp.q = q * coordSysCorrection; + // 90 degree rotation about x-axis to convert orientation to Celestia's + // coordinate system. + samp.q = q * celmath::XRot90; } diff --git a/src/celephem/spicerotation.cpp b/src/celephem/spicerotation.cpp index b757cd59805..70cddb538aa 100644 --- a/src/celephem/spicerotation.cpp +++ b/src/celephem/spicerotation.cpp @@ -182,9 +182,10 @@ SpiceRotation::computeSpin(double jd) const Eigen::Quaterniond q = Eigen::Quaterniond(Eigen::Map(matrixData)).conjugate(); // Transform into Celestia's coordinate system - static const Eigen::Quaterniond Rx90 = celmath::XRotation(celestia::numbers::pi / 2.0); - static const Eigen::Quaterniond Ry180 = celmath::YRotation(celestia::numbers::pi); - return Ry180 * Rx90.conjugate() * q.conjugate() * Rx90; + return celmath::YRot180 * + celmath::XRot90Conjugate * + q.conjugate() * + celmath::XRot90; } } diff --git a/src/celestia/gtk/dialog-eclipse.cpp b/src/celestia/gtk/dialog-eclipse.cpp index 84d812f1211..f8963597b89 100644 --- a/src/celestia/gtk/dialog-eclipse.cpp +++ b/src/celestia/gtk/dialog-eclipse.cpp @@ -332,7 +332,7 @@ static gint eclipseGoto(GtkButton*, EclipseData* ed) double distance = target.radius() * 4.0; sim->gotoLocation(UniversalCoord::Zero().offsetKm(Vector3d::UnitX() * distance), - (YRotation(-celestia::numbers::pi / 2) * XRotation(-celestia::numbers::pi / 2)), 2.5); + (celmath::YRot90Conjugate * celmath::XRot90Conjugate), 2.5); return TRUE; } diff --git a/src/celestia/win32/wineclipses.cpp b/src/celestia/win32/wineclipses.cpp index 94f632ba34c..d72f181db97 100644 --- a/src/celestia/win32/wineclipses.cpp +++ b/src/celestia/win32/wineclipses.cpp @@ -390,7 +390,7 @@ BOOL APIENTRY EclipseFinderProc(HWND hDlg, double distance = target.radius() * 4.0; sim->gotoLocation(UniversalCoord::Zero().offsetKm(Vector3d::UnitX() * distance), - YRotation(-celestia::numbers::pi / 2) * XRotation(-celestia::numbers::pi / 2), + celmath::YRot90Conjugate * celmath::XRot90Conjugate, 2.5); } break; diff --git a/src/celmath/geomutil.h b/src/celmath/geomutil.h index d604c7627a3..7a762e41896 100644 --- a/src/celmath/geomutil.h +++ b/src/celmath/geomutil.h @@ -15,35 +15,64 @@ #include #include +#include #include "mathlib.h" namespace celmath { +template const inline Eigen::Quaternion XRot90{ celestia::numbers::sqrt2_v * T{0.5}, + celestia::numbers::sqrt2_v * T{0.5}, + T{0}, + T{0} }; + +template const inline Eigen::Quaternion XRot180{ 0, 1, 0, 0 }; + +// Conjugate of 90 degree rotation = -90 degree rotation +template const inline Eigen::Quaternion XRot90Conjugate{ celestia::numbers::sqrt2_v * T{0.5}, + -(celestia::numbers::sqrt2_v * T{0.5}), + T{0}, + T{0} }; + +template const inline Eigen::Quaternion YRot90{ celestia::numbers::sqrt2_v * T{0.5}, + T{0}, + celestia::numbers::sqrt2_v * T{0.5}, + T{0} }; + +template const inline Eigen::Quaternion YRot180{ 0, 0, 1, 0 }; + +template const inline Eigen::Quaternion YRot90Conjugate{ celestia::numbers::sqrt2_v * T{0.5}, + T{0}, + -(celestia::numbers::sqrt2_v * T{0.5}), + T{0} }; + template Eigen::Quaternion XRotation(T radians) { - using std::cos, std::sin; - T halfAngle = radians * static_cast(0.5); - return Eigen::Quaternion(cos(halfAngle), sin(halfAngle), 0, 0); + T s; + T c; + sincos(radians * static_cast(0.5), s, c); + return Eigen::Quaternion(c, s, 0, 0); } template Eigen::Quaternion YRotation(T radians) { - using std::cos, std::sin; - T halfAngle = radians * static_cast(0.5); - return Eigen::Quaternion(cos(halfAngle), 0, sin(halfAngle), 0); + T s; + T c; + sincos(radians * static_cast(0.5), s, c); + return Eigen::Quaternion(c, 0, s, 0); } template Eigen::Quaternion ZRotation(T radians) { - using std::cos, std::sin; - T halfAngle = radians * static_cast(0.5); - return Eigen::Quaternion(cos(halfAngle), 0, 0, sin(halfAngle)); + T s; + T c; + sincos(radians * static_cast(0.5), s, c); + return Eigen::Quaternion(c, 0, 0, s); }