三体运行模拟器

bodies/body.py

@@ -32,14 +32,19 @@ class Body:
         :param size_scale: 尺寸缩放
         :param distance_scale: 距离缩放
         """
+        self.__his_pos = []
+        self.__his_vel = []
+        self.__his_acc = []
+        self.__his_reserved_num = 20
+
         self.name = name
         self.__mass = mass
 
         self.init_position = np.array(init_position, dtype='float32')
         self.init_velocity = np.array(init_velocity, dtype='float32')
 
-        self.position = self.init_position
-        self.velocity = self.init_velocity
+        self.__position = self.init_position
+        self.__velocity = self.init_velocity
 
         self.__density = density
 
@@ -50,14 +55,40 @@ class Body:
         self.distance_scale = distance_scale
 
         # 初始化后，加速度为0，只有多个天体的引力才会影响到加速度
-        self.acceleration = np.array([0, 0, 0], dtype='float32')
+        # m/s²
+        self.__acceleration = np.array([0, 0, 0], dtype='float32')
+        self.__record_history()
+
 
-        self.__his_pos = []
-        self.__his_vel = []
-        self.__his_acc = []
-        self.__his_reserved_num = 10
 
-    def append_history(self, his_list, data):
+    @property
+    def position(self):
+        return self.__position
+
+    @position.setter
+    def position(self, value):
+        self.__position = value
+        self.__record_history()
+
+    @property
+    def acceleration(self):
+        return self.__acceleration
+
+    @acceleration.setter
+    def acceleration(self, value):
+        self.__acceleration = value
+        self.__record_history()
+
+    @property
+    def velocity(self):
+        return self.__velocity
+
+    @velocity.setter
+    def velocity(self, value):
+        self.__velocity = value
+        self.__record_history()
+
+    def __append_history(self, his_list, data):
         """
 
         :param his_list:
@@ -69,23 +100,22 @@ class Body:
                 np.sum(data == his_list[-1]) < len(data):
             his_list.append(data.copy())
 
-    def record_history(self):
+    def __record_history(self):
         """
         记录历史
         :return:
         """
         # 如果历史记录数超过了保留数量，则截断，只保留 __his_reserved_num 数量的历史
         if len(self.__his_pos) > self.__his_reserved_num:
-            self.__his_pos = self.__his_pos[len(self.__his_pos)-self.__his_reserved_num:]
-            self.__his_vel = self.__his_vel[len(self.__his_vel)-self.__his_reserved_num:]
-            self.__his_acc = self.__his_acc[len(self.__his_acc)-self.__his_reserved_num:]
+            self.__his_pos = self.__his_pos[len(self.__his_pos) - self.__his_reserved_num:]
+            self.__his_vel = self.__his_vel[len(self.__his_vel) - self.__his_reserved_num:]
+            self.__his_acc = self.__his_acc[len(self.__his_acc) - self.__his_reserved_num:]
 
         # 追加历史记录(位置、速度、加速度)
