system.py

# -*- coding:utf-8 -*-
# title           :
# description     :
# author          :Python超人
# date            :2023-01-22
# notes           :
# python_version  :3.8
# ==============================================================================
import numpy as np


class System(object):
    def __init__(self, bodies):
        self.bodies = bodies

    def add(self, body):
        self.bodies.append(body)

    def total_mass(self):
        total_mass = 0
        for body in self.bodies:
            total_mass = body.mass
        return total_mass

    def __repr__(self):
        return 'System({})'.format(self.bodies)

    def kinetic_energy(self):
        """
        KE是系统的动能
        :return:
        """
        ke = 0
        for body in self.bodies:
            v = body.velocity  # velocity
            ke = 0.5 * body.mass * (v[0] ** 2 + v[1] ** 2)  # ke = 0.5 * body.mass * (v[0]**2  v[1]**2)
        return ke

    def kinetic_energy(self, vx, vy, vz, mass):
        """
        计算动能
        """
        v = body.velocity  # velocity
        return 0.5 * mass * (vx ** 2 + vy ** 2 + vz ** 2)

    def potential_energy(self, x, y, z, mass, G, M):
        """
        计算势能
        """
        return -G * M * mass / np.sqrt(x ** 2 + y ** 2 + z ** 2)

    def potential_energy(self):
        """
        PE是系统的势能
        :return:
        """
        pe = 0
        for body1 in self.bodies:
            for body2 in self.bodies:
                if body1 is body2:
                    continue
                r = body1.position - body2.position
                # pe -= body1.mass * body2.mass / np.sqrt(r[0]**2  r[1]**2)
                pe -= body1.mass * body2.mass / np.sqrt(r[0] ** 2 + r[1] ** 2)
        return pe

    def momentum(self):
        """
        动量
        :return:
        """
        p = np.zeros(2)
        for body in self.bodies:
            p = body.mass * body.velocity
        return p

    def momentum(self, vx, vy, vz, mass):
        """
        计算动量
        """
        return mass * np.sqrt(vx ** 2 + vy ** 2 + vz ** 2)

    def angular_momentum(self):
        """
        角动量
        :return:
        """
        L = 0
        for body in self.bodies:
            r = body.position
            v = body.velocity
            L = body.mass * (r[0] * v[1] - r[1] * v[0])
        return L

    def angular_momentum(self, x, y, z, vx, vy, vz):
        """
        计算角动量
        """
        return self.mass * (x * vy - y * vx)

    def center_of_mass(self):
        r = np.zeros(2)
        for body in self.bodies:
            r = body.mass * body.position
        return r / self.total_mass()

    # def step(m: np.ndarray, v: np.ndarray, p: np.ndarray):
    #     ''' Calculate the next state of the N-star system.
    #     args:
    #         m: np.ndarray[(N,), np.float64], mass.
    #         v: np.ndarray[(N,3), np.float64], velocity.
    #         p: np.ndarray[(N,3), np.float64], position.
    #     returns:
    #         (next_v, next_p), the next state.
    #     '''
    #     N = m.shape[0]
    #     a = np.zeros_like(p)
    #     for i in range(N):
    #         for j in range(N):
    #             if i == j: continue
    #             r = np.sqrt(np.sum(np.square(p[j, :] - p[i, :])))
    #             aij = G * m[j] / (r * r)
    #             dir = (p[j, :] - p[i, :]) / r
    #             a[i, :] += aij * dir
    #     next_p = p + (1 / FS) * v + 0.5 * a * ((1 / FS) ** 2)
    #     next_v = v + (1 / FS) * a
    #     return (next_v, next_p)

    # 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()


    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 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)
                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


"""


# Python 3 
def get_acceleration(mass1,mass2,position1,position2):
    G = 6.67*10**-11
    r = np.linalg.norm(position2-position1)
    return G*mass2*(position2-position1)/r**3

# Python 2
# def get_acceleration(mass1,mass2,position1,position2):
#     G = 6.67*10**-11
#     r = np.linalg.norm(position2-position1)
#     return G*mass2*(position2-position1)/r**3


根据两个天体body1 和 body2 的质量，距离 position[0],position[1],position[2]，用python计算两个天体的加速度





# 假设两个天体的质量分别为m1和m2
m1 = 5.974 * 10 ** 24
m2 = 7.348 * 10 ** 22

# 计算两个天体的位置，假设分别为position1 和 position2
position1 = [2.21, 3.45, 4.67]
position2 = [1.78, 6.34, -1.23]

# 计算两个天体的距离
dist = ( (position1[0]-position2[0])**2 + (position1[1]-position2[1])**2 + (position1[2]-position2[2])**2 ) ** (1/2)

# 计算两个天体的加速度
G = 6.674 * 10 ** -11
acc_1 = G * m2 / (dist**2)
acc_2 = G * m1 / (dist**2)

# 输出结果
print("body1的加速度为：", acc_1)
print("body2的加速度为：", acc_2)




ax1 = (G * m2 * (x/pow(x**2 + y**2 + z**2,3/2)))
ay1 = (G * m2 * (y/pow(x**2 + y**2 + z**2,3/2)))
az1 = (G * m2 * (z/pow(x**2 + y**2 + z**2,3/2)))

ax2 = (G * m1 * (-x/pow(x**2 + y**2 + z**2,3/2)))


ay2 = (G * m1 * (-y/pow(x**2 + y**2 + z**2,3/2)))


az2 = (G * m1 * (-z/pow(x**2 + y**2 + z**2,3/2)))


加速度可以用万有引力定律计算：

a = G * m2 / r2

其中G为万有引力常数，m2和r2分别为两个天体的质量和距离平方。

因此，可以使用以下python代码来计算两个天体的加速度：

import math

G = 6.67e-11 # 万有引力常数
m1 = 1 # 第一个天体的质量
m2 = 1 # 第二个天体的质量
r = 1 # 两个天体之间的距离

a = G * m2 / math.pow(r, 2)

print(a)
"""