wip: store true anomaly and mean anomaly at epoch

todo: account for discrepancy in argument of periapsis

lib/skein/include/skein/orbit.h

@@ -17,28 +17,40 @@ public:
     double getInclination() const;
     double getArgumentOfPeriapsis() const;
     double getLongitudeOfAscendingNode() const;
+    double getTrueAnomaly() const;
+    double getMeanAnomalyAtEpoch() 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);
+    void setMeanAnomalyAtEpoch(double meanAnomalyAtEpoch);
 
-    glm::dvec3 getPosition(double gravitationalParameter, double time) const;
+    glm::dvec3 getPosition(double gravitationalParameter) 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;
 
-    //glm::dvec3 getTangent(const double meanAnomaly) const;
+    const Vector6 getCartesianCoordinates(const double gravitationalParameter) const;
-    //glm::mat4 getLookAlongMatrix(const double meanAnomaly) const;
-
-    const Vector6 getCartesianCoordinates(const double meanAnomaly, const double gravitationalParameter) const;
     const double getMeanAnomaly(const double gravitationalParameter, const double time) const;
 
-private:
+    void update(double time, double gravitationalParameter);
-    Vector6 _keplerianElements;
 
-    const double getEccentricAnomaly(const double meanAnomaly) const;
+private:
-    const double getTrueAnomaly(const double meanAnomaly) const;
+    /*
+     * N.B.: Order of elements is important!                       <br>
+     * keplerianElements(0) = semiMajorAxis          [m]           <br>
+     * keplerianElements(1) = eccentricity           [-]           <br>
+     * keplerianElements(2) = inclination            [rad]         <br>
+     * keplerianElements(3) = argument of periapsis  [rad]         <br>
+     * keplerianElements(4) = longitude of           [rad]
+     *                         ascending node        [rad]         <br>
+     * keplerianElements(5) = true anomaly           [rad]
+     */
+    Vector6 _keplerianElements;
+    double _meanAnomalyAtEpoch = 0;
+
+    const double getEccentricAnomaly() const;
 };

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<std::string, Particle> _particles;
         // TODO: improve on this very vague name
         std::map<std::string, std::string> _relationships;
         std::map<std::string, Orbit> _orbits;
+
+        //double _time;
 };

lib/skein/src/orbit.cpp

@@ -56,55 +56,45 @@ void Orbit::setLongitudeOfAscendingNode(double longitudeOfAscendingNode)
     _keplerianElements[astro::longitudeOfAscendingNodeIndex] = longitudeOfAscendingNode;
 }
 
-//glm::mat4 Orbit::getLookAlongMatrix(const double meanAnomaly) const
+double Orbit::getTrueAnomaly() 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
 {
-    const double eccentricity = _keplerianElements[astro::eccentricityIndex];
+    return _keplerianElements[astro::trueAnomalyIndex];
-    double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS(
+}
-        eccentricity,
+void Orbit::setTrueAnomaly(double trueAnomaly)
-        meanAnomaly,
+{
-        200);
+    _keplerianElements[astro::trueAnomalyIndex] = trueAnomaly;
+}
+
+double Orbit::getMeanAnomalyAtEpoch() const
+{
+    return _meanAnomalyAtEpoch;
+}
+void Orbit::setMeanAnomalyAtEpoch(double meanAnomalyAtEpoch)
+{
+    _meanAnomalyAtEpoch = meanAnomalyAtEpoch;
+}
+
+const double Orbit::getEccentricAnomaly() const
+{
+    const double trueAnomaly = getTrueAnomaly();
+    const double eccentricity = getEccentricity();
+    double eccentricAnomaly = astro::convertTrueAnomalyToEllipticalEccentricAnomaly(trueAnomaly, eccentricity);
 
     return eccentricAnomaly;
 }
 
-glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly, double gravitationalParameter) const
+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,27 +112,25 @@ const double Orbit::getMeanAnomaly(const double gravitationalParameter, const do
     return meanMotion * time;
 }
 
-const double Orbit::getTrueAnomaly(const double meanAnomaly) const
+const Vector6 Orbit::getCartesianCoordinates(const double gravitationalParameter) const
 {
-    const double eccentricAnomaly = getEccentricAnomaly(meanAnomaly);
+    return astro::convertKeplerianToCartesianElements(_keplerianElements, gravitationalParameter);
-    const double eccentricity = _keplerianElements[astro::eccentricityIndex];
-    return astro::convertEccentricAnomalyToTrueAnomaly(
-        eccentricAnomaly,
-        eccentricity);
 }
 
-const Vector6 Orbit::getCartesianCoordinates(const double meanAnomaly, const double gravitationalParameter) const
+void Orbit::update(double time, double gravitationalParameter)
 {
-    Vector6 kepler(_keplerianElements);
+    // compute mean anomaly from time and mean anomaly at epoch
-    kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly);
+    const double semiMajorAxis = getSemiMajorAxis();
-    return astro::convertKeplerianToCartesianElements(kepler, gravitationalParameter);
+    const double meanMotion = astro::computeKeplerMeanMotion(semiMajorAxis, gravitationalParameter);
+    const double meanAnomaly = meanMotion * time + _meanAnomalyAtEpoch;
+
+    // compute elliptical eccentric anomaly from mean anomaly
+    const double eccentricity = getEccentricity();
+    const double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS(
+        eccentricity,
+        meanAnomaly,
+        200);
+
+    setTrueAnomaly(astro::convertEccentricAnomalyToTrueAnomaly(eccentricAnomaly, eccentricity));
 }
 
