system.py 8.8 KB
Newer Older
M
march3 已提交
1
# -*- coding:utf-8 -*-
M
march3 已提交
2 3
# title           :天体系统
# description     :天体系统,多个天体就是一个系统
M
march3 已提交
4
# author          :Python超人
M
march3 已提交
5 6
# date            :2023-02-11
# link            :https://gitcode.net/pythoncr/
M
march3 已提交
7 8 9
# python_version  :3.8
# ==============================================================================
import numpy as np
三月三net's avatar
三月三net 已提交
10
import math
M
march3 已提交
11
from common.consts import AU, G
M
march3 已提交
12
from bodies import Body, Sun, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto
M
march3 已提交
13
from common.func import calculate_distance
M
march3 已提交
14 15 16


class System(object):
M
march3 已提交
17 18 19 20
    """
    天体系统
    """

三月三net's avatar
三月三net 已提交
21
    def __init__(self, bodies, max_distance=200 * AU):
M
march3 已提交
22 23 24 25 26
        """

        :param bodies:
        :param max_distance:系统的最大范围,超出范围的天体就不显示了
        """
M
march3 已提交
27
        self.bodies = bodies
三月三net's avatar
三月三net 已提交
28
        # self.adjust_distance_and_velocity()
M
march3 已提交
29
        self.max_distance = max_distance
M
march3 已提交
30

三月三net's avatar
三月三net 已提交
31 32 33 34 35
    @staticmethod
    def calc_body_new_velocity_position(body, sun_mass=1.9891e30, G=6.674e-11):
        old_velocity = body.init_velocity
        old_position = body.init_position

三月三net's avatar
三月三net 已提交
36
        old_distance = np.linalg.norm(old_position - [0, 0, 0], axis=-1)
三月三net's avatar
三月三net 已提交
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
        new_distance = old_distance * body.distance_scale
        new_position = old_position * body.distance_scale

        new_velocity = System.get_new_velocity(old_velocity, old_distance, new_distance, body.mass)

        return new_velocity, new_position

    @staticmethod
    def get_new_velocity(old_velocity, old_distance, new_distance, mass, sun_mass=1.9891e30, G=6.674e-11):
        # 计算原速度的模长
        old_speed = np.linalg.norm(old_velocity * 1000)
        # 计算原动能和原势能
        old_kinetic_energy = 0.5 * mass * old_speed ** 2
        old_potential_energy = - G * mass * sun_mass / old_distance

        new_potential_energy = - G * mass * sun_mass / new_distance

        # 计算新动能
        new_kinetic_energy = old_kinetic_energy
        # 计算新速度的模长
        new_speed = math.sqrt(2 * (new_kinetic_energy - old_potential_energy) / mass)
        # 计算新速度向量
三月三net's avatar
三月三net 已提交
59
        new_velocity = old_velocity / old_speed * new_speed / 1000
三月三net's avatar
三月三net 已提交
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
        return new_velocity

    def get_new_velocity1(old_velocity, old_distance, new_distance, mass, sun_mass=1.9891e30, G=6.674e-11):
        # 计算原来的速度
        old_speed = math.sqrt(G * sun_mass / old_distance)
        # 计算新的速度
        new_speed = math.sqrt(G * sun_mass / new_distance)
        # 计算原来的动能
        old_kinetic_energy = 0.5 * mass * old_velocity ** 2
        # 计算新的动能
        new_kinetic_energy = old_kinetic_energy * new_speed ** 2 / old_speed ** 2
        # 计算新的速度
        new_velocity = math.sqrt(2 * new_kinetic_energy / mass)
        return new_velocity

M
march3 已提交
75 76 77 78 79
    def add(self, body):
        self.bodies.append(body)

    def total_mass(self):
        """
M
march3 已提交
80
        总质量
M
march3 已提交
81 82
        :return:
        """
M
march3 已提交
83
        total_mass = 0.0
M
march3 已提交
84
        for body in self.bodies:
M
march3 已提交
85 86
            total_mass += body.mass
        return total_mass
M
march3 已提交
87

M
march3 已提交
88 89
    def __repr__(self):
        return 'System({})'.format(self.bodies)
M
march3 已提交
90

M
march3 已提交
91
    def center_of_mass(self):
M
march3 已提交
92
        """
M
march3 已提交
93
        质心
M
march3 已提交
94 95
        :return:
        """
M
march3 已提交
96
        r = np.zeros(2)
M
march3 已提交
97
        for body in self.bodies:
M
march3 已提交
98 99
            r = body.mass * body.position
        return r / self.total_mass()
M
march3 已提交
100

M
march3 已提交
101
    def evolve(self, dt):
M
march3 已提交
102 103
        """

M
march3 已提交
104
        :param dt:
M
march3 已提交
105 106
        :return:
        """
M
march3 已提交
107
        self.calc_bodies_acceleration()
M
march3 已提交
108 109

        for body in self.bodies:
M
march3 已提交
110 111 112 113
            # acceleration 加速度
            body.velocity += body.acceleration * dt
            # body.position += 0.5 * body.acceleration * (dt ** 2)
            body.position += body.velocity * dt
M
march3 已提交
114

三月三net's avatar
三月三net 已提交
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
    def save_to_json(self, json_file_name, params=None):
        """

        :param json_file_name:
        :param params:
        :return:
        """
        import json
        import os
        # json_file = os.path.join("../data", json_file_name)
        filed_names = ["name", "mass", "init_position", "init_velocity",
                       "density", "color", "texture",
                       "size_scale", "distance_scale",  # "parent"
                       "rotation_speed", "ignore_mass", "is_fixed_star"]
        bodies = []
        for b in self.bodies:
            body = {}
            for filed_name in filed_names:
                filed_value = getattr(b, filed_name)
                if type(filed_value) is np.ndarray:
                    filed_value = filed_value.tolist()
                body[filed_name] = filed_value
            bodies.append(body)
        data = {"bodies": bodies}
        if params is not None:
            data["params"] = params
        json_str = json.dumps(data, indent=2, ensure_ascii=False, separators=(',', ': '))
        with open(json_file_name, "w", encoding='utf-8') as f:
            f.write(json_str)

M
march3 已提交
145
    def calc_bodies_acceleration(self):
M
march3 已提交
146
        """
M
march3 已提交
147 148
        计算加速度
        :return:
M
march3 已提交
149
        """
M
march3 已提交
150 151 152

        def valid_body(body):
            """
M
march3 已提交
153
            判断是否为有效的天体
M
march3 已提交
154 155 156
            :param body:
            :return:
            """
M
march3 已提交
157
            if not body.appeared:  # 不显示
M
march3 已提交
158
                return False
三月三net's avatar
三月三net 已提交
159 160 161 162 163
            # if self.max_distance > 0:
            #     # 超过了 max_distance 距离,则不显示,并消失
            #     if calculate_distance(body.position) > self.max_distance:
            #         body.appeared = False
            #         return False
M
march3 已提交
164 165 166

            return True

M
march3 已提交
167
        # self.bodies = list(filter(valid_body, self.bodies))
M
march3 已提交
168

M
march3 已提交
169
        for body1 in self.bodies:
三月三net's avatar
三月三net 已提交
170 171
            if body1.ignore_mass:
                continue
M
march3 已提交
172 173
            if not valid_body(body1):
                continue
M
march3 已提交
174
            acceleration = np.zeros(3)
M
march3 已提交
175
            for body2 in self.bodies:
三月三net's avatar
三月三net 已提交
176 177
                if body2.ignore_mass:
                    continue
M
march3 已提交
178
                if self.max_distance > 0:
M
march3 已提交
179
                    if calculate_distance(body1.position) > self.max_distance:  # 超过了max_distance距离,则消失
M
march3 已提交
180
                        body1.appeared = False
