frag_lit.glsl

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

lib/skein/include/skein/orbit.h

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

lib/skein/src/orbit.cpp

@@ -11,12 +11,6 @@ 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];
@@ -41,7 +35,6 @@ double Orbit::getInclination() const
 }
 void Orbit::setInclination(double inclination)
 {
-    makeSafe(inclination);
     _keplerianElements[astro::inclinationIndex] = inclination;
 }
 
@@ -51,7 +44,6 @@ double Orbit::getArgumentOfPeriapsis() const
 }
 void Orbit::setArgumentOfPeriapsis(double argumentOfPeriapsis)
 {
-    makeSafe(argumentOfPeriapsis);
     _keplerianElements[astro::argumentOfPeriapsisIndex] = argumentOfPeriapsis;
 }
 
@@ -61,7 +53,6 @@ double Orbit::getLongitudeOfAscendingNode() const
 }
 void Orbit::setLongitudeOfAscendingNode(double longitudeOfAscendingNode)
 {
-    makeSafe(longitudeOfAscendingNode);
     _keplerianElements[astro::longitudeOfAscendingNodeIndex] = longitudeOfAscendingNode;
 }
 

lib/skein/src/particlemap.cpp

@@ -119,14 +119,7 @@ void ParticleMap::setParticleLocalVelocity(const std::string& id, double time, c
     newOrbit.setTrueAnomaly(keplerian[astro::trueAnomalyIndex]);
 
     const double newMeanAnomaly = newOrbit.getMeanAnomaly(gravitationalParameter, time);
-    const double meanAnomalyShift = newMeanAnomaly - originalMeanAnomaly;
+    const double newMeanAnomalyAtEpoch = originalMeanAnomalyAtEpoch - (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,8 +1,11 @@
 #include "icosphere.hpp"
-#include "gfx.hpp"
 
+#include <iostream>
 #include <array>
-#include <glm/gtc/matrix_transform.hpp>
+
+#include "glm/gtc/matrix_transform.hpp"
+
+#include "gfx.hpp"
 
 Icosphere::Icosphere(float radius, int subdivisions, GLuint shaderProgram) :
     _shaderProgram(shaderProgram),
@@ -45,6 +48,10 @@ 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,8 +35,7 @@
 struct Input
 {
     bool pauseTime;
-    bool prograde;
+    bool updateOrbit;
-    bool retrograde;
 } input;
 
 
@@ -45,16 +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.prograde = key == GLFW_KEY_W;
+        input.updateOrbit = key == GLFW_KEY_W;
-        input.retrograde = key == GLFW_KEY_S;
     }
 }
 
 void clearInput()
 {
     input.pauseTime = false;
-    input.prograde = false;
+    input.updateOrbit = false;
-    input.retrograde = false;
 }
 
 // now should always increase linearly with real world time, this should not be modified by input
@@ -91,43 +88,36 @@ 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);
+    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({"earth", 1});
+    map.setParticle({"station", 1e6});
-    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);
@@ -176,18 +166,16 @@ int main()
         }
         else
         {
-            //glm::dvec3 v = map.getParticleLocalVelocity("station", time);
+            if (input.updateOrbit)
-            glm::dvec3 v = map.getParticleLocalVelocity("moon", time);
-            if (input.prograde)
             {
-                v *= 1.01;
+
+                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);
             }
-            else if (input.retrograde)
-            {
-                v *= 0.99;
-            }
-            //map.setParticleLocalVelocity("station", time, v);
-            map.setParticleLocalVelocity("moon", time, v);
         }
 
         // rendering
@@ -198,17 +186,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,6 +1,5 @@
 #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)
@@ -24,32 +23,12 @@ 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++)
     {
@@ -57,8 +36,18 @@ 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 = getPositionOnOrbit(t);
+        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
+        // 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,6 +25,5 @@ class OrbitVisualizer
         GLuint _vao;
         std::vector<float> _vertices;
 
-        glm::dvec3 getPositionOnOrbit(double trueAnomaly) const;
         void regenerateVertices(const glm::vec3& basePos);
 };

src/particlevisualizer.cpp

@@ -9,7 +9,14 @@ 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,9 +42,5 @@ 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 swap sides. This would cause the particle to jump to the other side of the planet
-
-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,13 +9,11 @@ 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);