feat: calculate mean orbital motion from M

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:

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;

lib/skein/include/skein/particlemap.h

@@ -22,6 +22,7 @@ class ParticleMap
 
     private:
         std::map<std::string, Particle> _particles;
+        // TODO: improve on this very vague name
         std::map<std::string, std::string> _relationships;
         std::map<std::string, Orbit> _orbits;
 };

lib/skein/src/orbit.cpp

@@ -2,6 +2,7 @@
 
 #include "astro/stateVectorIndices.hpp"
 #include "astro/orbitalElementConversions.hpp"
+#include "astro/twoBodyMethods.hpp"
 
 #include <glm/gtc/matrix_transform.hpp>
 
@@ -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);
+//}
 

lib/skein/src/particle.cpp

@@ -1,5 +1,7 @@
 #include "skein/particle.h"
 
+#include <astro/constants.hpp>
+
 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;
+}

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)

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);
 

src/main.cpp

@@ -103,10 +103,10 @@ int main()
 
     // set parameters of moon's orbit around earth
     Orbit orbit;
-    double semiMajorAxis =  3.84748e9;
+    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;

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;