Создание 3D модели Солнечной системы на Python с использованием Matplotlib

Одно из применений программирования (в том числе и на Python) заключается в том, чтобы помочь нам понять реальный мир с помощью симуляции. Этот метод используется в науке, финансах и многих других областях. В этой статье вы смоделируете Солнечную систему в 3D на Python, используя популярную библиотеку визуализации Matplotlib.

К концу этой статьи вы сможете создать свою собственную 3D модель Cолнечной системы на Python с таким количеством звёзд и планет, каким пожелаете. Вот пример простой Солнечной системы с одним Солнцем и двумя планетами:

@python_job_interview – практические задачи с Python собеседований.

Вы также сможете превратить её в 2D-проекцию, если захотите разнообразия. Ниже вы можете увидеть 2D макет Солнечной системы:

Краткое содержание статьи

Вот краткое содержание этой статьи:

• Краткое обсуждение гравитационного притяжения между двумя телами, которое вам нужно будет использовать для 3D моделирования Солнечной системы на Python.
• Краткое введение в векторы в 3D.
• Определение классов для солнечной системы и вращающихся внутри неё тел, таких как Солнце и планеты. Вы напишете эти классы в пошаговом подходе и протестируете их с помощью простой Солнечной системы.
• Добавление возможности показывать 2D-проекцию орбитальных тел вместе с 3D-моделированием. Эта 2D-проекция помогает визуализировать движение в 3D.
• Создание двойной звёздной системы.

В этой статье вы будете использовать объектно-ориентированное программирование и Matplotlib. Если вы хотите узнать побольше в этих сферах, вам будет полезно ознакомиться с ресурсами ниже:

• Object-Oriented Programming
• Basics of Data Visualisation in Python Using Matplotlib

Давайте начнём создание 3D модели Солнечной системы на Python с использованием библиотеки Matplotlib.

Давайте поговорим о Гравитации

Солнце, планеты и другие объекты в Солнечной системе – это тела, которые находятся в движении и которые притягивают друг друга из-за гравитационной силы, действующей между любыми двумя объектами.

Если два объекта имеют массы m1 и m2 и находятся на расстоянии r, то вы можете рассчитать гравитационную силу между ними, используя следующее уравнение:

Константа G – это гравитационная постоянная. Вы увидите, как вы сможете игнорировать эту константу в версии моделирования, которую вы напишете в этой статье. В этой статье вы будете использовать произвольные единицы измерения массы и расстояния, а не килограммы и метры.

Как только вы узнаете силу притяжения между двумя объектами, вы сможете рассчитать ускорение, которому подвергается каждый объект из-за силы притяжения, используя следующую формулу:

Используя это ускорение, вы можете регулировать скорость движущегося объекта. При изменении скорости изменятся как скорость, так и направление движения.

Представление точек и векторов в 3D

При создании 3D модели Солнечной системы на Python, вам нужно будет представить Солнечную систему как область пространства, используя три измерения. Следовательно, каждая точка в этом 3D-пространстве может быть представлена с помощью трёх чисел – x-, y- и z-координат. Например, если вы хотите поместить Солнце в центр Солнечной системы, вы можете представить положение солнца как (0, 0, 0).

Вам также нужно будет представлять векторы в 3D-пространстве. Вектор имеет как величину, так и направление. Вам понадобятся векторы для таких величин, как скорость, ускорение и сила.

Я не буду подробно обсуждать векторную алгебру в этой статье. Вместо этого я сообщу о любых результатах, которые вам понадобятся, по мере необходимости. При желании вы можете прочитать больше о векторах и векторной алгебре.

Чтобы упростить работу с векторами в коде, вы можете создать класс для работы с ними. Написание этого класса послужит быстрым обновлением знаний о классах и объектно-ориентированном программировании. Вы можете прочитать об объектно-ориентированном программировании на Python, если чувствуете, что вам нужно более подробное объяснение. Хотя вы также можете создать класс для работы с точками в 3D-пространстве, в этом нет необходимости, и я не буду создавать его в этой статье.

Создание Векторов

Если вы знакомы с векторами и объектно-ориентированным программированием, вы можете пропустить этот раздел и просто просмотреть код в конце, определяющий класс Vector.

Создайте новый файл с именем vectors.py, в котором вы будете определять класс Vector. Вы будете использовать этот скрипт для определения класса и его тестирования. Затем вы можете удалить тестовый код и оставить только определение класса:

# vectors.py
class Vector:
    def __init__(self, x=0, y=0, z=0):
        self.x = x
        self.y = y
        self.z = z
    def __repr__(self):
        return f"Vector({self.x}, {self.y}, {self.z})"
    def __str__(self):
        return f"{self.x}i + {self.y}j + {self.z}k"
# Testing Vector Class - TO BE DELETED
test = Vector(3, 5, 9)
print(test)
print(repr(test))
test = Vector(2, 2)
print(test)
print(repr(test))
test = Vector(y=5, z=3)
print(test)
print(repr(test))

Метод __init__() для класса Vector имеет три параметра, представляющих значение вдоль каждой оси. Каждый параметр имеет значение по умолчанию 0, представляющее начало координат для этой оси. Хотя мы предпочитаем не использовать имена из одной буквы в Python, x, y и z подходят, поскольку они представляют термины, обычно используемые в математике для декартовой системы координат.

Вы также определили два метода для представления объекта в виде строки:

• __repr__() возвращает вывод, предназначенный для программиста, показывающий имя класса. Выходные данные из repr() могут быть использованы для воссоздания объекта.
• __str__() возвращает непрограммную версию строкового представления объекта. В этом случае он возвращает изображение, которое обычно используется в математике для представления векторов, используя единичные векторы i, j и k.

Вы можете прочитать больше о различиях между двумя типами строковых представлений в разделе Snippets в конце главы 9 в книге The Python Coding Book.

Вот, что получится в поле вывода:

3i + 5j + 9k
Vector(3, 5, 9)
2i + 2j + 0k
Vector(2, 2, 0)
0i + 5j + 3k
Vector(0, 5, 3)

Делаем класс Vector индексируемым

В этом проекте 3D модели Солнечной системы было бы удобно, если бы класс Vector был индексируемым, чтобы вы могли использовать обозначение [] с индексом для извлечения одного из значений. С классом в его текущей форме, если вы добавите print(test[0]) в свой скрипт, вы получите TypeError. Вы можете исправить это, добавив другой магический метод в определение класса:

# vectors.py
class Vector:
    def __init__(self, x=0, y=0, z=0):
        self.x = x
        self.y = y
        self.z = z
    def __repr__(self):
        return f"Vector({self.x}, {self.y}, {self.z})"
    def __str__(self):
        return f"{self.x}i + {self.y}j + {self.z}k"
    def __getitem__(self, item):
        if item == 0:
            return self.x
        elif item == 1:
            return self.y
        elif item == 2:
            return self.z
        else:
            raise IndexError("There are only three elements in the vector")
# Testing Vector Class - TO BE DELETED
test = Vector(3, 5, 9)
print(test[0])

Определив __getitem__(), вы сделали класс Vector индексируемым. Первый элемент в векторе – это значение x, второй – значение y, а третий – значение z. Любой другой индекс вызовет ошибку. Вывод из блока тестового кода будет следующий:

3

test[0] возвращает первый элемент в векторе, значение для x.

Определение сложения и вычитания в классе Vector

Вы можете определить сложение и вычитание для объектов класса, определив методы __add__() и __sub__(). Эти методы позволят вам использовать символы + и – для выполнения операций. Без этих магических методов использование + и – вызывает ошибку TypeError.

Чтобы добавить или вычесть два вектора, вы можете добавить или вычесть каждый элемент векторов отдельно:

# vectors.py
class Vector:
    def __init__(self, x=0, y=0, z=0):
        self.x = x
        self.y = y
        self.z = z
    def __repr__(self):
        return f"Vector({self.x}, {self.y}, {self.z})"
    def __str__(self):
        return f"{self.x}i + {self.y}j + {self.z}k"
    def __getitem__(self, item):
        if item == 0:
            return self.x
        elif item == 1:
            return self.y
        elif item == 2:
            return self.z
        else:
            raise IndexError("There are only three elements in the vector")
    def __add__(self, other):
        return Vector(
            self.x + other.x,
            self.y + other.y,
            self.z + other.z,
        )
    def __sub__(self, other):
        return Vector(
            self.x - other.x,
            self.y - other.y,
            self.z - other.z,
        )
# Testing Vector Class - TO BE DELETED
test = Vector(3, 5, 9) + Vector(1, -3, 2)
print(test)
test = Vector(3, 5, 9) - Vector(1, -3, 2)
print(test)

И __add__() , и __sub__() возвращают другой объект Vector, каждый элемент которого равен сложению или вычитанию соответствующих элементов в двух исходных векторах. Результатом является следующее:

4i + 2j + 11k
2i + 8j + 7k

Вы можете сделать то же самое для умножения и деления, хотя эти операции требуют большей осторожности при работе с векторами.

Определение скалярного умножения, точечного произведения и скалярного деления в классе Vector

Вы не можете просто ссылаться на “умножение”, когда имеете дело с векторами, поскольку существуют различные типы ‘умножения’. В этом проекте вам понадобится только скалярное умножение. Скалярное умножение – это когда вектор умножается на скаляр (который имеет величину, но не имеет направления). Однако в этом подразделе вы также определите точечное произведение двух векторов. Вы могли бы использовать оператор * как для скалярного умножения, так и для точечного произведения. Следовательно, вы можете определить метод __mul__():

# vectors.py
class Vector:
    def __init__(self, x=0, y=0, z=0):
        self.x = x
        self.y = y
        self.z = z
    def __repr__(self):
        return f"Vector({self.x}, {self.y}, {self.z})"
    def __str__(self):
        return f"{self.x}i + {self.y}j + {self.z}k"
    def __getitem__(self, item):
        if item == 0:
            return self.x
        elif item == 1:
            return self.y
        elif item == 2:
            return self.z
        else:
            raise IndexError("There are only three elements in the vector")
    def __add__(self, other):
        return Vector(
            self.x + other.x,
            self.y + other.y,
            self.z + other.z,
        )
    def __sub__(self, other):
        return Vector(
            self.x - other.x,
            self.y - other.y,
            self.z - other.z,
        )
    def __mul__(self, other):
        if isinstance(other, Vector):  # Vector dot product
            return (
                self.x * other.x
                + self.y * other.y
                + self.z * other.z
            )
        elif isinstance(other, (int, float)):  # Scalar multiplication
            return Vector(
                self.x * other,
                self.y * other,
                self.z * other,
            )
        else:
            raise TypeError("operand must be Vector, int, or float")
# Testing Vector Class - TO BE DELETED
test = Vector(3, 5, 9) * Vector(1, -3, 2)
print(test)
test = Vector(3, 5, 9) * 3
print(test)

Результат использования оператора * будет зависеть от того, является ли второй операнд, следующий за символом *, скаляром или вектором. Если второй операнд, представленный параметром other, имеет тип Vector, вычисляется точечное произведение. Однако, если other имеет тип int или float, возвращаемый результат представляет собой новый вектор, масштабированный соответствующим образом.

Вывод из приведенного выше кода будет следующий:

6
9i + 15j + 27k

Если вы хотите скалярное умножение, скаляр должен идти после символа *. Если вы попытаетесь вместо этого запустить оператор 3*Vector(3, 5, 9), будет вызвано TypeError, поскольку класс Vector не является допустимым операндом для использования * с объектами типа int.

Два вектора не могут быть разделены. Однако вы можете разделить вектор на скаляр. Вы можете использовать оператор / с классом Vector, если вы определяете метод __truediv__():

# vectors.py
class Vector:
    def __init__(self, x=0, y=0, z=0):
        self.x = x
        self.y = y
        self.z = z
    def __repr__(self):
        return f"Vector({self.x}, {self.y}, {self.z})"
    def __str__(self):
        return f"{self.x}i + {self.y}j + {self.z}k"
    def __getitem__(self, item):
        if item == 0:
            return self.x
        elif item == 1:
            return self.y
        elif item == 2:
            return self.z
        else:
            raise IndexError("There are only three elements in the vector")
    def __add__(self, other):
        return Vector(
            self.x + other.x,
            self.y + other.y,
            self.z + other.z,
        )
    def __sub__(self, other):
        return Vector(
            self.x - other.x,
            self.y - other.y,
            self.z - other.z,
        )
    def __mul__(self, other):
        if isinstance(other, Vector):  # Vector dot product
            return (
                self.x * other.x
                + self.y * other.y
                + self.z * other.z
            )
        elif isinstance(other, (int, float)):  # Scalar multiplication
            return Vector(
                self.x * other,
                self.y * other,
                self.z * other,
            )
        else:
            raise TypeError("operand must be Vector, int, or float")
    def __truediv__(self, other):
        if isinstance(other, (int, float)):
            return Vector(
                self.x / other,
                self.y / other,
                self.z / other,
            )
        else:
            raise TypeError("operand must be int or float")
# Testing Vector Class - TO BE DELETED
test = Vector(3, 6, 9) / 3
print(test)

И на выходе получается:

1.0i + 2.0j + 3.0k

