diff --git a/bodies/body.py b/bodies/body.py index 0803c4135dc94cf1257fb906741b9762b438beed..ebecbdd5286a96946aaea4c7566d33b37f76dd90 100644 --- a/bodies/body.py +++ b/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): """ diff --git a/bodies/earth.py b/bodies/earth.py index b87e94a2887e6472b0512334ddf54159693a8664..ce68f82b4b20f1717c6649c89674cdbd6d9fc168 100644 --- a/bodies/earth.py +++ b/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, diff --git a/bodies/sun.py b/bodies/sun.py index d698b3624b58a69befb6417627adc54c554636bc..daa29ec21b0a1d016ef8c2e4b2a5e7b2093d8d6b 100644 --- a/bodies/sun.py +++ b/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, diff --git a/common/simulator.py b/common/simulator.py index 77c2c07d03699788a91535bf8fd8216cebcd5f51..af608f875c8d5183cd8647d0de62238a8d5b80b6 100644 --- a/common/simulator.py +++ b/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) diff --git a/common/system.py b/common/system.py index d0ee4965f3776797a678977d484c0317921e7340..88f6480224eae9c57b8f26664c8f89e2d613bd9f 100644 --- a/common/system.py +++ b/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 """