]> git.localhorst.tv Git - blobs.git/commitdiff
more orbits and stuff
authorDaniel Karbach <daniel.karbach@localhorst.tv>
Sat, 11 Nov 2017 13:03:23 +0000 (14:03 +0100)
committerDaniel Karbach <daniel.karbach@localhorst.tv>
Sat, 11 Nov 2017 13:03:59 +0000 (14:03 +0100)
15 files changed:
src/app/Assets.hpp
src/app/MasterState.hpp
src/app/states.cpp
src/blobs.cpp
src/graphics/SunSurface.hpp [new file with mode: 0644]
src/graphics/shader.cpp
src/world/Body.hpp
src/world/Orbit.hpp [new file with mode: 0644]
src/world/Simulation.hpp
src/world/sim.cpp
src/world/world.cpp
tst/assert.hpp [new file with mode: 0644]
tst/test.cpp [new file with mode: 0644]
tst/world/OrbitTest.cpp [new file with mode: 0644]
tst/world/OrbitTest.hpp [new file with mode: 0644]

index ed3c2e2605d8dd6bda3c16f0d4962d34fcfb1fd7..c80726a07c05adbd50058e6a41d6eea65e18cca9 100644 (file)
@@ -3,6 +3,7 @@
 
 #include "../graphics/ArrayTexture.hpp"
 #include "../graphics/PlanetSurface.hpp"
