doc: update todo

todo: account for discrepancy in argument of periapsis

frag_lit.glsl

@@ -2,6 +2,7 @@
 
 in vec3 Normal;
 in vec3 FragPos;
+in mat4x4 ModelViewProjection;
 
 out vec4 FragColor;
 
@@ -17,7 +18,7 @@ void main()
 
     // Directional lighting
     vec3 directionalLightColor = vec3(1.0, 1.0, 1.0);
-    vec3 lightPos = vec3(10.0, 7.0, -8.0);
+    vec3 lightPos = (ModelViewProjection * vec4(6.0, -7.0, 8.0, 1.0)).xyz;
     vec3 normal = normalize(Normal);
     vec3 lightDir = normalize(lightPos - FragPos);
     float diff = max(dot(normal, lightDir), 0.0);

lib/skein/include/skein/orbit.h

@@ -17,22 +17,41 @@ 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) const;
+
-    //glm::dvec3 getTangent(const double meanAnomaly) const;
+    glm::dvec3 getVelocity(double gravitationalParameter, double time) const;
-    //glm::mat4 getLookAlongMatrix(const double meanAnomaly) const;
+    glm::dvec3 getVelocityFromMeanAnomaly(double meanAnomaly, 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;
+    double _meanAnomalyAtEpoch = 0;
 
-    const double getEccentricAnomaly(const double meanAnomaly) const;
+    const double getEccentricAnomaly() const;
-    const double getTrueAnomaly(const double meanAnomaly) const;
+    static void makeSafe(double& value);
-    const Vector6 getCartesianCoordinates(const double meanAnomaly) const;
 };

lib/skein/include/skein/particlemap.h

@@ -15,15 +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 getParticlePosition(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 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

@@ -11,6 +11,12 @@ Orbit::Orbit()
     _keplerianElements.resize(6);
 }
 
+// some values cause singularities when they are zero. so, make them not zero
+void Orbit::makeSafe(double& value)
+{
+    value = value == 0.0 ? 0.0000001 : value;
+}
+
 double Orbit::getSemiMajorAxis() const
 {
     return _keplerianElements[astro::semiMajorAxisIndex];
@@ -35,6 +41,7 @@ double Orbit::getInclination() const
 }
 void Orbit::setInclination(double inclination)
 {
+    makeSafe(inclination);
     _keplerianElements[astro::inclinationIndex] = inclination;
 }
 
@@ -44,6 +51,7 @@ double Orbit::getArgumentOfPeriapsis() const
 }
 void Orbit::setArgumentOfPeriapsis(double argumentOfPeriapsis)
 {
+    makeSafe(argumentOfPeriapsis);
     _keplerianElements[astro::argumentOfPeriapsisIndex] = argumentOfPeriapsis;
 }
 
@@ -53,78 +61,85 @@ double Orbit::getLongitudeOfAscendingNode() const
 }
 void Orbit::setLongitudeOfAscendingNode(double longitudeOfAscendingNode)
 {
+    makeSafe(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) const
+glm::dvec3 Orbit::getPosition(double gravitationalParameter) const
 {
-    Vector6 cartesian = getCartesianCoordinates(meanAnomaly);
+    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
+glm::dvec3 Orbit::getVelocityFromMeanAnomaly(const double meanAnomaly, double gravitationalParameter) const
+{
+    Vector6 cartesian = getCartesianCoordinates(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 meanMotion * time;
-    return getPositionFromMeanAnomaly(meanAnomaly);
 }
 
-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
+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, 1.0);
+    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

@@ -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,7 +18,27 @@ const Orbit& ParticleMap::getOrbit(const std::string& id) const
     return _orbits.at(id);
 }
 
-glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) const
+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);
+    pos += orbit.getPosition(gravitationalParameter);
+
+    return pos;
+}
+
+glm::dvec3 ParticleMap::getParticlePosition(const std::string& id) const
 {
     if (_orbits.find(id) == _orbits.end())
         return {0,0,0};
@@ -27,13 +48,15 @@ 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);
         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())
