feat: update orbit based on current velocity

wip: planet position is not conserved

lib/skein/include/skein/orbit.h

@@ -25,14 +25,20 @@ public:
     void setLongitudeOfAscendingNode(double longitudeOfAscendingNode);
 
     glm::dvec3 getPosition(double gravitationalParameter, double time) const;
-    glm::dvec3 getPositionFromMeanAnomaly(double meanAnomaly) 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;
     //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:
     Vector6 _keplerianElements;
 
     const double getEccentricAnomaly(const double meanAnomaly) const;
     const double getTrueAnomaly(const double meanAnomaly) const;
-    const Vector6 getCartesianCoordinates(const double meanAnomaly) const;
 };

lib/skein/include/skein/particlemap.h

@@ -15,11 +15,14 @@ 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);
 
     private:
         std::map<std::string, Particle> _particles;

lib/skein/src/orbit.cpp

@@ -88,9 +88,9 @@ const double Orbit::getEccentricAnomaly(const double meanAnomaly) const
     return eccentricAnomaly;
 }
 
-glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly) const
+glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly, double gravitationalParameter) const
 {
-    Vector6 cartesian = getCartesianCoordinates(meanAnomaly);
+    Vector6 cartesian = getCartesianCoordinates(meanAnomaly, gravitationalParameter);
     return glm::dvec3(
         cartesian[astro::xPositionIndex],
         cartesian[astro::yPositionIndex],
@@ -98,10 +98,28 @@ glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly) const
 }
 
 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);
+    return glm::dvec3(
+        cartesian[astro::xVelocityIndex],
+        cartesian[astro::yVelocityIndex],
+        cartesian[astro::zVelocityIndex]);
+}
+
+glm::dvec3 Orbit::getVelocity(double gravitationalParameter, double time) const
+{
+    return getVelocityFromMeanAnomaly(getMeanAnomaly(gravitationalParameter, time), gravitationalParameter);
+}
+
+const double Orbit::getMeanAnomaly(const double gravitationalParameter, const double time) const
 {
     double meanMotion = astro::computeKeplerMeanMotion(getSemiMajorAxis(), gravitationalParameter);
-    double meanAnomaly = meanMotion * time;
-    return getPositionFromMeanAnomaly(meanAnomaly);
+    return meanMotion * time;
 }
 
 const double Orbit::getTrueAnomaly(const double meanAnomaly) const
@@ -113,11 +131,11 @@ const double Orbit::getTrueAnomaly(const double meanAnomaly) const
         eccentricity);
 }
 
-const Vector6 Orbit::getCartesianCoordinates(const double meanAnomaly) const
+const Vector6 Orbit::getCartesianCoordinates(const double meanAnomaly, const double gravitationalParameter) const
 {
     Vector6 kepler(_keplerianElements);
     kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly);
-    return astro::convertKeplerianToCartesianElements(kepler, 1.0);
+    return astro::convertKeplerianToCartesianElements(kepler, gravitationalParameter);
 }
 
 //glm::dvec3 Orbit::getTangent(const double meanAnomaly) const

lib/skein/src/particlemap.cpp

@@ -1,5 +1,6 @@
 #include "skein/particlemap.h"
 #include "skein/particle.h"
+#include "astro/orbitalElementConversions.hpp"
 
 const Particle& ParticleMap::getParticle(const std::string& id) const
 {
@@ -17,6 +18,25 @@ const Orbit& ParticleMap::getOrbit(const std::string& id) const
     return _orbits.at(id);
 }
 
+glm::dvec3 ParticleMap::getParticleLocalPosition(const std::string& id, double time) const
+{
+    if (_orbits.find(id) == _orbits.end())
+        return {0,0,0};
+
+    const Particle* particle = &(_particles.at(id));
+    glm::dvec3 pos(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();
+    const Orbit& orbit = _orbits.at(id);
+
+    pos += orbit.getPosition(gravitationalParameter, time);
+
+    return pos;
+}
+
 glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) const
 {
     if (_orbits.find(id) == _orbits.end())
@@ -27,6 +47,7 @@ glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time)
 
     do
     {
+        // TODO: this is shadowing a parameter with the same name
         const std::string& id = particle->getId();
         const std::string& parentId = _relationships.at(id);
         const Particle& parent = _particles.at(parentId);
@@ -49,6 +70,54 @@ glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time)
     return pos;
 }
 
