三体世界的模拟

  1 ‘‘‘三体问题求解及可视化,去掉了动图模块‘‘‘
  2 ‘‘‘Coworker：聂嘉颖，张瀚文‘‘‘
  3 import numpy as np
  4 from numpy import arange
  5 import matplotlib.pyplot as plt
  6 from mpl_toolkits.mplot3d import Axes3D
  7 from math import sqrt
  8 
  9 
 10 # 得到太阳对行星夹角的cot值
 11 def sun_height(x0, y0, z0, x1, y1, z1):
 12     a0 = x1 - x0
 13     b0 = y1 - y0
 14     c0 = z1 - z0
 15     return a0 / sqrt(b0 ** 2 + c0 ** 2)
 16 
 17 
 18 # 得到两星体间的距离
 19 def distants(x0, y0, z0, x1, y1, z1):
 20     return sqrt((x0 - x1) ** 2 + (y0 - y1) ** 2 + (z0 - z1) ** 2)
 21 
 22 
 23 # 计算某一瞬时某一星体对行星的辐射强度
 24 def sun_heat(d, sun_temperature, x0, y0, z0, x1, y1, z1):
 25     theta = d/distants(x0, y0, z0, x1, y1, z1)
 26     earth_temperature = 0.5 * sqrt(theta) * sun_temperature
 27     # print(earth_temperature)
 28     return earth_temperature
 29 
 30 # 定义星体质量
 31 # star_4选定为行星，其余为恒星
 32 star_1_mass = 3e30  # kg
 33 star_2_mass = 2e30  # kg
 34 star_3_mass = 3e30  # kg
 35 star_4_mass = 2e30  # kg
 36 
 37 diameter0 = 1.0392e10  # 米
 38 diameter1 = 1.0392e10  # 米
 39 diameter2 = 1.00392e11  # 米
 40 
 41 # 定义恒星表面温度
 42 sun_temperature0 = 60.  # K
 43 sun_temperature1 = 600.  # K
 44 sun_temperature2 = 60.  # K
 45 
 46 # 引力常数
 47 gravitational_constant = 6.67e-11  # m3 / kg s2
 48 
 49 # 行星运动总时长（按地球年计算）
 50 # 每两小时计算一次星体位置
 51 end_time = 16 * 365.26 * 24. * 3600.  # s
 52 h = 2. * 3600  # s
 53 num_steps = int(end_time / h)
 54 tpoints = arange(0, end_time, h)
 55 
 56 
 57 def three_body_problem():
 58     ‘‘‘主程序，计算星体位置和表面温度‘‘‘
 59     positions = np.zeros([num_steps + 1, 4, 3])  # m
 60     velocities = np.zeros([num_steps + 1, 4, 3])  # m / s
 61 
 62     positions[0] = np.array([[1., 3., 2.], [6., -5., 4.], [7., 8., -7.], [7., -2., 9.]]) * 1e11  # m
 63     velocities[0] = np.array([[-2., 0.5, 5.], [7., 0.5, 2.], [4., -0.5, 3.], [-10., 4., -2.]]) * 1e3  # m/s
 64     positions[0] = np.array([[1., 3., 2.], [3., -5., 1.], [2., 8., -7.], [-3., -2., 9.]]) * 1e11  # m
 65     velocities[0] = np.array([[0, 0., 0.], [0., 0., 0.], [0., 0., 0.], [0., 0., 0.]]) * 1e3  # m/s
 66 
 67     def acceleration(positions):
 68         a = np.zeros([4, 3])
 69         a[0] = gravitational_constant * star_2_mass / np.linalg.norm(positions[0] - positions[1]) ** 3 * (
 70             positions[1] - positions[0]) + gravitational_constant * star_3_mass / np.linalg.norm(
 71             positions[0] - positions[2]) ** 3 * (positions[2] - positions[
 72             0]) + gravitational_constant * star_4_mass / np.linalg.norm(positions[0] - positions[3]) ** 3 * (
 73         positions[3] - positions[0])
 74         a[1] = gravitational_constant * star_1_mass / np.linalg.norm(positions[1] - positions[0]) ** 3 * (
 75             positions[0] - positions[1]) + gravitational_constant * star_3_mass / np.linalg.norm(
 76             positions[1] - positions[2]) ** 3 * (positions[2] - positions[
 77             1]) + gravitational_constant * star_4_mass / np.linalg.norm(positions[1] - positions[3]) ** 3 * (
 78         positions[3] - positions[1])
 79         a[2] = gravitational_constant * star_1_mass / np.linalg.norm(positions[2] - positions[0]) ** 3 * (
 80             positions[0] - positions[2]) + gravitational_constant * star_2_mass / np.linalg.norm(
 81             positions[2] - positions[1]) ** 3 * (positions[1] - positions[
 82             2]) + gravitational_constant * star_4_mass / np.linalg.norm(positions[2] - positions[3]) ** 3 * (
 83         positions[3] - positions[2])
 84         a[3] = gravitational_constant * star_1_mass / np.linalg.norm(positions[3] - positions[0]) ** 3 * (
 85         positions[0] - positions[3]) + gravitational_constant * star_2_mass / np.linalg.norm(
 86             positions[3] - positions[1]) ** 3 * (positions[1] - positions[
 87             3]) + gravitational_constant * star_3_mass / np.linalg.norm(positions[3] - positions[2]) ** 3 * (
 88         positions[2] - positions[3])
 89         return a
 90 
 91     def velocity(p, t, velo):
 92         v = np.zeros([4, 3])
 93         v[0] = acceleration(p)[0] * t + velo[0]
 94         v[1] = acceleration(p)[1] * t + velo[1]
 95         v[2] = acceleration(p)[2] * t + velo[2]
 96         v[3] = acceleration(p)[3] * t + velo[3]
 97         return velocities[step]
 98 
 99     heat_sum, t = [0.1], [0]