@@ -49,6 +72,66 @@ 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& 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, double time, const glm::dvec3& velocity)
+{
+    // convert a position and velocity into a new set of keplerian elements.
+
+    const Orbit& orbit = _orbits.at(id);
+    const std::string& parentId = _relationships.at(id);
+    const Particle& parent = _particles.at(parentId);
+    const double gravitationalParameter = parent.getGravitationalParameter();
+    Vector6 cartesian(orbit.getCartesianCoordinates(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;
+    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]);
+    newOrbit.setTrueAnomaly(keplerian[astro::trueAnomalyIndex]);
+
+    const double newMeanAnomaly = newOrbit.getMeanAnomaly(gravitationalParameter, time);
+    const double meanAnomalyShift = newMeanAnomaly - originalMeanAnomaly;
+    const double periapsisShift = newOrbit.getArgumentOfPeriapsis() - orbit.getArgumentOfPeriapsis();
+    //const double ascendingNodeShift = newOrbit.getLongitudeOfAscendingNode() - orbit.getLongitudeOfAscendingNode();
+
+    const double newMeanAnomalyAtEpoch = originalMeanAnomalyAtEpoch
+        - meanAnomalyShift
+        - periapsisShift;
+        //- ascendingNodeShift;
+    newOrbit.setMeanAnomalyAtEpoch(newMeanAnomalyAtEpoch);
+
+    setRelationship(parentId, id, newOrbit);
+}
+
 void ParticleMap::setParticle(const Particle& particle)
 {
     _particles.insert({particle.getId(), particle});
@@ -58,6 +141,7 @@ 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});
+    _relationships.insert_or_assign(childId, parentId);
-    _orbits.insert({childId, orbit});
+    _orbits.insert_or_assign(childId, orbit);
 }
+

src/icosphere.cpp

@@ -1,12 +1,9 @@
 #include "icosphere.hpp"
-
-#include <iostream>
-#include <array>
-
-#include "glm/gtc/matrix_transform.hpp"
-
 #include "gfx.hpp"
 
+#include <array>
+#include <glm/gtc/matrix_transform.hpp>
+
 Icosphere::Icosphere(float radius, int subdivisions, GLuint shaderProgram) :
     _shaderProgram(shaderProgram),
     _position({})
@@ -48,10 +45,6 @@ void Icosphere::generateVertices(float radius, int subdivisions)
     {
         triangles = subdivide(vertices, triangles);
     }
-    std::cout <<
-        "subdivisions: " << subdivisions <<
-        " vertices: " << vertices.size() <<
-        " triangles: " << triangles.size() << std::endl;
 
     // Scale vertices by radius after subdivision as subdivision happens on a
     // unit sphere

src/main.cpp

@@ -35,20 +35,26 @@
 struct Input
 {
     bool pauseTime;
+    bool prograde;
+    bool retrograde;
 } 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.prograde = key == GLFW_KEY_W;
+        input.retrograde = key == GLFW_KEY_S;
     }
 }
 
 void clearInput()
 {
     input.pauseTime = false;
+    input.prograde = false;
+    input.retrograde = false;
 }
 
 // now should always increase linearly with real world time, this should not be modified by input
@@ -69,6 +75,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;
@@ -80,34 +91,43 @@ int main()
 
     // set parameters of moon's orbit around earth
     Orbit moonOrbit;
-    double moonOrbitSemiMajorAxis =  3.84748e8;
+    //double moonOrbitSemiMajorAxis =  3.84748e8;
+    double moonOrbitSemiMajorAxis =  1;
     moonOrbit.setSemiMajorAxis(moonOrbitSemiMajorAxis); // metres
-    moonOrbit.setEccentricity(0.055);
+    //moonOrbit.setEccentricity(0.055);
-    moonOrbit.setInclination(glm::radians(5.15));       // radians
+    //moonOrbit.setInclination(glm::radians(5.15));       // radians
-    moonOrbit.setArgumentOfPeriapsis(318.15);           // in the case of the moon these last two values are
+    //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.setLongitudeOfAscendingNode(60.0);        // pretty much constantly changing so use whatever
+    moonOrbit.setEccentricity(0.01);
+    moonOrbit.setInclination(0.0);
+    moonOrbit.setArgumentOfPeriapsis(0.0);
+    moonOrbit.setLongitudeOfAscendingNode(0.0);
+    moonOrbit.setTrueAnomaly(0.0);
 
     // set parameters of satellite orbit around moon
