From 412387cc71a74470acdb09ee2f3807488a3d3ed5 Mon Sep 17 00:00:00 2001
From: Quantum <quantum2048@gmail.com>
Date: Fri, 25 Oct 2013 19:14:32 -0400
Subject: [PATCH] Added proper orbit calculations, and limited orbital plane
 calculations (only inclination is used).

---
 punyverse/entity.py  | 18 +++++++++++++-----
 punyverse/orbit.py   | 36 ++++++++++++++++++++++++++++++++++++
 punyverse/world.json |  3 +++
 punyverse/world.py   |  5 +++--
 4 files changed, 55 insertions(+), 7 deletions(-)
 create mode 100644 punyverse/orbit.py

diff --git a/punyverse/entity.py b/punyverse/entity.py
index c3a8c8a..b6d28e9 100644
--- a/punyverse/entity.py
+++ b/punyverse/entity.py
@@ -1,5 +1,7 @@
 from punyverse import framedata
 from math import radians, sin, cos
+from punyverse.orbit import KeplerOrbit
+
 
 class Entity(object):
     def __init__(self, id, location, rotation=(0, 0, 0), direction=(0, 0, 0), background=False):
@@ -48,9 +50,17 @@ class Satellite(Planet):
     def __init__(self, *args, **kwargs):
         self.parent = kwargs.pop('parent')
         self.orbit_speed = kwargs.pop('orbit_speed', 1)
-        self.inclination = kwargs.pop('inclination', 0)
+
+        # Semi-major axis and eccentricity defines orbit
         self.distance = kwargs.pop('distance', 100)
+        self.eccentricity = kwargs.pop('eccentricity', 0)
+
+        # Inclination and longitude of ascending node defines orbital plane
+        self.inclination = kwargs.pop('inclination', 0)
+
         self.theta = 0
+        # OpenGL's z-axis is reversed
+        self.orbit = KeplerOrbit(self.distance, self.eccentricity, -self.inclination)
         super(Satellite, self).__init__(*args, **kwargs)
 
     def update(self):
@@ -63,10 +73,8 @@ class Satellite(Planet):
             self.rotation = pitch, yaw, roll
 
             self.parent.update()
-            x, y, z = self.location
             px, py, pz = self.parent.location
             self.theta += self.orbit_speed
-            x = cos(radians(self.theta)) * self.distance + px
-            z = sin(radians(self.theta)) * self.distance + pz
-            self.location = (x, y, z)
+            x, z, y = self.orbit.orbit(self.theta)
+            self.location = (x + px, y + py, z + pz)
 
diff --git a/punyverse/orbit.py b/punyverse/orbit.py
new file mode 100644
index 0000000..fb21589
--- /dev/null
+++ b/punyverse/orbit.py
@@ -0,0 +1,36 @@
+from math import sin, cos, tan, atan, sqrt, radians
+
+
+class KeplerOrbit(object):
+    def __init__(self, sma, eccentricity, inclination=0):
+        self.sma = sma
+        self.eccentricity = eccentricity
+        self.inclination = radians(inclination)
+        self.__true_anomaly_factor = sqrt((1 + eccentricity)/(1 - eccentricity))
+        self.__distance_factor = sma * (1 - eccentricity ** 2)
+        self.__sin_inclination = sin(self.inclination)
+        self.__cos_inclination = cos(self.inclination)
+
+    def eccentric_anomaly(self, mean_anomaly):
+        e1 = mean_anomaly
+        e2 = mean_anomaly + self.eccentricity * sin(e1)
+        while abs(e1 - e2) > 0.000001:
+            e1, e2 = e2, mean_anomaly + self.eccentricity * sin(e2)
+        return e2
+
+    def true_anomaly(self, mean_anomaly):
+        eccentric_anomaly = self.eccentric_anomaly(mean_anomaly)
+        return 2 * atan(self.__true_anomaly_factor * tan(eccentric_anomaly / 2))
+
+    def orbit(self, mean_anomaly):
+        mean_anomaly = radians(mean_anomaly)
+        phi = self.true_anomaly(mean_anomaly)
+        r = self.__distance_factor / (1 + self.eccentricity * cos(phi))
+        x = r * cos(phi)
+        y = r * sin(phi)
+        z = 0
+
+        x, z = (x * self.__cos_inclination - z * self.__sin_inclination,
+                x * self.__sin_inclination + z * self.__cos_inclination)
+
+        return x, y, z
diff --git a/punyverse/world.json b/punyverse/world.json
index 130dc56..c11a798 100644
--- a/punyverse/world.json
+++ b/punyverse/world.json
@@ -17,6 +17,9 @@
           "texture": ["moon.jpg", "moon_medium.jpg", "moon_small.jpg", [0.53, 0.53, 0.53, 1]],
           "radius": 25,
           "distance": 200,
+          "orbit_speed": 1,
+          "eccentricity": 0.0549,
+          "inclination": 5.145,
           "pitch": -90,
           "yaw": 6.68
         },
diff --git a/punyverse/world.py b/punyverse/world.py
index e25bd38..4346abf 100644
--- a/punyverse/world.py
+++ b/punyverse/world.py
@@ -107,8 +107,9 @@ def load_world(file):
                 x, y, z = parent.location
                 distance = e(info.get('distance', 100))
                 x -= distance
-                type = partial(Satellite, parent=parent, distance=distance, inclination=e(info.get('inclination', 0)),
-                               orbit_speed=e(info.get('orbit_speed', 1)))
+                type = partial(Satellite, parent=parent, orbit_speed=info.get('orbit_speed', 1),
+                               distance=distance, eccentricity=info.get('eccentricity', 0),
+                               inclination=info.get('inclination', 0))
 
             atmosphere_id = 0
             cloudmap_id = 0