doc: update todo

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

@@ -53,4 +53,5 @@ private:
     double _meanAnomalyAtEpoch = 0;
 
     const double getEccentricAnomaly() const;
+    static void makeSafe(double& value);
 };

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,6 +61,7 @@ double Orbit::getLongitudeOfAscendingNode() const
 }
 void Orbit::setLongitudeOfAscendingNode(double longitudeOfAscendingNode)
 {
+    makeSafe(longitudeOfAscendingNode);
     _keplerianElements[astro::longitudeOfAscendingNodeIndex] = longitudeOfAscendingNode;
 }
 

lib/skein/src/particlemap.cpp

@@ -119,7 +119,14 @@ void ParticleMap::setParticleLocalVelocity(const std::string& id, double time, c
     newOrbit.setTrueAnomaly(keplerian[astro::trueAnomalyIndex]);
 
     const double newMeanAnomaly = newOrbit.getMeanAnomaly(gravitationalParameter, time);
-    const double newMeanAnomalyAtEpoch = originalMeanAnomalyAtEpoch - (newMeanAnomaly - originalMeanAnomaly);
+    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);

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,7 +35,8 @@
 struct Input
 {
     bool pauseTime;
-    bool updateOrbit;
+    bool prograde;
+    bool retrograde;
 } input;
 
 
@@ -44,14 +45,16 @@ void keyCallback(GLFWwindow* window, int key, int scancode, int action, int mods
     if (action == GLFW_PRESS)
     {
         input.pauseTime = key == GLFW_KEY_SPACE;
-        input.updateOrbit = key == GLFW_KEY_W;
+        input.prograde = key == GLFW_KEY_W;
+        input.retrograde = key == GLFW_KEY_S;
     }
 }
 
 void clearInput()
 {
     input.pauseTime = false;
-    input.updateOrbit = false;
+    input.prograde = false;
+    input.retrograde = false;
 }
 
 // now should always increase linearly with real world time, this should not be modified by input
@@ -88,36 +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.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.setEccentricity(0.055);
+    //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.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;
-    stationOrbit.setSemiMajorAxis(5e7);
-    stationOrbit.setEccentricity(0.6);
-    stationOrbit.setInclination(glm::radians(89.0));
-    stationOrbit.setArgumentOfPeriapsis(43.2);
-    stationOrbit.setLongitudeOfAscendingNode(239.7);
-    stationOrbit.setTrueAnomaly(0.0);
+    //Orbit stationOrbit;
+    //stationOrbit.setSemiMajorAxis(5e7);
+    //stationOrbit.setEccentricity(0.6);
+    //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});
-    map.setParticle({"moon", 7.3e22});
-    map.setParticle({"station", 1e6});
+    //map.setParticle({"earth", 5.9e24});
+    //map.setParticle({"moon", 7.3e22});
+    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);
@@ -166,16 +176,18 @@ int main()
         }
         else
         {
-            if (input.updateOrbit)
+            //glm::dvec3 v = map.getParticleLocalVelocity("station", time);
+            glm::dvec3 v = map.getParticleLocalVelocity("moon", time);
+            if (input.prograde)
             {
-
-                const glm::dvec3 p = map.getParticleLocalPosition("station", time);
-
-                const glm::dvec3 v = map.getParticleLocalVelocity("station", time);
-                const glm::dvec3 newVelocity = v * 0.99;
-
-                map.setParticleLocalVelocity("station", time, newVelocity);
+                v *= 1.01;
             }
+            else if (input.retrograde)
+            {
+                v *= 0.99;
+            }
+            //map.setParticleLocalVelocity("station", time, v);
+            map.setParticleLocalVelocity("moon", time, v);
         }
 
         // rendering
@@ -186,17 +198,17 @@ int main()
         const Particle& earth = map.getParticle("earth");
         const Particle& moon = map.getParticle("moon");
         moonOrbit = map.getOrbit("moon");
-        stationOrbit = map.getOrbit("station");
+        //stationOrbit = map.getOrbit("station");
         moonOrbit.update(time, earth.getGravitationalParameter());
-        stationOrbit.update(time, moon.getGravitationalParameter());
+        //stationOrbit.update(time, moon.getGravitationalParameter());
         map.setRelationship("earth", "moon", moonOrbit);
-        map.setRelationship("moon", "station", stationOrbit);
+        //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)
@@ -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)
 {
-    Orbit orbit(_map.getOrbit(_particleId));
-
     _vertices.clear();
     for (int i = 0; i < _vertexCount; i++)
     {
@@ -36,18 +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;
-        orbit.setTrueAnomaly(t);
-        //glm::vec3 pos = orbit.getPositionFromMeanAnomaly(t, 1);
-        glm::vec3 pos = orbit.getPosition(1);
+        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!!
-    //const float meanAnomaly = time;
     glm::vec3 pos = _map.getParticlePosition(_particleId);
-    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

@@ -42,5 +42,9 @@ 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
 
-Test: Update the velocity such that the apo- and periapses of the orbit swap sides. This would cause the particle to jump to the other side of the planet
+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);