Нахождение величины вектора и нормализация вектора

Если у вас есть вектор (x, y, z)(x, y, z), вы можете найти его величину, используя выражение “корень из(x ^ 2 +y ^ 2 + z ^ 2)”. Вы также можете нормализовать вектор. Нормализация даеё вектор с тем же направлением, но с величиной 1. Вы можете вычислить нормализованный вектор, разделив каждый элемент вектора на его величину.

Вы можете определить два новых метода для завершения векторного класса:

# vectors.py
import math
class Vector:
    def __init__(self, x=0, y=0, z=0):
        self.x = x
        self.y = y
        self.z = z
    def __repr__(self):
        return f"Vector({self.x}, {self.y}, {self.z})"
    def __str__(self):
        return f"{self.x}i + {self.y}j + {self.z}k"
    def __getitem__(self, item):
        if item == 0:
            return self.x
        elif item == 1:
            return self.y
        elif item == 2:
            return self.z
        else:
            raise IndexError("There are only three elements in the vector")
    def __add__(self, other):
        return Vector(
            self.x + other.x,
            self.y + other.y,
            self.z + other.z,
        )
    def __sub__(self, other):
        return Vector(
            self.x - other.x,
            self.y - other.y,
            self.z - other.z,
        )
    def __mul__(self, other):
        if isinstance(other, Vector):  # Vector dot product
            return (
                self.x * other.x
                + self.y * other.y
                + self.z * other.z
            )
        elif isinstance(other, (int, float)):  # Scalar multiplication
            return Vector(
                self.x * other,
                self.y * other,
                self.z * other,
            )
        else:
            raise TypeError("operand must be Vector, int, or float")
    def __truediv__(self, other):
        if isinstance(other, (int, float)):
            return Vector(
                self.x / other,
                self.y / other,
                self.z / other,
            )
        else:
            raise TypeError("operand must be int or float")
    def get_magnitude(self):
        return math.sqrt(self.x ** 2 + self.y ** 2 + self.z ** 2)
    def normalize(self):
        magnitude = self.get_magnitude()
        return Vector(
            self.x / magnitude,
            self.y / magnitude,
            self.z / magnitude,
        )
# Testing Vector Class - TO BE DELETED
test = Vector(3, 6, 9)
print(test.get_magnitude())
print(test.normalize())
print(test.normalize().get_magnitude())

Тестовый код выдаёт следующий результат:

11.224972160321824

0.2672612419124244i + 0.5345224838248488j + 0.8017837257372732k

1.0

Третий вывод даёт величину нормализованного вектора, показывая, что его величина равна 1.

В зависимости от того, какую среду разработки или другие инструменты вы используете, вы можете получить предупреждение при разделении self.x, self.y и self.z, например, в __truediv__() и normalize(). Вам не нужно беспокоиться об этом, но если вы хотите это исправить, вы можете сделать изменить метод __init__() на любой другой из следующих:

def __init__(self, x=0.0, y=0.0, z=0.0):

или

def __init__(self, x:float=0, y:float=0, z:float=0):

Оба варианта позволяют вашей среде разработки знать, что аргументы должны иметь тип данных float.

Теперь вы можете удалить тестовый код в конце этого скрипта, чтобы в файле vectors.py осталось только это определение класса.

Создание 3D модели Солнечной системы на Python

Теперь вы можете начать работать над 3D моделью Солнечной системой на Python. Вы создадите два основных класса:

• SolarSystem: этот класс относится к Солнечной системе, отслеживает количество тел внутри неё и взаимодействия между ними.
  SolarSystemBody: этот класс посвящён каждому отдельному телу в Солнечной системе и движению тела.

Вы будете использовать Matplotlib для создания и визуализации Солнечной системы. Вы можете установить Matplotlib, используя следующие команды в терминале:

$ pip install matplotlib

или

$ python -m pip install matplotlib

Объект Axes3D в Matplotlib будет хостить солнечную систему. Если вы использовали Matplotlib для создания 2D-графики, вы бы эксплуатировали (сознательно или неосознанно) объект Axes. Axes3D – это 3D-эквивалент Axes, как следует из названия!

Пришло время приступить к написанию и тестированию этих классов. Вы можете создать два новых файла:

• solar_system_3d.py будет содержать определения классов.
• simple_solar_system.py будет содержать код для создания Солнечной системы. Вы будете использовать этот файл для тестирования классов по мере их написания, что приведёт к созданию простой Солнечной системы с одним Солнцем и двумя вращающимися планетами.

Далее вы начнёте работать над классом SolarSystem.

Настройка класса SolarSystem

Вы будете использовать произвольные единицы измерения на протяжении всего этого проекта. Это означает, что вместо того, чтобы использовать метры для расстояний и килограммы для масс, вы будете использовать величины без единиц измерения. Параметр size используется для определения размера куба, который будет содержать Солнечную систему:

# solar_system_3d.py
class SolarSystem:
    def __init__(self, size):
        self.size = size
        self.bodies = []
    def add_body(self, body):
        self.bodies.append(body)

Вы определяете класс SolarSystem с помощью метода __init__(), который включает в себя параметр size. Вы также определяете атрибут bodies. Этот атрибут представляет собой пустой список, который будет содержать все тела в пределах Солнечной системы, когда вы создадите их позже. Метод add_body() может быть использован для добавления орбитальных тел в Солнечную систему.

Следующим шагом будет внедрение Matplotlib. Вы можете создать фигуру и набор осей, используя функцию subplots() в matplotlib.pyplot:

# solar_system_3d.py
import matplotlib.pyplot as plt
class SolarSystem:
    def __init__(self, size):
        self.size = size
        self.bodies = []
        self.fig, self.ax = plt.subplots(
            1,
            1,
            subplot_kw={"projection": "3d"},
            figsize=(self.size / 50, self.size / 50),
        )
        self.fig.tight_layout()
    def add_body(self, body):
        self.bodies.append(body)

Вы вызываете plt.subplots(), который возвращает фигуру и набор осей. Возвращаемым значениям присваиваются атрибутам fig и ax. Вы вызываете plt.subplots() со следующими аргументами:

• Первые два аргумента равны 1 и 1 для создания единого набора осей на рисунке.
• Параметр subplot_kw имеет словарь в качестве аргумента, который устанавливает проекцию в 3D. Это означает, что созданные оси являются объектом Axes3D.
• figsize задает общий размер фигуры, содержащей объект Axes3D.

Вы также вызываете метод tight_layout(). Это метод класса Figure в Matplotlib. Он уменьшает поля по краям рисунка.

Вы можете попробовать код, приведенный ниже в консоли / REPL:

>>> import matplotlib.pyplot as plt
>>> from solar_system_3d import SolarSystem
>>> solar_system = SolarSystem(400)
>>> plt.show()  # if not using interactive mode

Это дает фигуру с пустым набором 3D-осей:

Дальше вы будете использовать параметр size, чтобы задать размер этого куба. Вы вернётесь к классу SolarSystem позже. На данный момент вы можете обратить своё внимание на определение класса SolarSystemBody.

Настройка класса SolarSystemBody

Вы можете начать создавать класс SolarSystemBody и его метод __init__(). Я усекаю код в определении класса SolarSystem в приведённом ниже коде для целей отображения. В этом и последующих блоках кода строки, содержащие # ..., указывают на код, который вы уже написали ранее и который не отображается:

# solar_system_3d.py
import matplotlib.pyplot as plt
from vectors import Vector
# class SolarSystem:
# ...  
class SolarSystemBody:
    def __init__(
        self,
        solar_system,
        mass,
        position=(0, 0, 0),
        velocity=(0, 0, 0),
    ):
        self.solar_system = solar_system
        self.mass = mass
        self.position = position
        self.velocity = Vector(*velocity)
        self.solar_system.add_body(self)

Параметрами в методе __init__() являются:

• solar_system позволяет вам связать тело с Солнечной системой. Аргумент должен быть типа SolarSystem.
• mass– это целое число или число с плавающей точкой, которое определяет массу тела. В этом проекте вы будете использовать произвольные единицы измерения, поэтому вам не нужно использовать “реальные” массы для звёзд и планет.
• position – это точка в трёхмерном пространстве, определяющая положение тела. Это кортеж, содержащий x-, y- и z-координаты точки. По умолчанию используется исходное значение.
• velocity определяет скорость тела. Поскольку скорость движущегося тела имеет величину и направление, она должна быть вектором. Хотя аргумент, необходимый при создании экземпляра SolarSystemBody , является кортежем, вы можете преобразовать кортеж в векторный объект, присвоив ему атрибут self.velocity.

Вы также вызываете метод add_body(), который вы определили ранее в классе SolarSystem, чтобы добавить это тело в Солнечную систему. Позже вы добавите немного больше к методу __init__().

Вы можете определить другой метод в SolarSystemBody для перемещения тела, используя его текущее положение и скорость:

# solar_system_3d.py
import matplotlib.pyplot as plt
from vectors import Vector
# class SolarSystem:
# ... 
class SolarSystemBody:
    def __init__(
        self,
        solar_system,
        mass,
        position=(0, 0, 0),
        velocity=(0, 0, 0),
    ):
        self.solar_system = solar_system
        self.mass = mass
        self.position = position
        self.velocity = Vector(*velocity)
        self.solar_system.add_body(self)
    def move(self):
        self.position = (
            self.position[0] + self.velocity[0],
            self.position[1] + self.velocity[1],
            self.position[2] + self.velocity[2],
        )

Метод move() переопределяет атрибут position на основе атрибута velocity. Мы уже обсуждали, как вы используете произвольные единицы измерения расстояния и массы. Вы также используете произвольные единицы измерения времени. Каждая “единица времени” будет представлять собой одну итерацию цикла, который вы будете использовать для запуска моделирования. Следовательно, move() сместит тело на величину, необходимую для одной итерации, что составляет одну единицу времени.

Рисование тел солнечной системы

Вы уже создали структуры Matplotlib, которая будет содержать Солнечную систему и все её тела. Теперь вы можете добавить метод draw() в SolarSystemBody для отображения тела на графике Matplotlib.

Прежде чем вы это сделаете, вам нужно будет определить ещё несколько атрибутов в SolarSystemBody, чтобы управлять цветом и размером маркеров, которые вы будете рисовать для представления тел:

# solar_system_3d.py
import math
import matplotlib.pyplot as plt
from vectors import Vector
# class SolarSystem:
# ... 
class SolarSystemBody:
    min_display_size = 10
    display_log_base = 1.3
    def __init__(
        self,
        solar_system,
        mass,
        position=(0, 0, 0),
        velocity=(0, 0, 0),
    ):
        self.solar_system = solar_system
        self.mass = mass
        self.position = position
        self.velocity = Vector(*velocity)
        self.display_size = max(
            math.log(self.mass, self.display_log_base),
            self.min_display_size,
        )
        self.colour = "black"
        self.solar_system.add_body(self)
    def move(self):
        self.position = (
            self.position[0] + self.velocity[0],
            self.position[1] + self.velocity[1],
            self.position[2] + self.velocity[2],
        )
    def draw(self):
        self.solar_system.ax.plot(
            *self.position,
            marker="o",
            markersize=self.display_size,
            color=self.colour
        )

Атрибуты класса min_display_size и display_log_base задают параметры для определения размера маркеров, которые вы будете отображать на 3D-графике. Вы устанавливаете минимальный размер, чтобы отображаемый маркер не был слишком маленьким даже для небольших тел. Вы будете использовать логарифмическую шкалу для преобразования массы в размер маркера, и вы задаете базу для этого логарифма в качестве другого атрибута класса.

Атрибут экземпляра display_size в методе  __init__() выбирает между рассчитанным размером маркера и установленным вами минимальным размером маркера. Чтобы определить размер отображения тела в этом проекте, вы используете его массу.

Вы также добавляете атрибут color в  __init__(), который на данный момент по умолчанию имеет значение black.

Чтобы протестировать новые дополнения, вы можете попробовать следующее в консоли / REPL:

>>> import matplotlib.pyplot as plt
>>> from solar_system_3d import SolarSystem, SolarSystemBody
>>> solar_system = SolarSystem(400)
>>> plt.show()  # if not using interactive mode
>>> body = SolarSystemBody(solar_system, 100, velocity=(1, 1, 1))
>>> body.draw()
>>> body.move()
>>> body.draw()

Первый вызов body.draw() рисует тело в начале координат, поскольку вы используете положение по умолчанию для объектов Солнечной системы. Вызов body.move() перемещает тело на величину, необходимую для одной “единицы времени”. Поскольку скорость тела равна (1, 1, 1), тело будет двигаться на одну единицу вдоль каждой из трёх осей. Второй вызов body.draw() рисует тело Солнечной системы во второй позиции. Обратите внимание, что при этом оси автоматически изменят масштаб. Вскоре вы позаботитесь об этом в основном коде.

Движение звёзд и планет

Вы можете вернуться к классу SolarSystem и связать Солнечную систему и её тела, добавив в класс два новых метода: update_all() и draw_all():