+glm::dvec3 ParticleMap::getParticleLocalVelocity(const std::string& id, double time) const
+{
+    if (_orbits.find(id) == _orbits.end())
+        return {0,0,0};
+
+    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();
+    const Orbit& orbit = _orbits.at(id);
+
+    vel += orbit.getVelocity(gravitationalParameter, time);
+
+    return vel;
+}
+
+void ParticleMap::setParticleLocalVelocity(const std::string& id, const double& time, const glm::dvec3& velocity)
+{
+    // we want to convert a local 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(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 Vector6 keplerian = astro::convertCartesianToKeplerianElements(cartesian, gravitationalParameter);
+
+    Orbit newOrbit;
+    newOrbit.setSemiMajorAxis(keplerian[astro::semiMajorAxisIndex]);
+    newOrbit.setEccentricity(keplerian[astro::eccentricityIndex]);
+    newOrbit.setInclination(keplerian[astro::inclinationIndex]);
+    newOrbit.setArgumentOfPeriapsis(keplerian[astro::argumentOfPeriapsisIndex]);
+    newOrbit.setLongitudeOfAscendingNode(keplerian[astro::longitudeOfAscendingNodeIndex]);
+
+    setRelationship(parentId, id, newOrbit);
+}
+
 void ParticleMap::setParticle(const Particle& particle)
 {
     _particles.insert({particle.getId(), particle});
@@ -58,6 +127,6 @@ void ParticleMap::setRelationship(const std::string& parentId, const std::string
 {
     // map children to parent - children can only have one parent, so we use the map to
     // identify the parent directly using the child as the key
-    _relationships.insert({childId, parentId});
-    _orbits.insert({childId, orbit});
+    _relationships.insert_or_assign(childId, parentId);
+    _orbits.insert_or_assign(childId, orbit);
 }

src/main.cpp

@@ -35,20 +35,23 @@
 struct Input
 {
     bool pauseTime;
+    bool printCartesian;
 } input;
 
 
 void keyCallback(GLFWwindow* window, int key, int scancode, int action, int mods)
 {
-    if (key == GLFW_KEY_SPACE && action == GLFW_PRESS)
+    if (action == GLFW_PRESS)
     {
-        input.pauseTime = true;
+        input.pauseTime = key == GLFW_KEY_SPACE;
+        input.printCartesian = key == GLFW_KEY_W;
     }
 }
 
 void clearInput()
 {
     input.pauseTime = false;
+    input.printCartesian = false;
 }
 
 // now should always increase linearly with real world time, this should not be modified by input
@@ -69,6 +72,11 @@ void updateTime()
     engineTime.now = now;
 }
 
+void printVector(const glm::dvec3& vector)
+{
+    std::cout << "(" << vector.x << ", " << vector.y << ", " << vector.z << ")";
+}
+
 int main()
 {
     GLFWwindow* window = nullptr;
@@ -150,7 +158,24 @@ int main()
         }
         else
         {
-            printf("maneouvre mode!\n");
+            if (input.printCartesian)
+            {
+
+                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;
+
+                printVector(v);
+                std::cout << " -> ";
+                printVector(newVelocity);
+                std::cout << std::endl;
+
+                map.setParticleLocalVelocity("station", time, newVelocity);
+            }
         }
 
         // rendering

src/orbitvisualizer.cpp

@@ -36,7 +36,7 @@ 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);
+        glm::vec3 pos = orbit.getPositionFromMeanAnomaly(t, 1);
         pos += basePos;
 
         // Vertices come out of the library with X and Y being in the 'flat' plane. Re-order them