-        self.append_history(self.__his_pos, self.position)
-        self.append_history(self.__his_vel, self.velocity)
-        self.append_history(self.__his_acc, self.acceleration)
-
-        print(self.name, "his pos->", self.__his_pos)
+        self.__append_history(self.__his_pos, self.position)
+        self.__append_history(self.__his_vel, self.velocity)
+        self.__append_history(self.__his_acc, self.acceleration)
+        # print(self.name, "his pos->", self.__his_pos)
 
     def his_position(self):
         """

bodies/earth.py

@@ -23,7 +23,7 @@ class Earth(Body):
 
     def __init__(self, name="earth", mass=5.97237e24,
                  init_position=[149597870.700, 0, 0],
-                 init_velocity=[29.79, 0, 0],
+                 init_velocity=[0, 29.79, 0],
                  texture="", size_scale=1.0, distance_scale=1.0):
         params = {
             "name": name,

bodies/sun.py

@@ -10,7 +10,9 @@ from bodies.body import Body
 
 
 class Sun(Body):
-    def __init__(self, name="sun", mass=1.9891e30, init_position=[0, 0, 0], init_velocity=[0, 0, 0],
+    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):
         params = {
             "name": name,

common/simulator.py

@@ -17,27 +17,47 @@ X_MIN, X_MAX = -2e+14, 2e+14  # the x range of the bounding box (m)
 Y_MIN, Y_MAX = -2e+14, 2e+14  # the y range of the bounding box (m)
 Z_MIN, Z_MAX = -2e+14, 2e+14  # the z range of the bounding box (m)
 
+X_MIN, X_MAX = -2e+12, 2e+12  # the x range of the bounding box (m)
+Y_MIN, Y_MAX = -2e+12, 2e+12  # the y range of the bounding box (m)
+Z_MIN, Z_MAX = -2e+12, 2e+12  # the z range of the bounding box (m)
+
+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)
+
 
 def show(bodies, idx=0):
+    # from scipy.interpolate import make_interp_spline
+
     # Creating figures for the plot
     fig = plt.figure(figsize=(16, 12))
     ax = plt.axes(projection="3d")
     ax.set_xlim(X_MIN, X_MAX)
     ax.set_ylim(Y_MIN, Y_MAX)
     ax.set_zlim(Z_MIN, Z_MAX)
+    ax.set_xlabel('X')
+    ax.set_ylabel('Y')
+    ax.set_zlabel('Z')
     # Creating a plot using the random datasets
-    colors = ['red', 'blue']
-    sizes = [800, 500]
+    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
         pos = body.position
-        ax.scatter(pos[0], pos[1], pos[2], color=colors[idx], s=sizes[idx])
-        for _his_pos in body.his_position():
-            ax.scatter3D(_his_pos[0], _his_pos[1], _his_pos[2], color=colors[idx], alpha=0.5)
-            # ax.scatter(his_pos[0], his_pos[1], his_pos[2], color=colors[idx], s=10)
+        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)
+        # for _his_pos in body.his_position():
+        #     ax.scatter3D(_his_pos[0], _his_pos[1], _his_pos[2], color=color, alpha=0.5)
+        # ax.scatter(his_pos[0], his_pos[1], his_pos[2], color=colors[idx], s=10)
         his_pos = body.his_position()
-        if len(his_pos) > 2:
+        if len(his_pos) > 1:
             _his_pos = list(zip(*his_pos))
-            ax.plot3D(_his_pos[0], _his_pos[1], _his_pos[2], color=colors[idx], alpha=0.5)
+            ax.plot3D(_his_pos[0], _his_pos[1], _his_pos[2], color=color, alpha=0.8)
     plt.title("3D scatter plot %s" % idx)
     # display the  plot
     plt.show()
@@ -52,6 +72,7 @@ class Simulator:
         self.dt = dt
 
     def evolve(self, dt):
+        show(self.system.bodies)
         self.system.update_acceleration()
 
         # next_p = p + dt * v + 0.5 * a * (dt ** 2)
@@ -68,30 +89,42 @@ class Simulator:
         #         body1.update_source_data()
 
         for body in self.system.bodies:
-            body.position += body.velocity * dt + 0.5 * body.acceleration * (dt ** 2)
+            # body.vx = body.vx + ax * dt
+            # body.x = body.x + body.vx * dt
             # acceleration 加速度
             body.velocity += body.acceleration * dt
-
-            print(body.name, body.position)
-            body.record_history()
-
-        show(self.system.bodies)
+            # body.position += body.velocity * dt  # - 0.5 * body.acceleration * (dt ** 2)
+            body.position += body.velocity * dt
+            print(body)
+            # print(body.name, body.position)
 
     def run(self, t):
         n = int(t / self.dt)
         for i in range(n):
-            time.sleep(1)
             self.evolve(self.dt)
+            time.sleep(1)
 
     def run_always(self):
         dt = 1  # 时间变化1秒
         while True:
-            time.sleep(1)
+            # time.sleep(1)
             self.evolve(dt)
 
 
 if __name__ == '__main__':
-    t = 60 * 60 * 24
-    _sys = System([Sun(), Earth()])
+    t = 60 * 60 * 24 * 100
+
+    # _sys = System([Sun(init_position=[0, 0, 149597870.700]), Earth()])
+
+    # _sys = System([Sun(init_position=[849597870.700, 0, 0], init_velocity=[0, 9.79, 0]),
+    #                Sun(init_position=[0, 0, 0], init_velocity=[0, -9.79, 0])])
+
+    # _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])])
+
     _sim = Simulator(_sys, t)
     _sim.run(t * 100)

common/system.py

@@ -139,12 +139,30 @@ class System(object):
     #         body1.position = body1.position + (body1.momentum / body1.mass) * dt
     #         body1.update_source_data()
 
+
+    def acceleration(self, body1, body2, G=6.67e-11):
+        """
+        计算两个天体之间的加速度
+        :param body1: 天体1
+        :param body2: 天体2
+        :param G: 引力常数
+        :return: 加速度
+        """
+        import math
+        dx = body2.position[0] - body1.position[0]
+        dy = body2.position[1] - body1.position[1]
+        dz = body2.position[2] - body1.position[2]
+        distance = math.sqrt(dx ** 2 + dy ** 2 + dz ** 2)
+        mag = G * body2.mass / (distance ** 2)
+        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)
+        # for body in self.bodies:
+        #     body.acceleration = np.zeros(3)
 
         for body1 in self.bodies:
-            body1.record_history()
+            acceleration = np.zeros(3)
             for body2 in self.bodies:
                 if body1 is body2:
                     continue
@@ -154,7 +172,9 @@ class System(object):
                 # body.acceleration1 += aij * dir
 
                 r = body2.position - body1.position
-                body.acceleration += -6.67e-11 * body2.mass * r / np.linalg.norm(r) ** 3
+                # m/s²
+                # body1.acceleration += self.acceleration(body1,body2)
+                acceleration += (6.67e-11 * body2.mass * r / np.linalg.norm(r) ** 3) / 1e9
 
                 # body1.acceleration = body.acceleration2
                 # r = body2.position - body1.position
@@ -168,6 +188,7 @@ class System(object):
                 # r = 1 # 两个天体之间的距离
                 #
                 # a = G * m2 / math.pow(r, 2)
+            body1.acceleration = acceleration
 
 
 """