# solar_system_3d.py
import math
import matplotlib.pyplot as plt
from vectors import Vector
class SolarSystem:
    def __init__(self, size):
        self.size = size
        self.bodies = []
        self.fig, self.ax = plt.subplots(
            1,
            1,
            subplot_kw={"projection": "3d"},
            figsize=(self.size / 50, self.size / 50),
        )
        self.fig.tight_layout()
    def add_body(self, body):
        self.bodies.append(body)
    def update_all(self):
        for body in self.bodies:
            body.move()
            body.draw()
    def draw_all(self):
        self.ax.set_xlim((-self.size / 2, self.size / 2))
        self.ax.set_ylim((-self.size / 2, self.size / 2))
        self.ax.set_zlim((-self.size / 2, self.size / 2))
        plt.pause(0.001)
        self.ax.clear()
# class SolarSystemBody:
# ...

Метод update_all() проходит через каждое тело в Солнечной системе и перемещает и рисует каждое тело. Метод draw_all() устанавливает ограничения для трёх осей, используя размер Солнечной системы, и обновляет график с помощью функции pause(). Этот метод также очищает оси, подготавливая их к следующему графику.

Вы можете начать создавать простую Солнечную систему и протестировать код, который вы написали ранее, создав новый скрипт под названием simple_solar_system.py:

# simple_solar_system.py
from solar_system_3d import SolarSystem, SolarSystemBody
solar_system = SolarSystem(400)
body = SolarSystemBody(solar_system, 100, velocity=(1, 1, 1))
for _ in range(100):
    solar_system.update_all()
    solar_system.draw_all()

Когда вы запустите этот сценарий, вы увидите черное тело, удаляющееся от центра куба:

Вы можете изменить перспективу 3D-графика таким образом, чтобы вы просматривали 3D-оси непосредственно вдоль одной из осей. Вы можете сделать это, установив как азимут, так и высоту обзора равными 0 в SolarSystem.__init__():

# solar_system_3d.py
import math
import matplotlib.pyplot as plt
from vectors import Vector
class SolarSystem:
    def __init__(self, size):
        self.size = size
        self.bodies = []
        self.fig, self.ax = plt.subplots(
            1,
            1,
            subplot_kw={"projection": "3d"},
            figsize=(self.size / 50, self.size / 50),
        )
        self.fig.tight_layout()
        self.ax.view_init(0, 0)
    def add_body(self, body):
        self.bodies.append(body)
    def update_all(self):
        for body in self.bodies:
            body.move()
            body.draw()
    def draw_all(self):
        self.ax.set_xlim((-self.size / 2, self.size / 2))
        self.ax.set_ylim((-self.size / 2, self.size / 2))
        self.ax.set_zlim((-self.size / 2, self.size / 2))
        plt.pause(0.001)
        self.ax.clear()
# class SolarSystemBody:
# ...

Запущенный simple_solar_system.py теперь даёт следующий результат:

Ось x теперь перпендикулярна вашему экрану. Поскольку вы отображаете 3D-вид на 2D-дисплее, у вас всегда будет одно направление, перпендикулярное 2D-плоскости, которую вы используете для отображения графика. Это ограничение может затруднить распознавание того, когда объект движется вдоль этой оси. Вы можете увидеть это, изменив скорость тела в simple_solar_system.py (1, 0, 0) и снова запустив скрипт. Тело кажется неподвижным, так как оно движется только вдоль оси, выходящей из вашего экрана!

Помощь с 3D-перспективой

Вы можете улучшить 3D-визуализацию, изменив размер маркера в зависимости от его координаты x. Объекты, расположенные ближе к вам, кажутся больше, а объекты, расположенные дальше, кажутся меньше. Вы можете внести изменения в метод draw() в классе SolarSystemBody:

# solar_system_3d.py
# ...
class SolarSystemBody:
# ...
    def draw(self):
        self.solar_system.ax.plot(
            *self.position,
            marker="o",
            markersize=self.display_size + self.position[0] / 30,
            color=self.colour
        )

self.position[0] представляет положение тела вдоль оси x, которая является перпендикулярной экрану. Коэффициент 30, на который вы делите – это произвольный коэффициент, который вы можете использовать, чтобы контролировать, насколько сильным вы хотите, чтобы этот эффект был.

Позже в этом руководстве вы также добавите ещё одну функцию, которая поможет визуализировать 3D-движение звёзд и планет.

Добавление Эффектов Гравитации

У вас есть Солнечная система с телами, которые могут перемещаться внутри неё. Код до сих пор работает нормально, если у вас есть одно тело. Но это не очень интересная Солнечная система! Если у вас есть два или более тел, они будут взаимодействовать через их взаимное гравитационное притяжение.

В начале этой статьи я кратко рассмотрел физику, которая вам понадобится для работы с гравитационной силой между двумя объектами. Поскольку в этом проекте вы используете произвольные единицы измерения, вы можете игнорировать гравитационную постоянную G и просто вычислить силу, обусловленную гравитацией между двумя объектами, используя формулу:

Как только вы узнаете силу между двумя объектами, поскольку F = maF = ma, вы можете рассчитать ускорение, которому подвержен каждый объект, используя:

И как только вы узнаете ускорение, вы сможете изменить скорость объекта.

Вы можете добавить два новых метода, один в SolarSystemBody и другой в SolarSystem, чтобы вычислить силу и ускорение между любыми двумя телами и пройти через все тела в Солнечной системе, дабы выяснить взаимодействия между ними.

Отработка ускорения за счет силы тяжести

Первый из этих методов вычисляет силу тяжести между двумя телами, вычисляет ускорение каждого из тел и изменяет скорости двух тел. Вы можете разделить эти задачи на три метода, если хотите, но в этом примере я объединю их в один метод в SolarSystemBody:

# solar_system_3d.py
import math
import matplotlib.pyplot as plt
from vectors import Vector
# class SolarSystem:
# ...
class SolarSystemBody:
# ...
    def accelerate_due_to_gravity(self, other):
        distance = Vector(*other.position) - Vector(*self.position)
        distance_mag = distance.get_magnitude()
        force_mag = self.mass * other.mass / (distance_mag ** 2)
        force = distance.normalize() * force_mag
        reverse = 1
        for body in self, other:
            acceleration = force / body.mass
            body.velocity += acceleration * reverse
            reverse = -1

accelerate_due_to_gravity() вызывается для объекта типа SolarSystemBody и нуждается в другом теле SolarSystemBody в качестве аргумента. Параметры self и other представляют два тела, взаимодействующих друг с другом. Шаги в этом методе следующие:

• Положения двух тел используются для определения расстояния между ними. Вы представляете это как вектор, поскольку важны как его величина, так и направление. Вы извлекаете значения x, y и z из атрибута position с помощью оператора распаковки * и преобразуете их в объекты типа Vector, которые вы определили ранее. Поскольку вы определили метод sub() для класса Vector, вы можете вычесть один вектор из другого, чтобы получить расстояние между ними в виде другого вектора.
• Вы также вычисляете величину вектора расстояния, используя метод get_magnitude() класса Vector.
• Затем вы вычисляете величину силы между двумя телами, используя приведённое выше уравнение.
• Однако сила имеет не только величину, но и направление. Следовательно, вам нужно представить её в виде вектора. Направление силы совпадает с направлением вектора, соединяющего два объекта. Вы получаете вектор силы, сначала нормализуя вектор расстояния. Эта нормализация даёт единичный вектор с тем же направлением, что и вектор, соединяющий два тела, но с величиной 1. Затем вы умножаете единичный вектор на величину силы. В данном случае вы используете скалярное умножение вектора, которое вы определили, когда включили mul() в класс Vector.
• Для каждого из двух тел вы вычисляете ускорение, используя уравнение, показанное выше. force – это вектор. Следовательно, когда вы делите на body.mass, вы используете скалярное деление, которое вы определили, когда включили truediv() в класс Vector. acceleration  – это объект, возвращаемый вектором.truediv(), который также является векторным объектом.
• Наконец, вы увеличиваете скорость, используя ускорение. Этот метод вычисляет значения, относящиеся к одной единице времени, которая в данном моделировании представляет собой время, необходимое для одной итерации цикла, который будет управлять моделированием. Параметр reverse гарантирует, что ко второму телу прикладывается противоположное ускорение, поскольку два тела притягиваются друг к другу. Оператор * снова вызывает Vector.mul() и приводит к скалярному умножению.

Вычисление взаимодействий между всеми телами в Cолнечной системе

