三体运行模拟器

bodies/__init__.py

 from bodies.body import Body
 from bodies.earth import Earth
-from bodies.sun import Sun
\ No newline at end of file
+from bodies.sun import Sun
+from bodies.jupiter import Jupiter
\ No newline at end of file

bodies/body.py

@@ -35,7 +35,7 @@ class Body:
         self.__his_pos = []
         self.__his_vel = []
         self.__his_acc = []
-        self.__his_reserved_num = 20
+        self.__his_reserved_num = 100
 
         self.name = name
         self.__mass = mass

bodies/earth.py

@@ -6,7 +6,7 @@
 # notes           :
 # python_version  :3.8
 # ==============================================================================
-from bodies.body import Body
+from bodies.body import Body, AU
 
 
 class Earth(Body):
@@ -21,8 +21,8 @@ class Earth(Body):
      　平均密度: 5507.85 kg/m³
     """
 
-    def __init__(self, name="earth", mass=5.97237e24,
-                 init_position=[149597870.700, 0, 0],
+    def __init__(self, name="Earth", mass=5.97237e24,
+                 init_position=[1.12 * AU, 0, 0],
                  init_velocity=[0, 29.79, 0],
                  texture="", size_scale=1.0, distance_scale=1.0):
         params = {

bodies/jupiter.py

+# -*- coding:utf-8 -*-
+# title           :
+# description     :
+# author          :Python超人
+# date            :2023-01-22
+# notes           :
+# python_version  :3.8
+# ==============================================================================
+from bodies.body import Body, AU
+
+
+class Jupiter(Body):
+    """
+    木星
+    ------------------------
+    远日点距离: 5.4588 天文单位
+    近日点距离: 4.9501 天文单位
+    　逃逸速度: 59.5 km/s
+    　公转速度: 13.06 km/s
+    　天体质量: 1.8982✕10²⁷ kg(317.8 M⊕)
+    　平均密度: 1.326 g/cm³ -> -> 1.326✕10³ kg/m³
+    """
+
+    def __init__(self, name="Jupiter", mass=1.8982e27,
+                 init_position=[5.2 * AU, 0, 0],
+                 init_velocity=[0, 13.06, 0],
+                 texture="", size_scale=1.0, distance_scale=1.0):
+        params = {
+            "name": name,
+            "mass": mass,
+            "init_position": init_position,
+            "init_velocity": init_velocity,
+            "density": 1.326e3,
+            "color": [
+                125,
+                125,
+                125
+            ],
+            "texture": texture,
+            "size_scale": size_scale,
+            "distance_scale": distance_scale
+        }
+        super().__init__(**params)
+
+
+if __name__ == '__main__':
+    jupiter = Jupiter()
+    print(jupiter)

bodies/sun.py

@@ -10,7 +10,7 @@ from bodies.body import Body
 
 
 class Sun(Body):
-    def __init__(self, name="sun", mass=1.9891e30,
+    def __init__(self, name="Sun", mass=1.9891e30,
                  init_position=[0, 0, 0],
                  init_velocity=[0, 0, 0],
                  texture="", size_scale=1.0, distance_scale=1.0):

common/simulator.py

@@ -7,7 +7,7 @@
 # python_version  :3.8
 # ==============================================================================
 from common.system import System
-from bodies import Body, Sun, Earth
+from bodies import Body, Sun, Earth, Jupiter
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -25,6 +25,7 @@ X_MIN, X_MAX = -1e+9, 1e+9  # the x range of the bounding box (m)
 Y_MIN, Y_MAX = -1e+9, 1e+9  # the y range of the bounding box (m)
 Z_MIN, Z_MAX = -1e+9, 1e+9  # the z range of the bounding box (m)
 
+
 # X_MIN, X_MAX = -8e+8, 8e+8  # the x range of the bounding box (m)
 # Y_MIN, Y_MAX = -8e+8, 8e+8  # the y range of the bounding box (m)
 # Z_MIN, Z_MAX = -8e+8, 8e+8  # the z range of the bounding box (m)
@@ -46,8 +47,8 @@ def show(bodies, idx=0):
     colors = ['red', 'blue', 'red', 'red']
     sizes = [800, 500, 800, 800]
     for idx, body in enumerate(bodies):
-        color = 'red' if str(body.name).startswith("sun") else 'blue'
-        size = 800 if str(body.name).startswith("sun") else 500
+        color = 'red' if str(body.name).lower().startswith("sun") else 'blue'
+        size = 800 if str(body.name).lower().startswith("sun") else 500
         pos = body.position
         ax.text(pos[0], pos[1], pos[2] + 1e8, body.name, color=color, fontsize=20)
         ax.scatter(pos[0], pos[1], pos[2], color=color, s=size)
@@ -73,28 +74,13 @@ class Simulator:
 
     def evolve(self, dt):
         show(self.system.bodies)
-        self.system.update_acceleration()
-
-        # next_p = p + dt * v + 0.5 * a * (dt ** 2)
-        # next_v = v + dt * a
-        #     for body1 in bodies:
-        #         for body2 in bodies:
-        #             if body1 == body2:  # 相等说明是同一个天体，跳过计算
-        #                 continue
-        #             r = body2.position - body1.position
-        #             # F = G x (m1 x m2) / r²  万有引力定律
-        #             F = Configs.G * body1.mass * body2.mass * r / np.linalg.norm(r) / r.dot(r)
-        #             body1.momentum = body1.momentum + F * dt  # 动量定理
-        #         body1.position = body1.position + (body1.momentum / body1.mass) * dt
-        #         body1.update_source_data()
+        self.system.calc_acceleration()
 
         for body in self.system.bodies:
-            # body.vx = body.vx + ax * dt
-            # body.x = body.x + body.vx * dt
             # acceleration 加速度
             body.velocity += body.acceleration * dt
-            # body.position += body.velocity * dt  # - 0.5 * body.acceleration * (dt ** 2)
-            body.position += body.velocity * dt
+            # body.position += 0.5 * body.acceleration * (dt ** 2)
+            body.position += 0.5 * body.velocity * dt
             print(body)
             # print(body.name, body.position)
 
@@ -121,10 +107,12 @@ if __name__ == '__main__':
 
     # _sys = System([Sun(), Earth()])
 
-    _sys = System([Sun(name="sun1", init_position=[849597870.700, 0, 0], init_velocity=[0, 7.0, 0]),
-                   Sun(name="sun2", init_position=[0, 0, 0], init_velocity=[0, -8.0, 0]),
-                   Sun(name="sun3", init_position=[0, -849597870.700, 0], init_velocity=[18.0, 0, 0]),
-                   Earth(init_position=[0, -349597870.700, 0], init_velocity=[15.50, 0, 0])])
+    # _sys = System([Sun(name="sun1", init_position=[849597870.700, 0, 0], init_velocity=[0, 7.0, 0]),
+    #                Sun(name="sun2", init_position=[0, 0, 0], init_velocity=[0, -8.0, 0]),
+    #                Sun(name="sun3", init_position=[0, -849597870.700, 0], init_velocity=[18.0, 0, 0]),
+    #                Earth(init_position=[0, -349597870.700, 0], init_velocity=[15.50, 0, 0])])
+
+    _sys = System([Sun(), Earth(), Jupiter()])
 
     _sim = Simulator(_sys, t)
     _sim.run(t * 100)

common/system.py

@@ -157,37 +157,21 @@ class System(object):
         acc = [dx / distance * mag, dy / distance * mag, dz / distance * mag]
         return np.array(acc)/1000/1000/1000
 
-    def update_acceleration(self):
-        # for body in self.bodies:
-        #     body.acceleration = np.zeros(3)
-
+    def calc_acceleration(self):
         for body1 in self.bodies:
             acceleration = np.zeros(3)
             for body2 in self.bodies:
                 if body1 is body2:
                     continue
-                # r = np.sqrt(np.sum(np.square(body2.position - body1.position)))
-                # aij = 6.67e-11 * body2.mass / (r * r)
-                # dir = (body2.position - body1.position) / r
-                # body.acceleration1 += aij * dir
 
                 r = body2.position - body1.position
-                # m/s²
-                # body1.acceleration += self.acceleration(body1,body2)
+                #
+                # G = 6.67e-11 # 万有引力常数
+                # m/s² = kg * m / m**3
+                # km/s² = kg * m / m**3 / 1e9
+                # acceleration = G * body2.mass * dx / (d ** 3)
                 acceleration += (6.67e-11 * body2.mass * r / np.linalg.norm(r) ** 3) / 1e9
 
-                # body1.acceleration = body.acceleration2
-                # r = body2.position - body1.position
-                # # body1.acceleration = -body2.mass * r / np.linalg.norm(r) ** 3
-                # body1.acceleration += 6.67e-11 * body2.mass * r / np.linalg.norm(r) ** 3
-                # body1.acceleration = (G * body2.mass * (x / pow(x ** 2 + y ** 2 + z ** 2, 3 / 2)))
-                # pow(x ** 2 + y ** 2 + z ** 2, 3 / 2)
-                # G = 6.67e-11 # 万有引力常数
-                # m1 = 1 # 第一个天体的质量
-                # m2 = 1 # 第二个天体的质量
-                # r = 1 # 两个天体之间的距离
-                #
-                # a = G * m2 / math.pow(r, 2)
             body1.acceleration = acceleration