From ae3e5c201b5e7ff7012bf23eca2ffc262932ba6a Mon Sep 17 00:00:00 2001 From: Cat Flynn Date: Wed, 21 Aug 2024 14:51:37 +0100 Subject: [PATCH] feat: use doubles in lib --- lib/skein/include/skein/orbit.h | 42 +++++++------ lib/skein/include/skein/particlemap.h | 2 +- lib/skein/src/orbit.cpp | 90 +++++++++++++-------------- lib/skein/src/particlemap.cpp | 26 +++++--- 4 files changed, 87 insertions(+), 73 deletions(-) diff --git a/lib/skein/include/skein/orbit.h b/lib/skein/include/skein/orbit.h index dd16ed5..4b679af 100644 --- a/lib/skein/include/skein/orbit.h +++ b/lib/skein/include/skein/orbit.h @@ -4,7 +4,7 @@ #include "glm/glm.hpp" -typedef std::vector Vector6; +typedef std::vector Vector6; class Orbit { @@ -12,28 +12,34 @@ public: Orbit(); ~Orbit() = default; - float getSemiMajorAxis() const; - float getEccentricity() const; - float getInclination() const; - float getArgumentOfPeriapsis() const; - float getLongitudeOfAscendingNode() const; + double getSemiMajorAxis() const; + double getEccentricity() const; + double getInclination() const; + double getArgumentOfPeriapsis() const; + double getLongitudeOfAscendingNode() const; - void setSemiMajorAxis(float semiMajorAxis); - void setEccentricity(float eccentricity); - void setInclination(float inclination); - void setArgumentOfPeriapsis(float argumentOfPeriapsis); - void setLongitudeOfAscendingNode(float longitudeOfAscendingNode); + void setSemiMajorAxis(double semiMajorAxis); + void setEccentricity(double eccentricity); + void setInclination(double inclination); + void setArgumentOfPeriapsis(double argumentOfPeriapsis); + void setLongitudeOfAscendingNode(double longitudeOfAscendingNode); - // TODO: meanAnomaly in all these arguments actually means eccentricMeanAnomaly, + // TODO: meanAnomaly in all these arguments actually means ellipticalMeanAnomaly, // will have to change that when adding non-ellipctical orbits - don't get confused! - glm::vec3 getPosition(const float meanAnomaly) const; - glm::vec3 getTangent(const float meanAnomaly) const; - glm::mat4 getLookAlongMatrix(const float meanAnomaly) const; + // + // do we actual want to use mean anomaly for this? we could instead pass in time and + // the gravitational parameter which are easier numbers to come by, i think? at least, + // the mean motion is determined from the mean motion and time time + //glm::vec3 getPosition(double gravitationalParameter, double time) const; + + glm::dvec3 getPosition(const double meanAnomaly) const; + glm::dvec3 getTangent(const double meanAnomaly) const; + //glm::mat4 getLookAlongMatrix(const double meanAnomaly) const; private: Vector6 _keplerianElements; - const float getEccentricAnomaly(const float meanAnomaly) const; - const float getTrueAnomaly(const float meanAnomaly) const; - const Vector6 getCartesianCoordinates(const float meanAnomaly) const; + const double getEccentricAnomaly(const double meanAnomaly) const; + const double getTrueAnomaly(const double meanAnomaly) const; + const Vector6 getCartesianCoordinates(const double meanAnomaly) const; }; diff --git a/lib/skein/include/skein/particlemap.h b/lib/skein/include/skein/particlemap.h index b9b5aa8..9f5e5ff 100644 --- a/lib/skein/include/skein/particlemap.h +++ b/lib/skein/include/skein/particlemap.h @@ -14,7 +14,7 @@ class ParticleMap // providing these as two methods keeps things unambiguous - manipulating particles // is just done for setup const Particle& getParticle(const std::string& id) const; - glm::vec3 getParticlePosition(const std::string& id, double time) const; + glm::dvec3 getParticlePosition(const std::string& id, double time) const; const Orbit& getOrbit(const std::string& id) const; void setParticle(const Particle& particle); diff --git a/lib/skein/src/orbit.cpp b/lib/skein/src/orbit.cpp index c7f4e97..7bdfe61 100644 --- a/lib/skein/src/orbit.cpp +++ b/lib/skein/src/orbit.cpp @@ -10,115 +10,115 @@ Orbit::Orbit() _keplerianElements.resize(6); } -float Orbit::getSemiMajorAxis() const +double Orbit::getSemiMajorAxis() const { return _keplerianElements[astro::semiMajorAxisIndex]; } -void Orbit::setSemiMajorAxis(float semiMajorAxis) +void Orbit::setSemiMajorAxis(double semiMajorAxis) { _keplerianElements[astro::semiMajorAxisIndex] = semiMajorAxis; } -float Orbit::getEccentricity() const +double Orbit::getEccentricity() const { return _keplerianElements[astro::eccentricityIndex]; } -void Orbit::setEccentricity(float eccentricity) +void Orbit::setEccentricity(double eccentricity) { _keplerianElements[astro::eccentricityIndex] = eccentricity; } -float Orbit::getInclination() const +double Orbit::getInclination() const { return _keplerianElements[astro::inclinationIndex]; } -void Orbit::setInclination(float inclination) +void Orbit::setInclination(double inclination) { _keplerianElements[astro::inclinationIndex] = inclination; } -float Orbit::getArgumentOfPeriapsis() const +double Orbit::getArgumentOfPeriapsis() const { return _keplerianElements[astro::argumentOfPeriapsisIndex]; } -void Orbit::setArgumentOfPeriapsis(float argumentOfPeriapsis) +void Orbit::setArgumentOfPeriapsis(double argumentOfPeriapsis) { _keplerianElements[astro::argumentOfPeriapsisIndex] = argumentOfPeriapsis; } -float Orbit::getLongitudeOfAscendingNode() const +double Orbit::getLongitudeOfAscendingNode() const { return _keplerianElements[astro::longitudeOfAscendingNodeIndex]; } -void Orbit::setLongitudeOfAscendingNode(float longitudeOfAscendingNode) +void Orbit::setLongitudeOfAscendingNode(double longitudeOfAscendingNode) { _keplerianElements[astro::longitudeOfAscendingNodeIndex] = longitudeOfAscendingNode; } -glm::mat4 Orbit::getLookAlongMatrix(const float meanAnomaly) const +//glm::mat4 Orbit::getLookAlongMatrix(const double meanAnomaly) const +//{ +// glm::vec3 position = getPosition(meanAnomaly); +// +// // get the tangent of the orbital ellipse +// glm::vec3 tan = getTangent(meanAnomaly); +// // we want to point along the orbit +// glm::vec3 target = position + tan; +// // TODO: this is not 'up' with respect to the orbited body! +// // 'up' is just the normalized position vector because the orbit is centred at the origin +// glm::vec3 up = glm::normalize(position); +// // easy peasy lookAt matrix +// glm::mat4 look = glm::lookAt(position, target, up); +// +// // invert the lookat matrix because it's meant for cameras, cameras work backwards and +// // we are not a camera +// glm::mat4 lookAlong = glm::inverse(look); +// +// return lookAlong; +//} + +const double Orbit::getEccentricAnomaly(const double meanAnomaly) const { - glm::vec3 position = getPosition(meanAnomaly); - - // get the tangent of the orbital ellipse - glm::vec3 tan = getTangent(meanAnomaly); - // we want to point along the orbit - glm::vec3 target = position + tan; - // TODO: this is not 'up' with respect to the orbited body! - // 'up' is just the normalized position vector because the orbit is centred at the origin - glm::vec3 up = glm::normalize(position); - // easy peasy lookAt matrix - glm::mat4 look = glm::lookAt(position, target, up); - - // invert the lookat matrix because it's meant for cameras, cameras work backwards and - // we are not a camera - glm::mat4 lookAlong = glm::inverse(look); - - return lookAlong; -} - -const float Orbit::getEccentricAnomaly(const float meanAnomaly) const -{ - const float eccentricity = _keplerianElements[astro::eccentricityIndex]; - float eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS( + const double eccentricity = _keplerianElements[astro::eccentricityIndex]; + double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS( eccentricity, meanAnomaly, - 100); + 200); return eccentricAnomaly; } // Interpolate a position around the orbit. // t is in range 0..1 and wraps. -glm::vec3 Orbit::getPosition(const float meanAnomaly) const +glm::dvec3 Orbit::getPosition(const double meanAnomaly) const { Vector6 cartesian = getCartesianCoordinates(meanAnomaly); - return glm::vec3( + return glm::dvec3( cartesian[astro::xPositionIndex], cartesian[astro::yPositionIndex], cartesian[astro::zPositionIndex]); } -const float Orbit::getTrueAnomaly(const float meanAnomaly) const +const double Orbit::getTrueAnomaly(const double meanAnomaly) const { - const float eccentricAnomaly = getEccentricAnomaly(meanAnomaly); - const float eccentricity = _keplerianElements[astro::eccentricityIndex]; + const double eccentricAnomaly = getEccentricAnomaly(meanAnomaly); + const double eccentricity = _keplerianElements[astro::eccentricityIndex]; return astro::convertEccentricAnomalyToTrueAnomaly( eccentricAnomaly, eccentricity); } -const Vector6 Orbit::getCartesianCoordinates(const float meanAnomaly) const +const Vector6 Orbit::getCartesianCoordinates(const double meanAnomaly) const { Vector6 kepler(_keplerianElements); kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly); return astro::convertKeplerianToCartesianElements(kepler, 1.0); } -glm::vec3 Orbit::getTangent(const float meanAnomaly) const +glm::dvec3 Orbit::getTangent(const double meanAnomaly) const { - float epsilon = 0.01; - glm::vec3 ahead = getPosition(meanAnomaly + epsilon); - glm::vec3 behind = getPosition(meanAnomaly - epsilon); + double epsilon = 0.01; + glm::dvec3 ahead = getPosition(meanAnomaly + epsilon); + glm::dvec3 behind = getPosition(meanAnomaly - epsilon); return glm::normalize(ahead - behind); } diff --git a/lib/skein/src/particlemap.cpp b/lib/skein/src/particlemap.cpp index 119fddc..842139b 100644 --- a/lib/skein/src/particlemap.cpp +++ b/lib/skein/src/particlemap.cpp @@ -12,17 +12,25 @@ const Orbit& ParticleMap::getOrbit(const std::string& id) const return _orbits.at(id); } -glm::vec3 ParticleMap::getParticlePosition(const std::string& id, double time) const +glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) const { + if (_orbits.find(id) == _orbits.end()) + return {0,0,0}; + + // TODO: actually get stuff based on physics + //const Particle& parent = _relationships.at(id); + + // how do we get the gravitational parameter? + // is it just G * M + // where M is the mass of the parent + //const double u = getGravitationalParameter(parent.getId()); + //const double mass = parent.getMass(); + + // how do we get the particle position from the time + gravitational parameter? + // TODO: actually nest stuff so position is determined from all parents - - if (_orbits.find(id) != _orbits.end()) - { - const Orbit& orbit = _orbits.at(id); - return orbit.getPosition(time); - } - - return {0,0,0}; + const Orbit& orbit = _orbits.at(id); + return orbit.getPosition(time); } void ParticleMap::setParticle(const Particle& particle)