Теперь, когда вы можете определить взаимодействие между любыми двумя телами, вы можете определить взаимодействие между всеми телами, присутствующими в Солнечной системе. Для этого вы можете переключить сво` внимание обратно на класс SolarSystem:

# solar_system_3d.py
import math
import matplotlib.pyplot as plt
from vectors import Vector
class SolarSystem:
# ...
    def calculate_all_body_interactions(self):
        bodies_copy = self.bodies.copy()
        for idx, first in enumerate(bodies_copy):
            for second in bodies_copy[idx + 1:]:
                first.accelerate_due_to_gravity(second)
class SolarSystemBody:
# ...
    def accelerate_due_to_gravity(self, other):
        distance = Vector(*other.position) - Vector(*self.position)
        distance_mag = distance.get_magnitude()
        force_mag = self.mass * other.mass / (distance_mag ** 2)
        force = distance.normalize() * force_mag
        reverse = 1
        for body in self, other:
            acceleration = force / body.mass
            body.velocity += acceleration * reverse
            reverse = -1

Метод calculate_all_body_interactions() проходит через все тела в Солнечной системе. Каждое тело взаимодействует со всеми другими телами в Cолнечной системе:

• Вы используете копию self.bodies, чтобы учесть возможность того, что тела будут удалены из Cолнечной системы во время цикла. В версии, которую вы пишете в этой статье, вы не будете удалять какие-либо тела из Солнечной системы. Однако вам может понадобиться сделать это в будущем, если вы будете расширять этот проект дальше.
• Чтобы убедиться, что ваш код не вычисляет взаимодействия между одними и теми же двумя телами дважды, вы вычисляете только взаимодействия между телом и теми телами, которые следуют за ним в списке. Вот почему вы используете срез idx + 1: во втором цикле for.
• Последняя строка вызывает accelerate_due_to_gravity() для первого тела и включает второе тело в качестве аргумента метода.

Теперь вы готовы создать простую Солнечную систему и протестировать код, который вы написали ранее.

Создание простой Солнечной системы

В этом проекте вы сосредоточитесь на создании одного из двух типов тел: Солнца и планет. Вы можете создать два класса для этих тел. Новые классы наследуются от SolarSystemBody:

# solar_system_3d.py
import itertools
import math
import matplotlib.pyplot as plt
from vectors import Vector
# class SolarSystem:
# ...
# class SolarSystemBody:
# ...
class Sun(SolarSystemBody):
    def __init__(
        self,
        solar_system,
        mass=10_000,
        position=(0, 0, 0),
        velocity=(0, 0, 0),
    ):
        super(Sun, self).__init__(solar_system, mass, position, velocity)
        self.colour = "yellow"
class Planet(SolarSystemBody):
    colours = itertools.cycle([(1, 0, 0), (0, 1, 0), (0, 0, 1)])
    def __init__(
        self,
        solar_system,
        mass=10,
        position=(0, 0, 0),
        velocity=(0, 0, 0),
    ):
        super(Planet, self).__init__(solar_system, mass, position, velocity)
        self.colour = next(Planet.colours)

Класс Sun использует массу по умолчанию 10 000 единиц и устанавливает цвет на жёлтый. Вы используете строку ‘yellow‘, которая является допустимым цветом в Matplotlib.

В классе Planet вы создаёте объект itertools.cycle с тремя цветами. В данном случае три цвета – красный, зеленый и синий. Вы можете использовать любые цвета RGB, которые пожелаете. В этом классе вы определяете цвета, используя кортеж со значениями RGB вместо строки с названием цвета. Это также допустимый способ определения цветов в Matplotlib. Вы циклически меняете их с помощью функции next() каждый раз, когда создаёте новую планету.

Вы также устанавливаете массу по умолчанию равной 10 единицам.

Теперь вы можете создать Солнечную систему с одним Солнцем и двумя планетами в simple_solar_system.py:

# simple_solar_system.py
from solar_system_3d import SolarSystem, Sun, Planet
solar_system = SolarSystem(400)
sun = Sun(solar_system)
planets = (
    Planet(
        solar_system,
        position=(150, 50, 0),
        velocity=(0, 5, 5),
    ),
    Planet(
        solar_system,
        mass=20,
        position=(100, -50, 150),
        velocity=(5, 0, 0)
    )
)
while True:
    solar_system.calculate_all_body_interactions()
    solar_system.update_all()
    solar_system.draw_all()

В этом сценарии вы создаёте Солнце и две планеты. Вы присваиваете Солнце и планеты переменным с именами sun и planets, но это не является строго обязательным, поскольку после создания объектов Sun и Planet они добавляются в solar_system, и вам не нужно ссылаться на них напрямую.

Вы используете цикл while для запуска симуляции. Цикл выполняет три операции на каждой итерации. Когда вы запустите этот скрипт, вы получите следующую анимацию:

Вроде как работает. Вы можете видеть Солнце, закрепленное в центре этой Солнечной системы, и планеты, на которые влияет гравитационное притяжение Солнца. В дополнение к движению планет в плоскости, содержащей экран вашего компьютера (это оси y и z), вы также можете видеть, как планеты становятся всё больше или меньше, поскольку они также движутся по оси x, которая перпендикулярна вашему экрану.

Однако вы, возможно, заметили некоторое своеобразное поведение планет. Когда они должны находиться за Солнцем, планеты по-прежнему отображаются перед солнцем. Это не проблема с математикой — если вы отследите положение планет, вы увидите, что их координаты x показывают, что они на самом деле находятся за Солнцем, как и следовало ожидать.

Настройка показа тел, находящихся за другими телами

Данная проблема возникает из-за того, как Matplotlib рисует объекты на графике. Matplotlib отображает объекты слоями в том порядке, в котором вы их отображаете. Поскольку вы создали солнце перед планетами, объект Sun появляется первым в solar_system.bodies и рисуется как нижний слой. Вы можете проверить этот факт, создав Солнце после планет, и вы увидите, что в этом случае планеты всегда будут появляться позади Солнца.

Вы бы хотели, чтобы Matplotlib отображал тела Солнечной системы в правильном порядке, начиная с тех, которые находятся дальше всего. Чтобы достичь этого, вы можете сортировать список Solar System.bodies на основе значения координаты x каждый раз, когда вы хотите обновить 3D-график. Вот как вы можете сделать это с помощью метода update_all() в SolarSystem:

# solar_system_3d.py
import itertools
import math
import matplotlib.pyplot as plt
from vectors import Vector
class SolarSystem:
# ...
    def update_all(self):
        self.bodies.sort(key=lambda item: item.position[0])
        for body in self.bodies:
            body.move()
            body.draw()
# ...
# class SolarSystemBody:
# ...
# class Sun(SolarSystemBody):
# ...
# class Planet(SolarSystemBody):
# ...

Вы используете метод sort списка с параметром key, чтобы определить правило, которое вы хотели бы использовать для сортировки списка. Лямбда-функция устанавливает это правило. В этом случае вы используете значение position[0] каждого тела, которое представляет координату x. Следовательно, каждый раз, когда вы вызываете update_all() в цикле симуляции while, список тел переупорядочивается в зависимости от их положения вдоль оси x.

Результат выполнения simple_solar_system.py сценарий сейчас выглядит следующим образом:

Теперь вы можете визуализировать орбиты планет по мере их обращения вокруг Солнца. Изменяющийся размер показывает их положение по оси x, и когда планеты находятся за Солнцем, они скрыты из виду!

Наконец, вы также можете удалить оси и сетку, чтобы были видны только Солнце и планеты. Вы можете сделать это, добавив вызов метода Matplotlib axis() в SolarSystem.draw_all():

# solar_system_3d.py
import itertools
import math
import matplotlib.pyplot as plt
from vectors import Vector
class SolarSystem:
# ...
    def draw_all(self):
        self.ax.set_xlim((-self.size / 2, self.size / 2))
        self.ax.set_ylim((-self.size / 2, self.size / 2))
        self.ax.set_zlim((-self.size / 2, self.size / 2))
        self.ax.axis(False)
        plt.pause(0.001)
        self.ax.clear()
# ...
# class SolarSystemBody:
# ...
# class Sun(SolarSystemBody):
# ...
# class Planet(SolarSystemBody):
# ...

И симуляция теперь выглядит следующим образом:

3D моделирование Солнечной системы на Python с использованием Matplotlib теперь завершено. В следующем разделе вы добавите функцию, которая позволит вам просматривать 2D-проекцию плоскости xy в нижней части моделирования. Это может помочь в визуализации 3D-динамики тел в Солнечной системе.

Добавление 2D-проекции xy-плоскости

Чтобы помочь визуализировать движение тел при 3D-моделировании Солнечной системы на Python, вы можете добавить 2D-проекцию на “floor” анимации. Эта 2D-проекция покажет положение тел в плоскости xy. Чтобы достичь этого, вам нужно будет добавить другой график к тем же осям, по которым вы показываете анимацию, и показывать только изменения в координатах x и y. Вы можете привязать z-координату к нижней части графика, чтобы 2D-проекция отображалась в нижней части анимации.

Вы можете начать с добавления нового параметра в метод __init__() для класса SolarSystem:

# solar_system_3d.py
import itertools
import math
import matplotlib.pyplot as plt
from vectors import Vector
class SolarSystem:
    def __init__(self, size, projection_2d=False):
        self.size = size
        self.projection_2d = projection_2d
        self.bodies = []
        self.fig, self.ax = plt.subplots(
            1,
            1,
            subplot_kw={"projection": "3d"},
            figsize=(self.size / 50, self.size / 50),
        )
        self.ax.view_init(0, 0)
        self.fig.tight_layout()
# ...
# class SolarSystemBody:
# ...
# class Sun(SolarSystemBody):
# ...
# class Planet(SolarSystemBody):
# ...

Новый параметр projection_2d, значение которого по умолчанию равно False, позволит вам переключаться между двумя вариантами визуализации. Если projection_2d имеет значение False, анимация будет показывать только тела, движущиеся в 3D, без осей и сетки, как в последнем результате, который вы видели.

Давайте начнём вносить некоторые изменения, когда projection_2d имеет значение True:

# solar_system_3d.py
import itertools
import math
import matplotlib.pyplot as plt
from vectors import Vector
class SolarSystem:
    def __init__(self, size, projection_2d=False):
        self.size = size
        self.projection_2d = projection_2d
        self.bodies = []
        self.fig, self.ax = plt.subplots(
            1,
            1,
            subplot_kw={"projection": "3d"},
            figsize=(self.size / 50, self.size / 50),
        )
        self.fig.tight_layout()
        if self.projection_2d:
            self.ax.view_init(10, 0)
        else:
            self.ax.view_init(0, 0)
    def add_body(self, body):
        self.bodies.append(body)
    def update_all(self):
        self.bodies.sort(key=lambda item: item.position[0])
        for body in self.bodies:
            body.move()
            body.draw()
    def draw_all(self):
        self.ax.set_xlim((-self.size / 2, self.size / 2))
        self.ax.set_ylim((-self.size / 2, self.size / 2))
        self.ax.set_zlim((-self.size / 2, self.size / 2))
        if self.projection_2d:
            self.ax.xaxis.set_ticklabels([])
            self.ax.yaxis.set_ticklabels([])
            self.ax.zaxis.set_ticklabels([])
        else:
            self.ax.axis(False)
        plt.pause(0.001)
        self.ax.clear()
    def calculate_all_body_interactions(self):
        bodies_copy = self.bodies.copy()
        for idx, first in enumerate(bodies_copy):
            for second in bodies_copy[idx + 1:]:
                first.accelerate_due_to_gravity(second)
class SolarSystemBody:
    min_display_size = 10
    display_log_base = 1.3
    def __init__(
        self,
        solar_system,
        mass,
        position=(0, 0, 0),
        velocity=(0, 0, 0),
    ):
        self.solar_system = solar_system
        self.mass = mass
        self.position = position
        self.velocity = Vector(*velocity)
        self.display_size = max(
            math.log(self.mass, self.display_log_base),
            self.min_display_size,
        )
        self.colour = "black"
        self.solar_system.add_body(self)
    def move(self):
        self.position = (
            self.position[0] + self.velocity[0],
            self.position[1] + self.velocity[1],
            self.position[2] + self.velocity[2],
        )
    def draw(self):
        self.solar_system.ax.plot(
            *self.position,
            marker="o",
            markersize=self.display_size + self.position[0] / 30,
            color=self.colour
        )
        if self.solar_system.projection_2d:
            self.solar_system.ax.plot(
                self.position[0],
                self.position[1],
                -self.solar_system.size / 2,
                marker="o",
                markersize=self.display_size / 2,
                color=(.5, .5, .5),
            )
    def accelerate_due_to_gravity(self, other):
        distance = Vector(*other.position) - Vector(*self.position)
        distance_mag = distance.get_magnitude()
        force_mag = self.mass * other.mass / (distance_mag ** 2)
        force = distance.normalize() * force_mag
        reverse = 1
        for body in self, other:
            acceleration = force / body.mass
            body.velocity += acceleration * reverse
            reverse = -1
class Sun(SolarSystemBody):
    def __init__(
        self,
        solar_system,
        mass=10_000,
        position=(0, 0, 0),
        velocity=(0, 0, 0),
    ):
        super(Sun, self).__init__(solar_system, mass, position, velocity)
        self.colour = "yellow"
class Planet(SolarSystemBody):
    colours = itertools.cycle([(1, 0, 0), (0, 1, 0), (0, 0, 1)])
    def __init__(
        self,
        solar_system,
        mass=10,
        position=(0, 0, 0),
        velocity=(0, 0, 0),
    ):
        super(Planet, self).__init__(solar_system, mass, position, velocity)
        self.colour = next(Planet.colours)

Внесённые вами изменения заключаются в следующем:

• В SolarSystem.__init__() 3D-вид устанавливается в view_init(0, 0), когда 2D-проекция выключена, как и раньше. Однако высота изменяется на 10 градусов, когда включена опция 2D-проекции, чтобы нижняя часть была видна.
• В SolarSystem.draw_all() сетка и оси отключаются только тогда, когда отсутствует 2D-проекция. Когда включена 2D-проекция, они отображаются. Однако галочки заменяются пробелами, поскольку цифры на трёх осях произвольны и не нужны.
• В SolarSystemBody.draw() добавляется второй график, когда projection_2d имеет значение True. Первые два аргумента в plot() – это положения тел x и y. Однако вместо того, чтобы использовать z-позицию в качестве третьего аргумента, вы используете минимальное значение z, которое представляет “floor” куба, содержащего три оси. Затем вы наносите серый маркер вдвое меньшего размера, чем основные маркеры в анимации.

Вам также нужно будет внести небольшое изменение в simple_solar_system.py чтобы включить 2D-проекцию:

# simple_solar_system.py
from solar_system_3d import SolarSystem, Sun, Planet
solar_system = SolarSystem(400, projection_2d=True)
sun = Sun(solar_system)
planets = (
    Planet(
        solar_system,
        position=(150, 50, 0),
        velocity=(0, 5, 5),
    ),
    Planet(
        solar_system,
        mass=20,
        position=(100, -50, 150),
        velocity=(5, 0, 0)
    )
)
while True:
    solar_system.calculate_all_body_interactions()
    solar_system.update_all()
    solar_system.draw_all()

Симуляция теперь выглядит следующим образом:

2D-проекция плоскости xy облегчает отслеживание траекторий вращающихся тел.

Создание двойной Звездной системы

Мы закончим созданием ещё одной 3D моделью Солнечной системы на Python. Вы будете моделировать двойную Звездную систему, используя те же классы, которые вы уже определили. Создайте новый файл с именем binary_star_system.py. Сейчас мы создадим два Солнца и две планеты:

# binary_star_system.py
from solar_system_3d import SolarSystem, Sun, Planet
solar_system = SolarSystem(400)
suns = (
    Sun(solar_system, position=(40, 40, 40), velocity=(6, 0, 6)),
    Sun(solar_system, position=(-40, -40, 40), velocity=(-6, 0, -6)),
)
planets = (
    Planet(
        solar_system,
        10,
        position=(100, 100, 0),
        velocity=(0, 5.5, 5.5),
    ),
    Planet(
        solar_system,
        20,
        position=(0, 0, 0),
        velocity=(-11, 11, 0),
    ),
)
while True:
    solar_system.calculate_all_body_interactions()
    solar_system.update_all()
    solar_system.draw_all()

Моделирование этой двойной звездной системы выглядит следующим образом:

Или вы можете включить 2D-проекцию при создании объекта SolarSystem:

# binary_star_system.py
from solar_system_3d import SolarSystem, Sun, Planet
solar_system = SolarSystem(400, projection_2d=True)
# ...

Эта версия даёт следующий результат:

Эта двойная Звездная система нестабильна, и обе планеты вскоре выбрасываются из системы двумя Солнцами!

При желании вы можете расширить определения классов, чтобы обнаруживать столкновения между двумя телами и удалять планету, если она сталкивается с Солнцем. Более простая 2D-версия этого проекта, которая имитирует орбитальные планеты в 2D, включает в себя эту функцию. Вы можете посмотреть, как это было реализовано в том более простом проекте, если хотите добавить его в этот проект.

Окончательные версии кода, использованного в этой статье, также доступны в этом репозитории GitHub.

Заключение

Теперь вы можете смоделировать 3D Солнечную систему на Python, используя Matplotlib. В этой статье вы узнали, как размещать объекты в 3D-пространстве с помощью векторов и графических возможностей Matplotlib. Вы можете прочитать больше о том, как использовать Matplotlib, включая создание более сложных анимаций с использованием подмодуля animations в Matplotlib, в главе Основы визуализации данных в Python с использованием Matplotlib из книги The Python Coding Book.

Теперь ваша очередь пытаться создавать простые и более сложные Солнечные системы. Можете ли вы создать стабильную двойную Звёздную систему?

Я надеюсь, вам понравилось создавать 3D модель Солнечной системы на Python с помощью Matplotlib. Теперь вы готовы попробовать создать свои собственные симуляции реальных процессов.