M
march3 已提交
181
                    if calculate_distance(body2.position) > self.max_distance:  # 超过了max_distance距离,则消失
M
march3 已提交
182 183
                        body2.appeared = False

三月三net's avatar
三月三net 已提交
184
                if not body1.appeared or not body2.appeared:
M
march3 已提交
185 186
                    continue

M
march3 已提交
187 188
                if body1 is body2:
                    continue
M
march3 已提交
189 190
                elif body1.ignore_gravity(body2) or body2.ignore_gravity(body1):
                    continue
M
march3 已提交
191 192

                r = body2.position - body1.position
M
march3 已提交
193 194 195 196
                # G = 6.67e-11 # 万有引力常数
                # m/s² = kg * m / m**3
                # km/s² = kg * m / m**3 / 1e9
                # acceleration = G * body2.mass * dx / (d ** 3)
M
march3 已提交
197
                acceleration += (G * body2.mass * r / np.linalg.norm(r) ** 3) / 1e9
M
march3 已提交
198

M
march3 已提交
199
            body1.acceleration = acceleration
M
march3 已提交
200 201


M
march3 已提交
202
if __name__ == '__main__':
三月三net's avatar
三月三net 已提交
203 204 205 206 207 208 209 210 211 212 213 214
    # body_sys = System([
    #     Sun(),  # 太阳
    #     Mercury(),  # 水星
    #     Venus(),  # 金星
    #     Earth(),  # 地球
    #     Mars(),  # 火星
    #     Jupiter(),  # 木星
    #     Saturn(),  # 土星
    #     Uranus(),  # 天王星
    #     Neptune(),  # 海王星
    #     Pluto()  # 冥王星(从太阳系的行星中排除)
    # ])
三月三net's avatar
三月三net 已提交
215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242
    # import math
    #
    # mass = 2e30
    # r = 2 * AU
    # # p = 14.9
    # p = 14.89
    # bodies = [
    #     Sun(name="太阳A红色", mass=mass,
    #         init_position=[0, r * math.sqrt(3), 0],  # 位置
    #         init_velocity=[-p, 0, 0],  # 速度(km/s)
    #         size_scale=5e1, texture="sun2.jpg", color=(255, 0, 0)),  # 太阳放大 100 倍
    #     Sun(name="太阳B绿色", mass=mass,
    #         init_position=[-r, 0, 0],
    #         init_velocity=[1 / 2 * p, -math.sqrt(3) / 2 * p, 0],
    #         size_scale=5e1, texture="sun2.jpg", color=(0, 255, 0)),  # 太阳放大 100 倍
    #     Sun(name="太阳C蓝色", mass=mass,
    #         init_position=[r, 0, 0],
    #         init_velocity=[1 / 2 * p, math.sqrt(3) / 2 * p, 0],
    #         size_scale=5e1, texture="sun2.jpg", color=(0, 0, 255)),  # 太阳放大 100 倍
    #     Earth(name="地球",
    #           # init_position=[0, -AU * -2, 5 * AU],
    #           init_position=[0, math.sqrt(3) * r / 6, 5 * AU],
    #           init_velocity=[0, 0, -10],
    #           size_scale=4e3, distance_scale=1),  # 地球放大 4000 倍,距离保持不变
    # ]
    # body_sys = System(bodies)
    # print(body_sys.save_to_json("../data/tri_bodies_sim_perfect_01.json"))
    earth = Earth(name="地球",
三月三net's avatar
三月三net 已提交
243 244 245 246
                  # init_position=[0, -AU * -2, 5 * AU],
                  init_position=[0, 1000000, 500000],
                  init_velocity=[0, 0, -10],
                  size_scale=4e3, distance_scale=1)
三月三net's avatar
三月三net 已提交
247 248 249
    new_velocity, new_position = System.calc_body_new_velocity_position(earth)

    print(new_velocity, new_position)
三月三net's avatar
三月三net 已提交
250
    print(earth.init_velocity, earth.init_position)