+#include "../graphics/SunSurface.hpp"
 
 
 namespace blobs {
@@ -16,6 +17,7 @@ struct Assets {
 
        struct {
                graphics::PlanetSurface planet_surface;
+               graphics::SunSurface sun_surface;
        } shaders;
 
        Assets();
index fb245cb62b2afb9e080cdd272e75e992a632508c..f62c77309ea7ea4a04133c5b1013e5accebded10 100644 (file)
@@ -33,6 +33,8 @@ public:
 private:
        void OnResize(int w, int h) override;
 
+       void OnKeyDown(const SDL_KeyboardEvent &) override;
+
        void OnUpdate(int dt) override;
        void OnRender(graphics::Viewport &) override;
 
@@ -48,6 +50,7 @@ private:
 
        int remain;
        int thirds;
+       bool paused;
 
 };
 
index ef8b5140141c5cbe8481712b86a78da95f37292b..67c073a0d5292b5c6548b0c04044d810f6e5d817 100644 (file)
@@ -16,8 +16,14 @@ MasterState::MasterState(Assets &assets, world::Simulation &sim) noexcept
 , reference(&sim.Root())
 , cam()
 , remain(0)
-, thirds(0) {
-       cam.View(glm::lookAt(glm::vec3(2.0f, 3.0f, 10.0f), glm::vec3(0.0f, 0.0f, 0.0f), glm::vec3(0.0f, 1.0f, 0.0f)));
+, thirds(0)
+, paused(false) {
+       // sunset view: standing in the center of surface 0 (+Z), looking west (-X)
+       //cam.View(glm::lookAt(glm::vec3(0.0f, 0.0f, 5.6f), glm::vec3(-1.0f, 0.0f, 5.6f), glm::vec3(0.0f, 0.0f, 1.0f)));
+       // sunrise view: standing in the center of surface 0 (+Z), looking east (+X)
+       cam.View(glm::lookAt(glm::vec3(0.0f, 0.0f, 5.6f), glm::vec3(1.0f, 0.0f, 5.6f), glm::vec3(0.0f, 0.0f, 1.0f)));
+       // far out, looking at planet
+       //cam.View(glm::lookAt(glm::vec3(10.0f, 10.0f, 50.0f), glm::vec3(0.0f, 0.0f, 0.0f), glm::vec3(0.0f, 1.0f, 0.0f)));
 }
 
 MasterState::~MasterState() noexcept {
@@ -36,7 +42,9 @@ void MasterState::OnUpdate(int dt) {
 }
 
 void MasterState::Tick() {
-       sim.Tick();
+       if (!paused) {
+               sim.Tick();
+       }
        remain -= FrameMS();
        thirds = (thirds + 1) % 3;
 }
@@ -46,13 +54,29 @@ int MasterState::FrameMS() const noexcept {
 }
 
 
+void MasterState::OnKeyDown(const SDL_KeyboardEvent &e) {
+       if (e.keysym.sym == SDLK_p) {
+               paused = !paused;
+       }
+}
+
 void MasterState::OnRender(graphics::Viewport &viewport) {
-       glm::mat4 ppos = reference->ToParent();
+       glm::dmat4 ppos = reference->InverseTransform() * reference->ToParent();
        assets.shaders.planet_surface.Activate();
-       assets.shaders.planet_surface.SetMVP(glm::mat4(1.0f), cam.View(), cam.Projection());
        assets.shaders.planet_surface.SetTexture(assets.textures.tiles);
-       assets.shaders.planet_surface.SetLight(glm::vec3(cam.View() * ppos[3]), glm::vec3(1.0f, 1.0f, 1.0f), 100.0f);
+
+       assets.shaders.planet_surface.SetMVP(glm::mat4(1.0f), cam.View(), cam.Projection());
+       assets.shaders.planet_surface.SetLight(glm::vec3(cam.View() * ppos[3]), glm::vec3(1.0f, 1.0f, 1.0f), 2.0e4f);
        reference->Draw(assets, viewport);
+
+       world::Body *child = reference->Children()[0];
+       assets.shaders.planet_surface.SetMVP(reference->InverseTransform() * child->FromParent() * child->LocalTransform(), cam.View(), cam.Projection());
+       child->Draw(assets, viewport);
+
+       assets.shaders.sun_surface.Activate();
+       assets.shaders.sun_surface.SetMVP(ppos * reference->Parent().LocalTransform(), cam.View(), cam.Projection());
+       assets.shaders.sun_surface.SetLight(glm::vec3(1.0f, 1.0f, 1.0f), 2.0e4f);
+       assets.shaders.sun_surface.Draw();
 }
 
 }
index 5f7bb7674dfb67caa51565040baad6c7c90d43d8..c302159c0124bab310d8a55820192b0024a2ac8a 100644 (file)
@@ -8,6 +8,7 @@
 #include "world/Sun.hpp"
 
 #include <cstdint>
+#include <iostream>
 
 
 using namespace blobs;
@@ -17,14 +18,34 @@ int main(int argc, char *argv[]) {
        app::Assets assets;
 
        world::Sun sun;
-       sun.Mass(1.0e13);
-       world::Simulation sim(sun);
-       world::Planet planet(3);
-       planet.Mass(1.0e7);
-       planet.Inclination(PI * 0.25);
+       sun.Mass(1.0e12);
+       sun.Radius(1.0);
+       sun.SurfaceTilt(glm::dvec2(PI * 0.25, PI * 0.25));
+       sun.AngularMomentum(1.0e9);
+
+       world::Planet planet(11);
        world::GenerateTest(planet);
        planet.SetParent(sun);
-       planet.SemiMajorAxis(10.0);
+       planet.Mass(1.0e9);
+       planet.GetOrbit().SemiMajorAxis(100.0);
+       planet.SurfaceTilt(glm::dvec2(PI * 0.25, PI * 0.25));
+       planet.AxialTilt(glm::dvec2(PI * 0.127, 0.0));
+       planet.AngularMomentum(3.0e9);
+
+       world::Planet moon(3);
+       world::GenerateTest(moon);
+       moon.SetParent(planet);
+       moon.Mass(1.0e6);
+       moon.GetOrbit().SemiMajorAxis(25.0);
+       moon.AngularMomentum(1.0e5);
+
+       world::Simulation sim(sun);
+       sim.AddBody(planet);
+       sim.AddBody(moon);
+
+       std::cout << "length of year: " << planet.OrbitalPeriod() << "s" << std::endl;
+       std::cout << "length of moon cycle: " << moon.OrbitalPeriod() << "s" << std::endl;
+       std::cout << "length of day: " << planet.RotationalPeriod() << "s" << std::endl;
 
        app::MasterState state(assets, sim);
        state.SetReference(planet);
diff --git a/src/graphics/SunSurface.hpp b/src/graphics/SunSurface.hpp
new file mode 100644 (file)
index 0000000..5ce91b1
--- /dev/null
@@ -0,0 +1,68 @@
+#ifndef BLOBS_GRAPHICS_SUNSURFACE_HPP_
+#define BLOBS_GRAPHICS_SUNSURFACE_HPP_
+
+#include "Program.hpp"
+#include "SimpleVAO.hpp"
+
+#include "glm.hpp"
+
+#include <cstdint>
+
+
+namespace blobs {
+namespace graphics {
+
+class ArrayTexture;
+
+class SunSurface {
+
+public:
+       SunSurface();
+       ~SunSurface();
+
+       SunSurface(const SunSurface &) = delete;
+       SunSurface &operator =(const SunSurface &) = delete;
+
+       SunSurface(SunSurface &&) = delete;
+       SunSurface &operator =(SunSurface &&) = delete;
+
+public:
+       void Activate() noexcept;
+
+       void SetMVP(const glm::mat4 &m, const glm::mat4 &v, const glm::mat4 &p) noexcept;
+       void SetLight(const glm::vec3 &color, float strength) noexcept;
+
+       const glm::mat4 &M() const noexcept { return m; }
+       const glm::mat4 &V() const noexcept { return v; }
+       const glm::mat4 &P() const noexcept { return p; }
+       const glm::mat4 &MV() const noexcept { return mv; }
+       const glm::mat4 &MVP() const noexcept { return mvp; }
+
+       void Draw() const noexcept;
+
+private:
+       struct Attributes {
+               glm::vec3 position;
+       };
+       SimpleVAO<Attributes, std::uint8_t> vao;
+       Program prog;
+
+       glm::mat4 m;
+       glm::mat4 v;
+       glm::mat4 p;
+       glm::mat4 mv;
+       glm::mat4 mvp;
+
+       GLuint m_handle;
+       GLuint mv_handle;
+       GLuint mvp_handle;
+
+       GLuint light_color_handle;
+       GLuint light_strength_handle;
+
+};
+
+}
+}
+
+#endif
index b184bb98d08d6009bde2af53c30027b623f67ecc..d771d9f943385f4709108b99c2f58540cd868db2 100644 (file)
@@ -1,6 +1,7 @@
 #include "PlanetSurface.hpp"
 #include "Program.hpp"
 #include "Shader.hpp"
+#include "SunSurface.hpp"
 
 #include "ArrayTexture.hpp"
 #include "../app/init.hpp"
@@ -277,5 +278,153 @@ void PlanetSurface::SetLight(const glm::vec3 &pos, const glm::vec3 &color, float
        prog.Uniform(light_strength_handle, strength);
 }
 
+
+SunSurface::SunSurface()
+: prog() {
+       prog.LoadShader(
+               GL_VERTEX_SHADER,
+               "#version 330 core\n"
+
+               "layout(location = 0) in vec3 vtx_position;\n"
+
+               "uniform mat4 M;\n"
+               "uniform mat4 MV;\n"
+               "uniform mat4 MVP;\n"
+
+               "out vec3 vtx_viewspace;\n"
+
+               "void main() {\n"
+                       "gl_Position = MVP * vec4(vtx_position, 1);\n"
+                       "vtx_viewspace = (MV * vec4(vtx_position, 1)).xyz;\n"
+               "}\n"
+       );
+       prog.LoadShader(
+               GL_FRAGMENT_SHADER,
+               "#version 330 core\n"
+
+               "in vec3 vtx_viewspace;\n"
+
+               "uniform vec3 light_color;\n"
+               "uniform float light_strength;\n"
+
+               "out vec3 color;\n"
+
+               "void main() {\n"
+                       "vec3 to_light = -vtx_viewspace;\n"
+                       "float distance = length(to_light);\n"
+                       //"vec3 light_dir = normalize(to_light);\n"
+                       "float attenuation = light_strength / (distance * distance);\n"
+                       "color = attenuation * light_color;\n"
+               "}\n"
+       );
+       prog.Link();
+       if (!prog.Linked()) {
+               prog.Log(std::cerr);
+               throw std::runtime_error("link program");
+       }
+       m_handle = prog.UniformLocation("M");
+       mv_handle = prog.UniformLocation("MV");
+       mvp_handle = prog.UniformLocation("MVP");
+       light_color_handle = prog.UniformLocation("light_color");
+       light_strength_handle = prog.UniformLocation("light_strength");
+
+       vao.Bind();
+       vao.BindAttributes();
+       vao.EnableAttribute(0);
+       vao.AttributePointer<glm::vec3>(0, false, offsetof(Attributes, position));
+       vao.ReserveAttributes(8, GL_STATIC_DRAW);
+       {
+               auto attrib = vao.MapAttributes(GL_WRITE_ONLY);
+               attrib[0].position = glm::vec3(-1.0f, -1.0f, -1.0f);
+               attrib[1].position = glm::vec3(-1.0f, -1.0f,  1.0f);
+               attrib[2].position = glm::vec3(-1.0f,  1.0f, -1.0f);
+               attrib[3].position = glm::vec3(-1.0f,  1.0f,  1.0f);
+               attrib[4].position = glm::vec3( 1.0f, -1.0f, -1.0f);
+               attrib[5].position = glm::vec3( 1.0f, -1.0f,  1.0f);
+               attrib[6].position = glm::vec3( 1.0f,  1.0f, -1.0f);
+               attrib[7].position = glm::vec3( 1.0f,  1.0f,  1.0f);
+       }
+       vao.BindElements();
+       vao.ReserveElements(36, GL_STATIC_DRAW);
+       {
+               auto element = vao.MapElements(GL_WRITE_ONLY);
+               // -X
+               element[ 0] = 0;
+               element[ 1] = 1;
+               element[ 2] = 2;
+               element[ 3] = 2;
+               element[ 4] = 1;
+               element[ 5] = 3;
+               // -Y
+               element[ 6] = 0;
+               element[ 7] = 4;
+               element[ 8] = 1;
+               element[ 9] = 1;
+               element[10] = 4;
+               element[11] = 5;
+               // -Z
+               element[12] = 0;
+               element[13] = 2;
+               element[14] = 4;
+               element[15] = 4;
+               element[16] = 2;
+               element[17] = 6;
+               // +Z
+               element[18] = 1;
+               element[19] = 5;
+               element[20] = 3;
+               element[21] = 3;
+               element[22] = 5;
+               element[23] = 7;
+               // +Y
+               element[24] = 3;
+               element[25] = 7;
+               element[26] = 2;
+               element[27] = 2;
+               element[28] = 7;
+               element[29] = 6;
+               // +X
+               element[30] = 5;
+               element[31] = 4;
+               element[32] = 7;
+               element[33] = 7;
+               element[34] = 4;
+               element[35] = 6;
+       }
+       vao.Unbind();
+}
+
+SunSurface::~SunSurface() {
+}
+
+void SunSurface::Activate() noexcept {
+       prog.Use();
+       glEnable(GL_DEPTH_TEST);
+       glDepthFunc(GL_LESS);
+       glEnable(GL_CULL_FACE);
+       glDisable(GL_BLEND);
+}
+
+void SunSurface::SetMVP(const glm::mat4 &mm, const glm::mat4 &vv, const glm::mat4 &pp) noexcept {
+       m = mm;
+       v = vv;
+       p = pp;
+       mv = v * m;
+       mvp = p * mv;
+       prog.Uniform(m_handle, m);
+       prog.Uniform(mv_handle, mv);
+       prog.Uniform(mvp_handle, mvp);
+}
+
+void SunSurface::SetLight(const glm::vec3 &color, float strength) noexcept {
+       prog.Uniform(light_color_handle, color);
+       prog.Uniform(light_strength_handle, strength);
+}
+
+void SunSurface::Draw() const noexcept {
+       vao.Bind();
+       vao.DrawTriangles(36);
+}
+
 }
 }
index 2f87d22723cd2239575355bf5153c006ef8d13af..1b9e2e3a2a65036b75f370aa47a58b3730374545 100644 (file)
@@ -1,6 +1,7 @@
 #ifndef BLOBS_WORLD_BODY_HPP_
 #define BLOBS_WORLD_BODY_HPP_
 
+#include "Orbit.hpp"
 #include "../graphics/glm.hpp"
 
 #include <vector>
@@ -40,35 +41,40 @@ public:
        void SetParent(Body &);
        void UnsetParent();
 
-       double Mass() const noexcept;
-       void Mass(double) noexcept;
+       const std::vector<Body *> &Children() const noexcept { return children; }
 
-       double Radius() const noexcept;
-       void Radius(double) noexcept;
+       double Mass() const noexcept { return mass; }
+       void Mass(double m) noexcept { mass = m; }
 
-       double SemiMajorAxis() const noexcept;
-       void SemiMajorAxis(double) noexcept;
+       double Radius() const noexcept { return radius; }
+       void Radius(double r) noexcept { radius = r; }
 
-       double Eccentricity() const noexcept;
-       void Eccentricity(double) noexcept;
+       Orbit &GetOrbit() noexcept { return orbit; }
+       const Orbit &GetOrbit() const noexcept { return orbit; }
 
-       double Inclination() const noexcept;
-       void Inclination(double) noexcept;
+       const glm::dvec2 &SurfaceTilt() const noexcept { return surface_tilt; }
+       void SurfaceTilt(const glm::dvec2 &t) noexcept { surface_tilt = t; }
 
-       double LongitudeAscending() const noexcept;
-       void LongitudeAscending(double) noexcept;
+       const glm::dvec2 &AxialTilt() const noexcept { return axis_tilt; }
+       void AxialTilt(const glm::dvec2 &t) noexcept { axis_tilt = t; }
 
-       double ArgumentPeriapsis() const noexcept;
-       void ArgumentPeriapsis(double) noexcept;
+       double Rotation() const noexcept { return rotation; }
+       void Rotation(double r) noexcept { rotation = r; }
 
-       double MeanAnomaly() const noexcept;
-       void MeanAnomaly(double) noexcept;
+       double AngularMomentum() const noexcept { return angular; }
+       void AngularMomentum(double m) noexcept { angular = m; }
+
+       double Inertia() const noexcept;
 
        double GravitationalParameter() const noexcept;
        double OrbitalPeriod() const noexcept;
+       double RotationalPeriod() const noexcept;
+
+       glm::dmat4 LocalTransform() const noexcept;
+       glm::dmat4 InverseTransform() const noexcept;
 
-       glm::mat4 ToParent() const noexcept;
-       glm::mat4 FromParent() const noexcept;
+       glm::dmat4 ToParent() const noexcept;
+       glm::dmat4 FromParent() const noexcept;
 
        virtual void Draw(app::Assets &, graphics::Viewport &) { }
 
@@ -82,14 +88,11 @@ private:
        std::vector<Body *> children;
        double mass;
        double radius;
-
-       // Orbit
-       double sma; // semi-major axis
-       double ecc; // eccentricity
-       double inc; // inclination
-       double asc; // longitude of ascending node
-       double arg; // argument of periapsis
-       double mna; // mean anomaly (at t=0)
+       Orbit orbit;
+       glm::dvec2 surface_tilt;
+       glm::dvec2 axis_tilt;
+       double rotation;
+       double angular;
 
 };
 
