From a570986485ccfe963b25dea33cc9cc3b795e75bc Mon Sep 17 00:00:00 2001 From: Cat Flynn Date: Wed, 21 Aug 2024 15:45:33 +0100 Subject: [PATCH] feat: calculate mean orbital motion from M --- lib/skein/include/skein/orbit.h | 13 +++---------- lib/skein/include/skein/particle.h | 1 + lib/skein/include/skein/particlemap.h | 1 + lib/skein/src/orbit.cpp | 26 ++++++++++++++++---------- lib/skein/src/particle.cpp | 7 +++++++ lib/skein/src/particlemap.cpp | 16 ++++------------ src/gfx.cpp | 4 +++- src/main.cpp | 11 +++++++---- src/orbitvisualizer.cpp | 6 +++++- 9 files changed, 47 insertions(+), 38 deletions(-) diff --git a/lib/skein/include/skein/orbit.h b/lib/skein/include/skein/orbit.h index 4b679af..f2e80fd 100644 --- a/lib/skein/include/skein/orbit.h +++ b/lib/skein/include/skein/orbit.h @@ -24,16 +24,9 @@ public: void setArgumentOfPeriapsis(double argumentOfPeriapsis); void setLongitudeOfAscendingNode(double longitudeOfAscendingNode); - // TODO: meanAnomaly in all these arguments actually means ellipticalMeanAnomaly, - // will have to change that when adding non-ellipctical orbits - don't get confused! - // - // 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::dvec3 getPosition(double gravitationalParameter, double time) const; + glm::dvec3 getPositionFromMeanAnomaly(double meanAnomaly) const; + //glm::dvec3 getTangent(const double meanAnomaly) const; //glm::mat4 getLookAlongMatrix(const double meanAnomaly) const; private: diff --git a/lib/skein/include/skein/particle.h b/lib/skein/include/skein/particle.h index ec00aa5..05d44e7 100644 --- a/lib/skein/include/skein/particle.h +++ b/lib/skein/include/skein/particle.h @@ -13,6 +13,7 @@ public: const std::string& getId() const; double getMass() const; + double getGravitationalParameter() const; private: const std::string _id; diff --git a/lib/skein/include/skein/particlemap.h b/lib/skein/include/skein/particlemap.h index 9f5e5ff..17bb2ff 100644 --- a/lib/skein/include/skein/particlemap.h +++ b/lib/skein/include/skein/particlemap.h @@ -22,6 +22,7 @@ class ParticleMap private: std::map _particles; + // TODO: improve on this very vague name std::map _relationships; std::map _orbits; }; diff --git a/lib/skein/src/orbit.cpp b/lib/skein/src/orbit.cpp index 7bdfe61..c1ab218 100644 --- a/lib/skein/src/orbit.cpp +++ b/lib/skein/src/orbit.cpp @@ -2,6 +2,7 @@ #include "astro/stateVectorIndices.hpp" #include "astro/orbitalElementConversions.hpp" +#include "astro/twoBodyMethods.hpp" #include @@ -87,9 +88,7 @@ const double Orbit::getEccentricAnomaly(const double meanAnomaly) const return eccentricAnomaly; } -// Interpolate a position around the orbit. -// t is in range 0..1 and wraps. -glm::dvec3 Orbit::getPosition(const double meanAnomaly) const +glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly) const { Vector6 cartesian = getCartesianCoordinates(meanAnomaly); return glm::dvec3( @@ -98,6 +97,13 @@ glm::dvec3 Orbit::getPosition(const double meanAnomaly) const cartesian[astro::zPositionIndex]); } +glm::dvec3 Orbit::getPosition(double gravitationalParameter, double time) const +{ + double meanMotion = astro::computeKeplerMeanMotion(getSemiMajorAxis(), gravitationalParameter); + double meanAnomaly = meanMotion * time; + return getPositionFromMeanAnomaly(meanAnomaly); +} + const double Orbit::getTrueAnomaly(const double meanAnomaly) const { const double eccentricAnomaly = getEccentricAnomaly(meanAnomaly); @@ -114,11 +120,11 @@ const Vector6 Orbit::getCartesianCoordinates(const double meanAnomaly) const return astro::convertKeplerianToCartesianElements(kepler, 1.0); } -glm::dvec3 Orbit::getTangent(const double meanAnomaly) const -{ - double epsilon = 0.01; - glm::dvec3 ahead = getPosition(meanAnomaly + epsilon); - glm::dvec3 behind = getPosition(meanAnomaly - epsilon); - return glm::normalize(ahead - behind); -} +//glm::dvec3 Orbit::getTangent(const double meanAnomaly) const +//{ +// 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/particle.cpp b/lib/skein/src/particle.cpp index bbdf5cc..5cff006 100644 --- a/lib/skein/src/particle.cpp +++ b/lib/skein/src/particle.cpp @@ -1,5 +1,7 @@ #include "skein/particle.h" +#include + Particle::Particle(const std::string& id, double mass) : _id(id), _mass(mass) { @@ -19,3 +21,8 @@ double Particle::getMass() const { return _mass; } + +double Particle::getGravitationalParameter() const +{ + return _mass * astro::ASTRO_GRAVITATIONAL_CONSTANT; +} diff --git a/lib/skein/src/particlemap.cpp b/lib/skein/src/particlemap.cpp index 842139b..55e12de 100644 --- a/lib/skein/src/particlemap.cpp +++ b/lib/skein/src/particlemap.cpp @@ -1,7 +1,6 @@ #include "skein/particlemap.h" #include "skein/particle.h" - const Particle& ParticleMap::getParticle(const std::string& id) const { return _particles.at(id); @@ -17,20 +16,13 @@ glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) 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? + const std::string& parentId = _relationships.at(id); + const Particle& parent = _particles.at(parentId); + const double gravitationalParameter = parent.getGravitationalParameter(); // TODO: actually nest stuff so position is determined from all parents const Orbit& orbit = _orbits.at(id); - return orbit.getPosition(time); + return orbit.getPosition(gravitationalParameter, time); } void ParticleMap::setParticle(const Particle& particle) diff --git a/src/gfx.cpp b/src/gfx.cpp index 4461a1a..d0bf87c 100644 --- a/src/gfx.cpp +++ b/src/gfx.cpp @@ -128,7 +128,9 @@ void updateViewMatrix(GLuint shaderProgram, float time) glm::mat4 view = glm::mat4(1.0); // Rotation - constexpr float angle = glm::radians(10.0); + // TODO: disconnect application and simulation time + constexpr float angle = 0; + //constexpr float angle = glm::radians(10.0); glm::vec3 axis = glm::vec3(0.0, 1.0, 0.0); view = glm::rotate(view, angle * time, axis); diff --git a/src/main.cpp b/src/main.cpp index 4df46cf..89937cc 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -103,10 +103,10 @@ int main() // set parameters of moon's orbit around earth Orbit orbit; - double semiMajorAxis = 3.84748e9; - orbit.setSemiMajorAxis(semiMajorAxis); // metres + double semiMajorAxis = 3.84748e8; + orbit.setSemiMajorAxis(semiMajorAxis); // metres orbit.setEccentricity(0.055); - orbit.setInclination(glm::radians(5.15)); // radians? + orbit.setInclination(glm::radians(5.15)); // radians orbit.setArgumentOfPeriapsis(318.15); // in the case of the moon these last two values are orbit.setLongitudeOfAscendingNode(60.0); // pretty much constantly changing so use whatever @@ -149,8 +149,11 @@ int main() } } + // the moon takes like 27 days to orbit earth. to see this in motion, then, we need to + // increase the speed of time by 60 * 60 * 24 to see 1 day per second. + const double speed = 60 * 60 * 24; + // only update time if playing the orbiting animation - const double speed = 0.5; if (animation == ANIM_ORBITING) { time = getTime() * speed; diff --git a/src/orbitvisualizer.cpp b/src/orbitvisualizer.cpp index 258f970..5f52a96 100644 --- a/src/orbitvisualizer.cpp +++ b/src/orbitvisualizer.cpp @@ -27,8 +27,12 @@ void OrbitVisualizer::regenerateVertices() for (int i = 0; i < _vertexCount; i++) { + // TODO: this method of getting ellipse vertices is a huge hack. it would be + // 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.getPosition(t); + glm::vec3 pos = _orbit.getPositionFromMeanAnomaly(t); + // Vertices come out of the library with X and Y being in the 'flat' plane. Re-order them // here such that Z is up. float y = pos.z;