-    Orbit stationOrbit;
+    //Orbit stationOrbit;
-    stationOrbit.setSemiMajorAxis(5e7);
+    //stationOrbit.setSemiMajorAxis(5e7);
-    stationOrbit.setEccentricity(0.6);
+    //stationOrbit.setEccentricity(0.6);
-    stationOrbit.setInclination(glm::radians(89.0));
+    //stationOrbit.setInclination(glm::radians(89.0));
-    stationOrbit.setArgumentOfPeriapsis(43.2);
+    //stationOrbit.setArgumentOfPeriapsis(43.2);
-    stationOrbit.setLongitudeOfAscendingNode(239.7);
+    //stationOrbit.setLongitudeOfAscendingNode(239.7);
+    //stationOrbit.setTrueAnomaly(0.0);
 
     ParticleMap map;
-    map.setParticle({"earth", 5.9e24});
+    //map.setParticle({"earth", 5.9e24});
-    map.setParticle({"moon", 7.3e22});
+    //map.setParticle({"moon", 7.3e22});
-    map.setParticle({"station", 1e6});
+    map.setParticle({"earth", 1});
+    map.setParticle({"moon", .5});
+    //map.setParticle({"station", 1e6});
     map.setRelationship("earth", "moon", moonOrbit);
-    map.setRelationship("moon", "station", stationOrbit);
+    //map.setRelationship("moon", "station", stationOrbit);
 
     float scale = moonOrbitSemiMajorAxis * 1.1;
     ParticleVisualizer earthVis(map, "earth", 0.1, litProgram, unlitProgram, scale);
     ParticleVisualizer moonVis(map, "moon", 0.02, litProgram, unlitProgram, scale);
-    ParticleVisualizer stationVis(map, "station", 0.01, litProgram, unlitProgram, scale);
+    //ParticleVisualizer stationVis(map, "station", 0.01, litProgram, unlitProgram, scale);
     OrbitVisualizer moonOrbitVis(map, "moon", unlitProgram, scale);
-    OrbitVisualizer stationOrbitVis(map, "station", unlitProgram, scale);
+    //OrbitVisualizer stationOrbitVis(map, "station", unlitProgram, scale);
 
     // register input
     glfwSetKeyCallback(window, keyCallback);
@@ -132,11 +152,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
@@ -150,18 +176,39 @@ int main()
         }
         else
         {
-            printf("maneouvre mode!\n");
+            //glm::dvec3 v = map.getParticleLocalVelocity("station", time);
+            glm::dvec3 v = map.getParticleLocalVelocity("moon", time);
+            if (input.prograde)
+            {
+                v *= 1.01;
+            }
+            else if (input.retrograde)
+            {
+                v *= 0.99;
+            }
+            //map.setParticleLocalVelocity("station", time, v);
+            map.setParticleLocalVelocity("moon", time, v);
         }
 
         // rendering
         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);
+        //stationVis.render(time);
         moonOrbitVis.render(time);
-        stationOrbitVis.render(time);
+        //stationOrbitVis.render(time);
 
         glfwSwapBuffers(window);
     }

src/orbitvisualizer.cpp

@@ -1,5 +1,6 @@
 #include "orbitvisualizer.hpp"
 #include "gfx.hpp"
+#include "icosphere.hpp"
 
 OrbitVisualizer::OrbitVisualizer(const ParticleMap& map, const std::string& particleId, const GLuint shaderProgram, float scale)
     : _map(map), _particleId(particleId), _shaderProgram(shaderProgram), _scale(scale)
@@ -12,7 +13,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);
 
@@ -23,12 +24,32 @@ void OrbitVisualizer::render(const float time)
     glBindBuffer(GL_ARRAY_BUFFER, _vbo);
 
     glDrawArrays(GL_LINE_LOOP, 0, _vertices.size() / 3);
