From c5556cf5f6813887a3503433c021ccd2e7fae865 Mon Sep 17 00:00:00 2001
From: Daniel Karbach <daniel.karbach@localhorst.tv>
Date: Sat, 11 Nov 2017 14:03:23 +0100
Subject: [PATCH] more orbits and stuff

---
 src/app/Assets.hpp          |   2 +
 src/app/MasterState.hpp     |   3 +
 src/app/states.cpp          |  36 ++++++--
 src/blobs.cpp               |  33 +++++--
 src/graphics/SunSurface.hpp |  68 ++++++++++++++
 src/graphics/shader.cpp     | 149 ++++++++++++++++++++++++++++++
 src/world/Body.hpp          |  55 ++++++------
 src/world/Orbit.hpp         |  53 +++++++++++
 src/world/Simulation.hpp    |   5 ++
 src/world/sim.cpp           |  17 +++-
 src/world/world.cpp         | 175 ++++++++++++++++++++++--------------
 tst/assert.hpp              |  37 ++++++++
 tst/test.cpp                |  16 ++++
 tst/world/OrbitTest.cpp     |  77 ++++++++++++++++
 tst/world/OrbitTest.hpp     |  32 +++++++
 15 files changed, 648 insertions(+), 110 deletions(-)
 create mode 100644 src/graphics/SunSurface.hpp
 create mode 100644 src/world/Orbit.hpp
 create mode 100644 tst/assert.hpp
 create mode 100644 tst/test.cpp
 create mode 100644 tst/world/OrbitTest.cpp
 create mode 100644 tst/world/OrbitTest.hpp

diff --git a/src/app/Assets.hpp b/src/app/Assets.hpp
index ed3c2e2..c80726a 100644
--- a/src/app/Assets.hpp
+++ b/src/app/Assets.hpp
@@ -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();
diff --git a/src/app/MasterState.hpp b/src/app/MasterState.hpp
index fb245cb..f62c773 100644
--- a/src/app/MasterState.hpp
+++ b/src/app/MasterState.hpp
@@ -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;
 
 };
 
diff --git a/src/app/states.cpp b/src/app/states.cpp
index ef8b514..67c073a 100644
--- a/src/app/states.cpp
+++ b/src/app/states.cpp
@@ -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();
 }
 
 }
diff --git a/src/blobs.cpp b/src/blobs.cpp
index 5f7bb76..c302159 100644
--- a/src/blobs.cpp
+++ b/src/blobs.cpp
@@ -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
index 0000000..5ce91b1
--- /dev/null
+++ b/src/graphics/SunSurface.hpp
@@ -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
diff --git a/src/graphics/shader.cpp b/src/graphics/shader.cpp
index b184bb9..d771d9f 100644
--- a/src/graphics/shader.cpp
+++ b/src/graphics/shader.cpp
@@ -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);
+}
+
 }
 }
diff --git a/src/world/Body.hpp b/src/world/Body.hpp
index 2f87d22..1b9e2e3 100644
--- a/src/world/Body.hpp
+++ b/src/world/Body.hpp
@@ -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
index 0000000..7def521
--- /dev/null
+++ b/src/world/Orbit.hpp
@@ -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
diff --git a/src/world/Simulation.hpp b/src/world/Simulation.hpp
index 6d22ed3..3025a74 100644
--- a/src/world/Simulation.hpp
+++ b/src/world/Simulation.hpp
@@ -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;
 
 };
diff --git a/src/world/sim.cpp b/src/world/sim.cpp
index abf2829..46454fd 100644
--- a/src/world/sim.cpp
+++ b/src/world/sim.cpp
@@ -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());
+	}
 }
 
 }
diff --git a/src/world/world.cpp b/src/world/world.cpp
index 5c11350..0b1b22d 100644
--- a/src/world/world.cpp
+++ b/src/world/world.cpp
@@ -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
index 0000000..85e4324
--- /dev/null
+++ b/tst/assert.hpp
@@ -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
index 0000000..ff1ce71
--- /dev/null
+++ b/tst/test.cpp
@@ -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
index 0000000..e4c96a0
--- /dev/null
+++ b/tst/world/OrbitTest.cpp
@@ -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
index 0000000..2710a44
--- /dev/null
+++ b/tst/world/OrbitTest.hpp
@@ -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
-- 
2.39.2