100 
101     per_0 = 0
102     for step in range(num_steps):
103         # 四阶龙格库塔法求星体位置
104         j1 = h * velocity(positions[step], h, velocities[step])
105         j2 = h * velocity(positions[step] + 0.5 * j1, h + 0.5 * h, velocities[step])
106         j3 = h * velocity(positions[step] + 0.5 * j2, h + 0.5 * h, velocities[step])
107         j4 = h * velocity(positions[step] + j3, h + h, velocities[step])
108         positions[step + 1] = positions[step] + (j1 + 2 * j2 + 2 * j3 + j4) / 6
109         velocities[step + 1] = velocities[step] + h * acceleration(positions[step + 1])
110 
111         # 从 positions 中取出坐标值
112         x0, y0, z0 = positions[step, 0, 0], positions[step, 0, 1], positions[step, 0, 2]
113         x1, y1, z1 = positions[step, 1, 0], positions[step, 1, 1], positions[step, 1, 2]
114         x2, y2, z2 = positions[step, 2, 0], positions[step, 2, 1], positions[step, 2, 2]
115         x3, y3, z3 = positions[step, 3, 0], positions[step, 3, 1], positions[step, 3, 2]
116 
117         # 计算三个太阳对行星表面的总辐射强度
118         heat0 = sun_heat(diameter0, sun_temperature0, x0, y0, z0, x3, y3, z3)
119         heat1 = sun_heat(diameter1, sun_temperature1, x1, y1, z1, x3, y3, z3)
120         heat2 = sun_heat(diameter2, sun_temperature2, x2, y2, z2, x3, y3, z3)
121         heat_sum.append(heat0 + heat1 + heat2)
122 
123         per = int(100 * step / num_steps)
124         if per > per_0:
125             print(per, ‘%‘)
126             per_0 = per
127 
128         t.append(step)
129 
130     return positions, t, heat_sum
131 
132 positions, t, heat_sum = three_body_problem()
133 
134 
135 empty, extinction_line, frozen_line = [], [], []
136 
137 for i in range(len(t)):
138     empty.append(0)
139     extinction_line.append(100)
140     frozen_line.append(40)
141 
142 
143 fig, ax = plt.subplots()
144 fig.set_tight_layout(True)
145 plt.plot(t, heat_sum, ‘b‘)
146 plt.plot(t, empty, ‘g‘)
147 plt.plot(t, extinction_line, ‘r‘)
148 plt.plot(t, frozen_line, ‘r‘)
149 
150 ax.set_xlabel(‘X‘)
151 ax.set_xlim(0, len(t)+10e3)
152 ax.set_ylabel(‘Y‘)
153 plt.show()
154 
155 print("\033[1;31;47m---Begin ploting image---\033[0m")
156 
157 fig = plt.figure()
158 ax2 = Axes3D(fig)
159 
160 ax2.plot(positions[:, 0, 0], positions[:, 0, 1], positions[:, 0, 2], color=‘b‘)
161 ax2.plot(positions[:, 1, 0], positions[:, 1, 1], positions[:, 1, 2], color=‘y‘)
162 ax2.plot(positions[:, 2, 0], positions[:, 2, 1], positions[:, 2, 2], color=‘b‘)
163 ax2.plot(positions[:, 3, 0], positions[:, 3, 1], positions[:, 3, 2], color=‘r‘)
164 
165 zoom = 1.2e12
166 ax2.set_xlabel(‘X‘)
167 # ax2.set_xlim3d(-zoom, zoom + 3.e12)
168 ax2.set_ylabel(‘Y‘)
169 # ax2.set_ylim3d(-zoom, zoom)
170 ax2.set_zlabel(‘Z‘)
171 # ax2.set_zlim3d(-zoom, zoom + 2.e12)
172 print("100 %")
173 
174 plt.show()