+
+    glm::dvec3 position;
+
+    Icosphere apoapsis(0.01, 2, _shaderProgram);
+    position = getPositionOnOrbit(3.142);
+    position /= _scale;
+    apoapsis.setPosition(position);
+    apoapsis.render(time);
+
+    Icosphere periapsis(0.01, 2, _shaderProgram);
+    position = getPositionOnOrbit(0.0);
+    position /= _scale;
+    periapsis.setPosition(position);
+    periapsis.render(time);
+
+}
+
+glm::dvec3 OrbitVisualizer::getPositionOnOrbit(double trueAnomaly) const
+{
+    Orbit orbit(_map.getOrbit(_particleId));
+    orbit.setTrueAnomaly(trueAnomaly);
+    return orbit.getPosition(1);
 }
 
 void OrbitVisualizer::regenerateVertices(const glm::vec3& basePos)
 {
-    const Orbit& orbit = _map.getOrbit(_particleId);
-
     _vertices.clear();
     for (int i = 0; i < _vertexCount; i++)
     {
@@ -36,16 +57,8 @@ 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 = getPositionOnOrbit(t);
         pos += basePos;
-
-        // 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;
-        pos.z = pos.y;
-        pos.y = y;
-        pos.z *= -1;
-
         pos /= _scale;
 
         _vertices.push_back(pos.x);

src/orbitvisualizer.hpp

@@ -25,5 +25,6 @@ class OrbitVisualizer
         GLuint _vao;
         std::vector<float> _vertices;
 
+        glm::dvec3 getPositionOnOrbit(double trueAnomaly) const;
         void regenerateVertices(const glm::vec3& basePos);
 };

src/particlevisualizer.cpp

@@ -9,14 +9,7 @@ ParticleVisualizer::ParticleVisualizer(const ParticleMap& map, const std::string
 
 void ParticleVisualizer::render(float time)
 {
-    // TODO: get mean anomly from particle which has the mass!!
+    glm::vec3 pos = _map.getParticlePosition(_particleId);
-    const float meanAnomaly = time;
-    glm::vec3 pos = _map.getParticlePosition(_particleId, meanAnomaly);
-    float y = pos.z;
-    pos.z = pos.y;
-    pos.y = y;
-    pos.z *= -1;
-
     pos /= _scale;
 
     // TODO: extract widget to its own visualizer since we know it wants an orbit but we

todo.md

@@ -0,0 +1,50 @@
+Can we determine orbits entirely from an orbital state vector?
+Can we use a state vector to generate initial orbital elements?
+Can we use a state vector to modify an existing orbit?
+
+The [perifocal reference frame](https://orbital-mechanics.space/classical-orbital-elements/perifocal-frame.html) has a right-handed coordinate system where p=x,q=y and w=z.
+
+We want to be able to represent three orbit types individually as elliptical orbits are fundamentally cyclical in a way that para/hyperbolic orbits are not.
+Hyperbolic orbits more useful than parabolic since they allow us to actually leave a planet.
+They need to be positioned at a unique point in time, where elliptical orbits do not.
+What additional orbital element information is required to position them in time?
+
+Orbits take place in a (non-inertial?) reference frame. How should we model these?
+We need a set of 3D axes and a clock for a reference frame.
+What elements of our existing implementation would change?
+It should be possible to get reference frames from a particle.
+Should we be constructing orbits with a reference frame?
+
+- [x] Make C properly pause time
+- [x] Make W print out the normalized Cartesian velocity
+- [x] Verify that the velocity is in opposite directions on opposite sides of the orbit
+- [x] Make W convert the orbital elements to cartesian 
+    - [x] add 10% to the velocity
+    - [x] then convert cartesian back to kepler elements and update the orbit
+
+Now we can update the shape of the orbit but the position of the planet is not preserved. This is probably because we are calculated mean anomaly as a function of time, instead of taking into account the planet's existing elements.
+
+Therefore when converting from cartesian to keplerian elements, we lose the information about the true anomaly. If we conserve this information when setting the planet's position in following frames we should maintain the planet's position within its orbit.
+
+Currently we get a planet's position by querying the map with a time parameter. Instead, we should update the map with a time parameter and then query it. This means that we can update the entire map to the same time value, and then operate on its objects in the context of a single instant. We still avoid directly storing a time value in the map, but we will have complete elements with which to orbital operations.
+
+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.
+
+Argument of Periapsis
+
+Hypothesis: The shifting in position is caused by a the change in the argument of periapsis which is not accounted for when updating the orbit
+
+What happens when the initial orbit is circular?
+
+Test: Update the velocity such that the apo- and periapses of the orbit should swap sides. This should cause the particle to jump to the other side of the planet.
+Test: Position the particle exactly at the apo- or periapsis before updating the velocity. The particle should remain stationary.
+
+

vert.glsl

@@ -9,11 +9,13 @@ uniform mat4x4 _Projection;
 
 out vec3 Normal;
 out vec3 FragPos;
+out mat4x4 ModelViewProjection;
 
 void main()
 {
     vec4 pos = vec4(aPos, 1.0);
     mat4x4 mvp = _Projection * _View * _Model;
+    ModelViewProjection = mvp;
 
     gl_Position = mvp * pos;
     FragPos = vec3(_Model * pos);