diff --git a/src/world/Orbit.hpp b/src/world/Orbit.hpp
new file mode 100644 (file)
index 0000000..7def521
--- /dev/null
@@ -0,0 +1,53 @@
+#ifndef BLOBS_WORLD_ORBIT_HPP_
+#define BLOBS_WORLD_ORBIT_HPP_
+
+#include "../graphics/glm.hpp"
+
+
+namespace blobs {
+namespace world {
+
+class Orbit {
+
+public:
+       Orbit();
+       ~Orbit();
+
+public:
+       double SemiMajorAxis() const noexcept;
+       Orbit &SemiMajorAxis(double) noexcept;
+
+       double Eccentricity() const noexcept;
+       Orbit &Eccentricity(double) noexcept;
+
+       double Inclination() const noexcept;
+       Orbit &Inclination(double) noexcept;
+
+       double LongitudeAscending() const noexcept;
+       Orbit &LongitudeAscending(double) noexcept;
+
+       double ArgumentPeriapsis() const noexcept;
+       Orbit &ArgumentPeriapsis(double) noexcept;
+
+       double MeanAnomaly() const noexcept;
+       Orbit &MeanAnomaly(double) noexcept;
+
+       /// calculate transformation matrix at position t in the
+       /// orbit, measured in degrees from mean anomaly at t=0
+       glm::dmat4 Matrix(double t) const noexcept;
+       glm::dmat4 InverseMatrix(double t) const noexcept;
+
+private:
+       double sma; // semi-major axis
+       double ecc; // eccentricity
+       double inc; // inclination
+       double asc; // longitude of ascending node
+       double arg; // argument of periapsis
+       double mna; // mean anomaly (at t=0)
+
+};
+
+}
+}
+
+#endif
index 6d22ed3ef9e5087f886e2f1c62e94ef34af57aea..3025a745b5bebb3665174e3599faf99a4c8cccdf 100644 (file)
@@ -1,6 +1,8 @@
 #ifndef BLOBS_WORLD_SIMULATION_HPP_
 #define BLOBS_WORLD_SIMULATION_HPP_
 
