wip: store true anomaly in kepler elements

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 getPosition(double gravitationalParameter) const;
-    glm::dvec3 getPositionFromMeanAnomaly(double meanAnomaly, 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!                       <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;
 
-    const double getEccentricAnomaly(const double meanAnomaly) const;
+    const double getEccentricAnomaly() const;
-    const double getTrueAnomaly(const double meanAnomaly) 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,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];
+    // TODO: get the eccentric anomaly from the true anomaly?
-    double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS(
+//Real convertTrueAnomalyToEllipticalEccentricAnomaly(const Real trueAnomaly,
-        eccentricity,
+                                                    //const Real eccentricity)
-        meanAnomaly,
+
-        200);
+    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);
+    //Vector6 kepler(_keplerianElements);
-    const double eccentricity = _keplerianElements[astro::eccentricityIndex];
+    //kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly);
-    return astro::convertEccentricAnomalyToTrueAnomaly(
+    return astro::convertKeplerianToCartesianElements(_keplerianElements, gravitationalParameter);
-        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 motion from time
-    kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly);
+    // we assume true anomaly is 0 at t=0 so we this motion is the mean anomaly
-    return astro::convertKeplerianToCartesianElements(kepler, gravitationalParameter);
+    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

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
+//}

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

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;