-//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/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())
@@ -78,7 +80,6 @@ glm::dvec3 ParticleMap::getParticleLocalVelocity(const std::string& id, double t
     const Particle* particle = &(_particles.at(id));
     glm::dvec3 vel(0,0,0);
 
-    //const std::string& id = particle->getId();
     const std::string& parentId = _relationships.at(id);
     const Particle& parent = _particles.at(parentId);
     const double gravitationalParameter = parent.getGravitationalParameter();
@@ -89,23 +90,24 @@ 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.
+    // convert a position and velocity into a new set of keplerian elements.
 
-    // first, get the local position of the e
     const Orbit& orbit = _orbits.at(id);
     const std::string& parentId = _relationships.at(id);
     const Particle& parent = _particles.at(parentId);
     const double gravitationalParameter = parent.getGravitationalParameter();
-    const double meanAnomaly = orbit.getMeanAnomaly(gravitationalParameter, time);
+    Vector6 cartesian(orbit.getCartesianCoordinates(gravitationalParameter));
-    Vector6 cartesian(orbit.getCartesianCoordinates(meanAnomaly, gravitationalParameter));
 
     // update velocity, leave position as it is
     cartesian[astro::xVelocityIndex] = velocity.x;
     cartesian[astro::yVelocityIndex] = velocity.y;
     cartesian[astro::zVelocityIndex] = velocity.z;
 
+    const double originalMeanAnomaly = orbit.getMeanAnomaly(gravitationalParameter, time);
+    const double originalMeanAnomalyAtEpoch = orbit.getMeanAnomalyAtEpoch();
+
     const Vector6 keplerian = astro::convertCartesianToKeplerianElements(cartesian, gravitationalParameter);
 
     Orbit newOrbit;
@@ -114,6 +116,11 @@ void ParticleMap::setParticleLocalVelocity(const std::string& id, const double&
     newOrbit.setInclination(keplerian[astro::inclinationIndex]);
     newOrbit.setArgumentOfPeriapsis(keplerian[astro::argumentOfPeriapsisIndex]);
     newOrbit.setLongitudeOfAscendingNode(keplerian[astro::longitudeOfAscendingNodeIndex]);
+    newOrbit.setTrueAnomaly(keplerian[astro::trueAnomalyIndex]);
+
+    const double newMeanAnomaly = newOrbit.getMeanAnomaly(gravitationalParameter, time);
+    const double newMeanAnomalyAtEpoch = originalMeanAnomalyAtEpoch - (newMeanAnomaly - originalMeanAnomaly);
+    newOrbit.setMeanAnomalyAtEpoch(newMeanAnomalyAtEpoch);
 
     setRelationship(parentId, id, newOrbit);
 }
@@ -130,3 +137,4 @@ void ParticleMap::setRelationship(const std::string& parentId, const std::string
     _relationships.insert_or_assign(childId, parentId);
     _orbits.insert_or_assign(childId, orbit);
 }
+

src/main.cpp

@@ -35,7 +35,7 @@
 struct Input
 {
     bool pauseTime;
-    bool printCartesian;
+    bool updateOrbit;
 } input;
 
 
@@ -44,14 +44,14 @@ void keyCallback(GLFWwindow* window, int key, int scancode, int action, int mods
     if (action == GLFW_PRESS)
     {
         input.pauseTime = key == GLFW_KEY_SPACE;
-        input.printCartesian = key == GLFW_KEY_W;
+        input.updateOrbit = key == GLFW_KEY_W;
     }
 }
 
 void clearInput()
 {
     input.pauseTime = false;
-    input.printCartesian = false;
+    input.updateOrbit = false;
 }
 
 // now should always increase linearly with real world time, this should not be modified by input
@@ -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
@@ -158,21 +166,13 @@ int main()
         }
         else
         {
-            if (input.printCartesian)
+            if (input.updateOrbit)
             {
 
                 const glm::dvec3 p = map.getParticleLocalPosition("station", time);
-                printVector(p);
-                std::cout << std::endl;
 
-                // 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 << " -> ";
-                printVector(newVelocity);
-                std::cout << std::endl;
 
                 map.setParticleLocalVelocity("station", time, newVelocity);
             }
@@ -182,6 +182,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);

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

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;
+    //const float meanAnomaly = time;
-    glm::vec3 pos = _map.getParticlePosition(_particleId, meanAnomaly);
+    glm::vec3 pos = _map.getParticlePosition(_particleId);
     float y = pos.z;
     pos.z = pos.y;
     pos.y = y;

todo.md

@@ -30,3 +30,10 @@ Currently we get a planet's position by querying the map with a time parameter.
 
 I think to resolve this we will have to update the particle map (particle simulation?) 
 
+We assume that true anomaly is 0 at t=0 to avoid having to consider the previous position. We can't actually make this assumption when we change the shape of the orbit, because the period of the orbit changes as well. As such the calculated position doesn't take into account motion that happened when the orbit was a different shape.
+
+Now we set the true anomaly when setting the particle velocity. However, the particle still jumps around due to us always calculating the particle position based on an absolute time, rather than taking into account the new period.
+
+We need to calculate the particle's position as a function of the absolute time AND the specific particle's true anomaly at t=0. We need to update the latter when we modify the orbit.
+
+We actually want to calculate the particle's mean anomaly at epoch when we update orbital parameters. This is because it is the mean anomaly which is computed as a function of time. When we update the orbital parameters, we can compute the mean anomly from the true anomaly on both occasions, which gives us the different in the mean anomaly at epoch before and after the update.