+#include <vector>
+
 
 namespace blobs {
 namespace world {
@@ -22,6 +24,8 @@ public:
 public:
        void Tick();
 
+       void AddBody(Body &);
+
        Body &Root() { return root; }
        const Body &Root() const { return root; }
 
@@ -29,6 +33,7 @@ public:
 
 private:
        Body &root;
+       std::vector<Body *> all_bodies;
        double time;
 
 };
index abf2829498f20356115a128801279bf5c53eb10e..46454fd2952197311b72db9bfeace430a44bebdc 100644 (file)
@@ -7,16 +7,27 @@ namespace blobs {
 namespace world {
 
 Simulation::Simulation(Body &r)
-: root(r) {
-       r.SetSimulation(*this);
+: root(r)
+, all_bodies()
+, time(0.0) {
+       AddBody(r);
 }
 
 Simulation::~Simulation() {
 }
 
 
+void Simulation::AddBody(Body &b) {
+       b.SetSimulation(*this);
+       all_bodies.push_back(&b);
+}
+
 void Simulation::Tick() {
-       time += 0.01666666666666666666666666666666;
+       constexpr double dt = 0.01666666666666666666666666666666;
+       time += dt;
+       for (auto body : all_bodies) {
+               body->Rotation(body->Rotation() + dt * body->AngularMomentum() / body->Inertia());
+       }
 }
 
 }
index 5c1135049612f63e1492fa225f867939cc6de0f2..0b1b22df6fca974c5657d9beeb745b270a4e1c23 100644 (file)
@@ -1,4 +1,5 @@
 #include "Body.hpp"
+#include "Orbit.hpp"
 #include "Planet.hpp"
 #include "Simulation.hpp"
 #include "Sun.hpp"
@@ -10,6 +11,8 @@
 
 #include <algorithm>
 #include <cmath>
+#include <glm/gtc/matrix_transform.hpp>
+#include <glm/gtx/euler_angles.hpp>
 #include <glm/gtx/transform.hpp>
 
 using blobs::G;
@@ -30,12 +33,11 @@ Body::Body()
 , children()
 , mass(1.0)
 , radius(1.0)
-, sma(1.0)
-, ecc(0.0)
-, inc(0.0)
-, asc(0.0)
-, arg(0.0)
-, mna(0.0) {
+, orbit()
+, surface_tilt(0.0, 0.0)
+, axis_tilt(0.0, 0.0)
+, rotation(0.0)
+, angular(0.0) {
 }
 
 Body::~Body() {
@@ -74,124 +76,158 @@ void Body::RemoveChild(Body &c) {
        }
 }
 
-double Body::Mass() const noexcept {
-       return mass;
+double Body::Inertia() const noexcept {
+       // assume solid sphere for now
+       return (2.0/5.0) * Mass() * pow(Radius(), 2);
+}
+
+double Body::GravitationalParameter() const noexcept {
+       return G * Mass();
+}
+
+double Body::OrbitalPeriod() const noexcept {
+       if (parent) {
+               return PI_2p0 * sqrt(pow(orbit.SemiMajorAxis(), 3) / (G * (parent->Mass() + Mass())));
+       } else {
+               return 0.0;
+       }
+}
+
+double Body::RotationalPeriod() const noexcept {
+       if (std::abs(angular) < std::numeric_limits<double>::epsilon()) {
+               return std::numeric_limits<double>::infinity();
+       } else {
+               return PI_2p0 * Inertia() / angular;
+       }
+}
+
+glm::dmat4 Body::LocalTransform() const noexcept {
+       glm::dmat4 srf = glm::eulerAngleXY(surface_tilt.x, surface_tilt.y);
+       glm::dmat4 rot = glm::eulerAngleY(rotation);
+       glm::dmat4 tilt = glm::eulerAngleXY(axis_tilt.x, axis_tilt.y);
+       return tilt * rot * srf;
 }
 
-void Body::Mass(double m) noexcept {
-       mass = m;
+glm::dmat4 Body::InverseTransform() const noexcept {
+       glm::dmat4 srf = glm::eulerAngleYX(-surface_tilt.y, -surface_tilt.x);
+       glm::dmat4 rot = glm::eulerAngleY(-rotation);
+       glm::dmat4 tilt = glm::eulerAngleYX(-axis_tilt.y, -axis_tilt.x);
+       return srf * rot * tilt;
+}
+
+glm::dmat4 Body::ToParent() const noexcept {
+       if (!parent) {
+               return glm::dmat4(1.0);
+       }
+       return orbit.InverseMatrix(PI_2p0 * (GetSimulation().Time() / OrbitalPeriod()));
+}
+
+glm::dmat4 Body::FromParent() const noexcept {
+       if (!parent) {
+               return glm::dmat4(1.0);
+       }
+       return orbit.Matrix(PI_2p0 * (GetSimulation().Time() / OrbitalPeriod()));
 }
 
-double Body::Radius() const noexcept {
-       return radius;
+
+Orbit::Orbit()
+: sma(1.0)
+, ecc(0.0)
+, inc(0.0)
+, asc(0.0)
+, arg(0.0)
+, mna(0.0) {
 }
 
-void Body::Radius(double r) noexcept {
-       radius = r;
+Orbit::~Orbit() {
 }
 
-double Body::SemiMajorAxis() const noexcept {
+double Orbit::SemiMajorAxis() const noexcept {
        return sma;
 }
 
-void Body::SemiMajorAxis(double s) noexcept {
+Orbit &Orbit::SemiMajorAxis(double s) noexcept {
        sma = s;
+       return *this;
 }
 
-double Body::Eccentricity() const noexcept {
+double Orbit::Eccentricity() const noexcept {
        return ecc;
 }
 
-void Body::Eccentricity(double e) noexcept {
+Orbit &Orbit::Eccentricity(double e) noexcept {
        ecc = e;
+       return *this;
 }
 
-double Body::Inclination() const noexcept {
+double Orbit::Inclination() const noexcept {
        return inc;
 }
 
-void Body::Inclination(double i) noexcept {
+Orbit &Orbit::Inclination(double i) noexcept {
        inc = i;
+       return *this;
 }
 
-double Body::LongitudeAscending() const noexcept {
+double Orbit::LongitudeAscending() const noexcept {
        return asc;
 }
 
-void Body::LongitudeAscending(double l) noexcept {
+Orbit &Orbit::LongitudeAscending(double l) noexcept {
        asc = l;
+       return *this;
 }
 
-double Body::ArgumentPeriapsis() const noexcept {
+double Orbit::ArgumentPeriapsis() const noexcept {
        return arg;
 }
 
-void Body::ArgumentPeriapsis(double a) noexcept {
+Orbit &Orbit::ArgumentPeriapsis(double a) noexcept {
        arg = a;
+       return *this;
 }
 
-double Body::MeanAnomaly() const noexcept {
+double Orbit::MeanAnomaly() const noexcept {
        return mna;
 }
 
-void Body::MeanAnomaly(double m) noexcept {
+Orbit &Orbit::MeanAnomaly(double m) noexcept {
        mna = m;
+       return *this;
 }
 
-double Body::GravitationalParameter() const noexcept {
-       return G * Mass();
-}
-
-double Body::OrbitalPeriod() const noexcept {
-       if (parent) {
-               return PI_2p0 * sqrt((sma * sma * sma) / (G * (parent->Mass() + Mass())));
-       } else {
-               return 0.0;
-       }
-}
-
-glm::mat4 Body::ToParent() const noexcept {
-       if (!parent) {
-               return glm::mat4(1.0f);
-       }
-
-       double T = OrbitalPeriod();
-
-       double M = mna + PI_2p0 * (GetSimulation().Time() / T); // + time
+namespace {
 
+double mean2eccentric(double M, double e) {
        double E = M; // eccentric anomaly, solve M = E - e sin E
-       while (true) {
-               double dE = (E - ecc * sin(E) - M) / (1 - ecc * cos(E));
+       // limit to 100 steps to prevent deadlocks in impossible situations
+       for (int i = 0; i < 100; ++i) {
+               double dE = (E - e * sin(E) - M) / (1 - e * cos(E));
                E -= dE;
                if (abs(dE) < 1.0e-6) break;
        }
+       return E;
+}
 
-       // coordinates in orbital plane
+}
+
+glm::dmat4 Orbit::Matrix(double t) const noexcept {
+       double M = mna + t;
+       double E = mean2eccentric(M, ecc);
+
+       // coordinates in orbital plane, P=x, Q=-z
        double P = sma * (cos(E) - ecc);
        double Q = sma * sin(E) * sqrt(1 - (ecc * ecc));
 
-       // tile by argument of periapsis, …
-       double x = cos(arg) * P - sin(arg) * Q;
-       double y = sin(arg) * P + cos(arg) * Q;
-       // …inclination, …
-       double z = sin(inc) * x;
-              x = cos(inc) * x;
-       // …and longitude of ascending node
-       glm::vec3 pos(
-               cos(asc) * x - sin(asc) * y,
-               sin(asc) * x + cos(asc) * y,
-               z);
-
-       // TODO: calculate complete matrix
-       return glm::translate(-pos);
+       return glm::translate(glm::yawPitchRoll(asc, inc, arg), glm::dvec3(P, 0.0, -Q));
 }
 
-glm::mat4 Body::FromParent() const noexcept {
-       if (!parent) {
-               return glm::mat4(1.0f);
-       }
-       // TODO: calculate real position
-       return glm::translate(glm::vec3(-sma, 0.0f, 0.0f));
+glm::dmat4 Orbit::InverseMatrix(double t) const noexcept {
+       double M = mna + t;
+       double E = mean2eccentric(M, ecc);
+       double P = sma * (cos(E) - ecc);
+       double Q = sma * sin(E) * sqrt(1 - (ecc * ecc));
+       return glm::transpose(glm::yawPitchRoll(asc, inc, arg)) * glm::translate(glm::dvec3(-P, 0.0, Q));
 }
 
 
@@ -200,6 +236,7 @@ Planet::Planet(int sidelength)
 , sidelength(sidelength)
 , tiles(new Tile[TilesTotal()])
 , vao() {
+       Radius(double(sidelength) / 2.0);
 }
 
 Planet::~Planet() {
diff --git a/tst/assert.hpp b/tst/assert.hpp
new file mode 100644 (file)
index 0000000..85e4324
--- /dev/null
@@ -0,0 +1,37 @@
+#ifndef BLOBS_TEST_ASSETS_HPP_
+#define BLOBS_TEST_ASSETS_HPP_
+
+#include "graphics/glm.hpp"
+
+#include <string>
+#include <limits>
+
+
+namespace blobs {
+namespace test {
+
+template<class T>
+void AssertEqual(
+       const std::string &msg,
+       const glm::tvec3<T> &expected,
+       const glm::tvec3<T> &actual,
+       T epsilon = std::numeric_limits<T>::epsilon()
+) {
+       CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
+               msg + " (X component)",
+               expected.x, actual.x, epsilon
+       );
+       CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
+               msg + " (Y component)",
+               expected.y, actual.y, epsilon
+       );
+       CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
+               msg + " (Z component)",
+               expected.z, actual.z, epsilon
+       );
+}
+
+}
+}
+
+#endif
diff --git a/tst/test.cpp b/tst/test.cpp
new file mode 100644 (file)
index 0000000..ff1ce71
--- /dev/null
@@ -0,0 +1,16 @@
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+using CppUnit::TestFactoryRegistry;
+using CppUnit::TextUi::TestRunner;
+
+
+int main(int, char **) {
+       TestRunner runner;
+       TestFactoryRegistry &registry = TestFactoryRegistry::getRegistry();
+       runner.addTest(registry.makeTest());
+       runner.run();
+
+       return 0;
+
+}
diff --git a/tst/world/OrbitTest.cpp b/tst/world/OrbitTest.cpp
new file mode 100644 (file)
index 0000000..e4c96a0
--- /dev/null
@@ -0,0 +1,77 @@
+#include "OrbitTest.hpp"
+
+#include "../assert.hpp"
+
+#include "const.hpp"
+#include "world/Orbit.hpp"
+
+CPPUNIT_TEST_SUITE_REGISTRATION(blobs::test::world::OrbitTest);
+
+using blobs::world::Orbit;
+
+
+namespace blobs {
+namespace test {
+namespace world {
+
+void OrbitTest::setUp() {
+}
+
+void OrbitTest::tearDown() {
+}
+
+
+void OrbitTest::testSMA() {
+       Orbit orbit;
+       orbit.SemiMajorAxis(1.0);
+       CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
+               "wrong semi-major axis on orbit",
+               1.0, orbit.SemiMajorAxis(), std::numeric_limits<double>::epsilon()
+       );
+
+       // reference direction is +X, so at t=0, the body should be
+       // at (sma,0,0) relative to its parent
+       glm::vec4 pos(orbit.Matrix(0.0) * glm::vec4(0.0f, 0.0f, 0.0f, 1.0f));
+       AssertEqual(
+               "wrong position at t=0",
+               glm::vec3(1.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 90° position should be (0,0,sma) since the zero inclination
+       // reference plane is XZ and rotates counter-clockwise
+       pos = orbit.Matrix(PI_0p5) * glm::vec4(0.0f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=90°",
+               glm::vec3(0.0f, 0.0f, -1.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 180° position should be (-sma,0,0)
+       pos = orbit.Matrix(PI) * glm::vec4(0.0f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=180°",
+               glm::vec3(-1.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 270° position should be (0,0,-sma)
+       pos = orbit.Matrix(PI_1p5) * glm::vec4(0.0f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=270°",
+               glm::vec3(0.0f, 0.0f, 1.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 360° position should be (sma,0,0), the initial position
+       pos = orbit.Matrix(PI_2p0) * glm::vec4(0.0f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=360°",
+               glm::vec3(1.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+}
+
+}
+}
+}
diff --git a/tst/world/OrbitTest.hpp b/tst/world/OrbitTest.hpp
new file mode 100644 (file)
index 0000000..2710a44
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef BLOBS_TEST_WORLD_ORBITTEST_H_
+#define BLOBS_TEST_WORLD_ORBITTEST_H_
+
+#include <cppunit/extensions/HelperMacros.h>
+
+
+namespace blobs {
+namespace test {
+namespace world {
+
+class OrbitTest
+: public CppUnit::TestFixture {
+
+CPPUNIT_TEST_SUITE(OrbitTest);
+
+CPPUNIT_TEST(testSMA);
+
+CPPUNIT_TEST_SUITE_END();
+
+public:
+       void setUp();
+       void tearDown();
+
+       void testSMA();
+
+};
+
+}
+}
+}
+
+#endif