From 02c768e1b4dd58e6aee028dcf9bd30c89f7bf803 Mon Sep 17 00:00:00 2001 From: Cat Flynn Date: Wed, 25 Sep 2024 20:49:46 +0100 Subject: [PATCH] wip: store true anomaly in kepler elements --- lib/skein/include/skein/orbit.h | 24 +++++-- lib/skein/include/skein/particlemap.h | 12 ++-- lib/skein/src/orbit.cpp | 92 ++++++++++++++++++++------- lib/skein/src/particlemap.cpp | 23 +++++-- src/main.cpp | 26 ++++++-- src/orbitvisualizer.cpp | 8 ++- src/particlevisualizer.cpp | 4 +- 7 files changed, 142 insertions(+), 47 deletions(-) diff --git a/lib/skein/include/skein/orbit.h b/lib/skein/include/skein/orbit.h index 89000df..050977c 100644 --- a/lib/skein/include/skein/orbit.h +++ b/lib/skein/include/skein/orbit.h @@ -17,15 +17,18 @@ public: double getInclination() const; double getArgumentOfPeriapsis() const; double getLongitudeOfAscendingNode() const; + double getTrueAnomaly() const; void setSemiMajorAxis(double semiMajorAxis); void setEccentricity(double eccentricity); void setInclination(double inclination); void setArgumentOfPeriapsis(double argumentOfPeriapsis); void setLongitudeOfAscendingNode(double longitudeOfAscendingNode); + void setTrueAnomaly(double trueAnomaly); - glm::dvec3 getPosition(double gravitationalParameter, double time) const; - glm::dvec3 getPositionFromMeanAnomaly(double meanAnomaly, double gravitationalParameter) const; + glm::dvec3 getPosition(double gravitationalParameter) const; + //glm::dvec3 getPosition(double gravitationalParameter, double time) const; + //glm::dvec3 getPositionFromMeanAnomaly(double meanAnomaly, double gravitationalParameter) const; glm::dvec3 getVelocity(double gravitationalParameter, double time) const; glm::dvec3 getVelocityFromMeanAnomaly(double meanAnomaly, double gravitationalParameter) const; @@ -33,12 +36,23 @@ public: //glm::dvec3 getTangent(const double meanAnomaly) const; //glm::mat4 getLookAlongMatrix(const double meanAnomaly) const; - const Vector6 getCartesianCoordinates(const double meanAnomaly, const double gravitationalParameter) const; + const Vector6 getCartesianCoordinates(const double gravitationalParameter) const; const double getMeanAnomaly(const double gravitationalParameter, const double time) const; + void update(double time, double gravitationalParameter); + private: + /* + * N.B.: Order of elements is important!
+ * keplerianElements(0) = semiMajorAxis [m]
+ * keplerianElements(1) = eccentricity [-]
+ * keplerianElements(2) = inclination [rad]
+ * keplerianElements(3) = argument of periapsis [rad]
+ * keplerianElements(4) = longitude of [rad] + * ascending node [rad]
+ * keplerianElements(5) = true anomaly [rad] + */ Vector6 _keplerianElements; - const double getEccentricAnomaly(const double meanAnomaly) const; - const double getTrueAnomaly(const double meanAnomaly) const; + const double getEccentricAnomaly() const; }; diff --git a/lib/skein/include/skein/particlemap.h b/lib/skein/include/skein/particlemap.h index afbb1b5..bc5c55f 100644 --- a/lib/skein/include/skein/particlemap.h +++ b/lib/skein/include/skein/particlemap.h @@ -15,18 +15,22 @@ class ParticleMap // is just done for setup const Particle& getParticle(const std::string& id) const; const Particle& getParent(const std::string& id) const; - glm::dvec3 getParticleLocalPosition(const std::string& id, double time) const; - glm::dvec3 getParticlePosition(const std::string& id, double time) const; - glm::dvec3 getParticleLocalVelocity(const std::string& id, double time) const; const Orbit& getOrbit(const std::string& id) const; void setParticle(const Particle& particle); void setRelationship(const std::string& parentId, const std::string& childId, const Orbit& orbit); - void setParticleLocalVelocity(const std::string& id, const double& time, const glm::dvec3& velocity); + //void update(double time); + void setParticleLocalVelocity(const std::string& id, double time, const glm::dvec3& velocity); + + glm::dvec3 getParticlePosition(const std::string& id) const; + glm::dvec3 getParticleLocalPosition(const std::string& id, double time) const; + glm::dvec3 getParticleLocalVelocity(const std::string& id, double time) const; private: std::map _particles; // TODO: improve on this very vague name std::map _relationships; std::map _orbits; + + //double _time; }; diff --git a/lib/skein/src/orbit.cpp b/lib/skein/src/orbit.cpp index c251162..de2d632 100644 --- a/lib/skein/src/orbit.cpp +++ b/lib/skein/src/orbit.cpp @@ -56,6 +56,15 @@ void Orbit::setLongitudeOfAscendingNode(double longitudeOfAscendingNode) _keplerianElements[astro::longitudeOfAscendingNodeIndex] = longitudeOfAscendingNode; } +double Orbit::getTrueAnomaly() const +{ + return _keplerianElements[astro::trueAnomalyIndex]; +} +void Orbit::setTrueAnomaly(double trueAnomaly) +{ + _keplerianElements[astro::trueAnomalyIndex] = trueAnomaly; +} + //glm::mat4 Orbit::getLookAlongMatrix(const double meanAnomaly) const //{ // glm::vec3 position = getPosition(meanAnomaly); @@ -77,34 +86,48 @@ void Orbit::setLongitudeOfAscendingNode(double longitudeOfAscendingNode) // return lookAlong; //} -const double Orbit::getEccentricAnomaly(const double meanAnomaly) const +const double Orbit::getEccentricAnomaly() const { - const double eccentricity = _keplerianElements[astro::eccentricityIndex]; - double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS( - eccentricity, - meanAnomaly, - 200); + // TODO: get the eccentric anomaly from the true anomaly? +//Real convertTrueAnomalyToEllipticalEccentricAnomaly(const Real trueAnomaly, + //const Real eccentricity) + + const double trueAnomaly = getTrueAnomaly(); + const double eccentricity = getEccentricity(); + //double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS( + // eccentricity, + // meanAnomaly, + // 200); + double eccentricAnomaly = astro::convertTrueAnomalyToEllipticalEccentricAnomaly(trueAnomaly, eccentricity); return eccentricAnomaly; } -glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly, double gravitationalParameter) const +//glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly, double gravitationalParameter) const +//{ +// Vector6 cartesian = getCartesianCoordinates(meanAnomaly, gravitationalParameter); +// return glm::dvec3( +// cartesian[astro::xPositionIndex], +// cartesian[astro::yPositionIndex], +// cartesian[astro::zPositionIndex]); +//} +// +//glm::dvec3 Orbit::getPosition(double gravitationalParameter, double time) const +//{ +// return getPositionFromMeanAnomaly(getMeanAnomaly(gravitationalParameter, time), gravitationalParameter); +//} +glm::dvec3 Orbit::getPosition(double gravitationalParameter) const { - Vector6 cartesian = getCartesianCoordinates(meanAnomaly, gravitationalParameter); + const Vector6 cartesian = getCartesianCoordinates(gravitationalParameter); return glm::dvec3( cartesian[astro::xPositionIndex], cartesian[astro::yPositionIndex], cartesian[astro::zPositionIndex]); } -glm::dvec3 Orbit::getPosition(double gravitationalParameter, double time) const -{ - return getPositionFromMeanAnomaly(getMeanAnomaly(gravitationalParameter, time), gravitationalParameter); -} - glm::dvec3 Orbit::getVelocityFromMeanAnomaly(const double meanAnomaly, double gravitationalParameter) const { - Vector6 cartesian = getCartesianCoordinates(meanAnomaly, gravitationalParameter); + Vector6 cartesian = getCartesianCoordinates(gravitationalParameter); return glm::dvec3( cartesian[astro::xVelocityIndex], cartesian[astro::yVelocityIndex], @@ -122,20 +145,41 @@ const double Orbit::getMeanAnomaly(const double gravitationalParameter, const do return meanMotion * time; } -const double Orbit::getTrueAnomaly(const double meanAnomaly) const +//const double Orbit::getTrueAnomaly(const double meanAnomaly) const +//{ +// const double eccentricAnomaly = getEccentricAnomaly(meanAnomaly); +// const double eccentricity = _keplerianElements[astro::eccentricityIndex]; +// return astro::convertEccentricAnomalyToTrueAnomaly( +// eccentricAnomaly, +// eccentricity); +//} + +const Vector6 Orbit::getCartesianCoordinates(const double gravitationalParameter) const { - const double eccentricAnomaly = getEccentricAnomaly(meanAnomaly); - const double eccentricity = _keplerianElements[astro::eccentricityIndex]; - return astro::convertEccentricAnomalyToTrueAnomaly( - eccentricAnomaly, - eccentricity); + //Vector6 kepler(_keplerianElements); + //kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly); + return astro::convertKeplerianToCartesianElements(_keplerianElements, gravitationalParameter); } -const Vector6 Orbit::getCartesianCoordinates(const double meanAnomaly, const double gravitationalParameter) const +void Orbit::update(double time, double gravitationalParameter) { - Vector6 kepler(_keplerianElements); - kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly); - return astro::convertKeplerianToCartesianElements(kepler, gravitationalParameter); + // compute mean motion from time + // we assume true anomaly is 0 at t=0 so we this motion is the mean anomaly + const double semiMajorAxis = getSemiMajorAxis(); + const double meanMotion = astro::computeKeplerMeanMotion(semiMajorAxis, gravitationalParameter); + const double meanAnomaly = meanMotion * time; + + // compute elliptical eccentric anomaly from mean anomaly + const double eccentricity = getEccentricity(); + const double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS( + eccentricity, + meanAnomaly, + 200); + + // compute true anomaly from eccentric anomaly + const double trueAnomaly = astro::convertEccentricAnomalyToTrueAnomaly(eccentricAnomaly, eccentricity); + + setTrueAnomaly(trueAnomaly); } //glm::dvec3 Orbit::getTangent(const double meanAnomaly) const diff --git a/lib/skein/src/particlemap.cpp b/lib/skein/src/particlemap.cpp index 20ed187..f9086af 100644 --- a/lib/skein/src/particlemap.cpp +++ b/lib/skein/src/particlemap.cpp @@ -32,12 +32,13 @@ glm::dvec3 ParticleMap::getParticleLocalPosition(const std::string& id, double t const double gravitationalParameter = parent.getGravitationalParameter(); const Orbit& orbit = _orbits.at(id); - pos += orbit.getPosition(gravitationalParameter, time); + //pos += orbit.getPosition(gravitationalParameter, time); + pos += orbit.getPosition(gravitationalParameter); return pos; } -glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) const +glm::dvec3 ParticleMap::getParticlePosition(const std::string& id) const { if (_orbits.find(id) == _orbits.end()) return {0,0,0}; @@ -54,7 +55,8 @@ glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) const double gravitationalParameter = parent.getGravitationalParameter(); const Orbit& orbit = _orbits.at(id); - pos += orbit.getPosition(gravitationalParameter, time); + //pos += orbit.getPosition(gravitationalParameter, time); + pos += orbit.getPosition(gravitationalParameter); auto it = _relationships.find(parentId); if (it != _relationships.end()) @@ -89,7 +91,7 @@ glm::dvec3 ParticleMap::getParticleLocalVelocity(const std::string& id, double t return vel; } -void ParticleMap::setParticleLocalVelocity(const std::string& id, const double& time, const glm::dvec3& velocity) +void ParticleMap::setParticleLocalVelocity(const std::string& id, double time, const glm::dvec3& velocity) { // we want to convert a local position and velocity into a new set of keplerian elements. @@ -99,7 +101,8 @@ void ParticleMap::setParticleLocalVelocity(const std::string& id, const double& const Particle& parent = _particles.at(parentId); const double gravitationalParameter = parent.getGravitationalParameter(); const double meanAnomaly = orbit.getMeanAnomaly(gravitationalParameter, time); - Vector6 cartesian(orbit.getCartesianCoordinates(meanAnomaly, gravitationalParameter)); + //Vector6 cartesian(orbit.getCartesianCoordinates(meanAnomaly, gravitationalParameter)); + Vector6 cartesian(orbit.getCartesianCoordinates(gravitationalParameter)); // update velocity, leave position as it is cartesian[astro::xVelocityIndex] = velocity.x; @@ -130,3 +133,13 @@ void ParticleMap::setRelationship(const std::string& parentId, const std::string _relationships.insert_or_assign(childId, parentId); _orbits.insert_or_assign(childId, orbit); } + +//void ParticleMap::update(double time) +//{ +// //const double interval = time - _time; +// //_time = time; +// +// // what should happen when we update the time? +// +// // assume that all particles' true anomaly is 0 at t=0 +//} diff --git a/src/main.cpp b/src/main.cpp index 990c2a9..0048622 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -94,6 +94,7 @@ int main() moonOrbit.setInclination(glm::radians(5.15)); // radians moonOrbit.setArgumentOfPeriapsis(318.15); // in the case of the moon these last two values are moonOrbit.setLongitudeOfAscendingNode(60.0); // pretty much constantly changing so use whatever + moonOrbit.setTrueAnomaly(0.0); // set parameters of satellite orbit around moon Orbit stationOrbit; @@ -102,6 +103,7 @@ int main() stationOrbit.setInclination(glm::radians(89.0)); stationOrbit.setArgumentOfPeriapsis(43.2); stationOrbit.setLongitudeOfAscendingNode(239.7); + stationOrbit.setTrueAnomaly(0.0); ParticleMap map; map.setParticle({"earth", 5.9e24}); @@ -140,11 +142,17 @@ int main() // apply input if (input.pauseTime) { - mode++; - if (mode > MODE_MANEOUVRE) + if (mode == MODE_MANEOUVRE) { - mode = 0; + mode = MODE_ORBIT; + std::cout << "play"; } + else if (mode == MODE_ORBIT) + { + mode = MODE_MANEOUVRE; + std::cout << "pause"; + } + std::cout << std::endl; } // the moon takes like 27 days to orbit earth. to see this in motion, then, we need to @@ -167,7 +175,7 @@ int main() // TODO: what unit is this in lol const glm::dvec3 v = map.getParticleLocalVelocity("station", time); - const glm::dvec3 newVelocity = v * 0.9; + const glm::dvec3 newVelocity = v * 0.99; printVector(v); std::cout << " -> "; @@ -182,6 +190,16 @@ int main() glClearColor(0.2, 0.3, 0.3, 1.0); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + // update orbits with time + const Particle& earth = map.getParticle("earth"); + const Particle& moon = map.getParticle("moon"); + moonOrbit = map.getOrbit("moon"); + stationOrbit = map.getOrbit("station"); + moonOrbit.update(time, earth.getGravitationalParameter()); + stationOrbit.update(time, moon.getGravitationalParameter()); + map.setRelationship("earth", "moon", moonOrbit); + map.setRelationship("moon", "station", stationOrbit); + earthVis.render(time); moonVis.render(time); stationVis.render(time); diff --git a/src/orbitvisualizer.cpp b/src/orbitvisualizer.cpp index 79bf4ea..a45c4a1 100644 --- a/src/orbitvisualizer.cpp +++ b/src/orbitvisualizer.cpp @@ -12,7 +12,7 @@ void OrbitVisualizer::render(const float time) { // Orbit needs to be drawn relative to its parent const Particle& parent = _map.getParent(_particleId); - const glm::vec3 parentPos = _map.getParticlePosition(parent.getId(), time); + const glm::vec3 parentPos = _map.getParticlePosition(parent.getId()); regenerateVertices(parentPos); @@ -27,7 +27,7 @@ void OrbitVisualizer::render(const float time) void OrbitVisualizer::regenerateVertices(const glm::vec3& basePos) { - const Orbit& orbit = _map.getOrbit(_particleId); + Orbit orbit(_map.getOrbit(_particleId)); _vertices.clear(); for (int i = 0; i < _vertexCount; i++) @@ -36,7 +36,9 @@ void OrbitVisualizer::regenerateVertices(const glm::vec3& basePos) // better to actually create a first-class ellipse object and use that to generate // a nice continuous mesh, instead of using orbital positions. float t = (float)i / (float)_vertexCount * 2.0 * _pi; - glm::vec3 pos = orbit.getPositionFromMeanAnomaly(t, 1); + orbit.setTrueAnomaly(t); + //glm::vec3 pos = orbit.getPositionFromMeanAnomaly(t, 1); + glm::vec3 pos = orbit.getPosition(1); pos += basePos; // Vertices come out of the library with X and Y being in the 'flat' plane. Re-order them diff --git a/src/particlevisualizer.cpp b/src/particlevisualizer.cpp index 82fa1bf..3f1fd22 100644 --- a/src/particlevisualizer.cpp +++ b/src/particlevisualizer.cpp @@ -10,8 +10,8 @@ ParticleVisualizer::ParticleVisualizer(const ParticleMap& map, const std::string void ParticleVisualizer::render(float time) { // TODO: get mean anomly from particle which has the mass!! - const float meanAnomaly = time; - glm::vec3 pos = _map.getParticlePosition(_particleId, meanAnomaly); + //const float meanAnomaly = time; + glm::vec3 pos = _map.getParticlePosition(_particleId); float y = pos.z; pos.z = pos.y; pos.y = y;