Документация Engee
Notebook

Расчёт траектории и визуализация движения космического аппарата вокруг Земли

Открыть пример в Engee

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

Подключение необходимых библиотек:

In [ ]: 
using Plots
using LinearAlgebra

Расчёт траектории космического аппарата

Определение констант и шага расчёта:

In [ ]: 
G = 6.67430e-11    # Гравитационная постоянная
M = 5.972e24       # Масса Земли, кг
R = 6371000.0      # Радиус Земли, м
dt = 50;           # Шаг по времени, с

Определение функции, рассчитывающей ускорение, скорость и координаты спутника:

In [ ]: 
function calculate_trajectory(initial_velocity, initial_position, time)
    # Инициализация массивов для хранения результатов
    x = [initial_position[1]]
    y = [initial_position[2]]
    z = [initial_position[3]]
    
    # Инициализация массива для хранения абсолютной скорости
    absolute_velocity = []

    # Инициализация массива для хранения расстояния от центра Земли
    absolute_altitude = []

    # Инициализация массивов для хранения скоростей и ускорений
    velocities = []
    accelerations = []

    # Вычисление новых скоростей по формуле скорости
    v = initial_velocity

    for t in 0:dt:time
        # Вычисление ускорений по закону тяготения
        r = norm([x[end], y[end], z[end]])
        a = -G * M / r^3 * [x[end], y[end], z[end]]
        
        # Вычисление новых скоростей по формуле скорости
        v += a * dt
        
        # Вычисление абсолютной скорости
        absolute_v = norm(v)
        
        # Вычисление расстояния от поверхности Земли
        absolute_a = norm(r) - R
        
        # Вычисление новых координат по формулам
        x_new = x[end] + v[1] * dt
        y_new = y[end] + v[2] * dt
        z_new = z[end] + v[3] * dt
        
        # Добавление новых координат в массивы
        push!(x, x_new)
        push!(y, y_new)
        push!(z, z_new)
        
        # Добавление новой абсолютной скорости в массив
        push!(absolute_velocity, absolute_v)
        
        # Добавление нового расстояния от центра Земли в массив
        push!(absolute_altitude, absolute_a)

        # Добавление новых скоростей и ускорений в массивы
        push!(velocities, v)
        push!(accelerations, a)
    end

    return x, y, z, velocities, accelerations, absolute_velocity, absolute_altitude
end
Out[0]:
calculate_trajectory (generic function with 1 method)

Определение начальных условий:

In [ ]: 
# Задание начальной скорости спутника
initial_velocity = [0, 7660.0, 0]  # м/с
initial_velocity1 = [0, 8660.0, 0]  # м/с
initial_velocity2 = [0, 9660.0, 0]  # м/с

# Задание начальных координат спутника
initial_position = [R + 408000.0, 0, 0]  # м

# Задание времени движения спутника
total_time = 11200.0;  # сек

Расчёт параметров движения космического аппарата:

In [ ]: 
times = 0:dt:total_time
x, y, z, velocities, accelerations, absolute_velocity, absolute_altitude = calculate_trajectory(initial_velocity, initial_position, total_time)
x1, y1, z1, velocities1, accelerations1, absolute_velocity1, absolute_altitude1 = calculate_trajectory(initial_velocity1, initial_position, total_time)
x2, y2, z2, velocities2, accelerations2, absolute_velocity2, absolute_altitude2 = calculate_trajectory(initial_velocity2, initial_position, total_time);

Получение координат для определения поверхности Земли на графике:

Расчёт точек в сферической системе координат:

In [ ]: 
phi = range(0, stop=2*pi, length=100)
theta = range(0, stop=pi, length=50)

x_earth = zeros(length(phi), length(theta))
y_earth = zeros(length(phi), length(theta))
z_earth = zeros(length(phi), length(theta));

Для определения точек, характеризующих поверхность Земли, используется параметризация сферических координат.

Параметр phi представляет собой угол между осью x и проекцией радиус-вектора точки на плоскость x-y. Он изменяется от 0 до 2π, что соответствует полному обороту вокруг оси z.

Параметр theta представляет собой угол между радиус-вектором точки и положительным направлением оси z.

Для каждой комбинации значений phi и theta вычисляются координаты точки, используя преобразование из сферических в декартовы координаты.

Каждый элемент полученного массива x_earth-y_earth-z_earth соответствует одной точке поверхности Земли.

Определение функции для перехода из сферической системы координат в декартову для отображения на графике:

In [ ]: 
function spherical_to_cartesian(phi::Float64, theta::Float64)
    x = R * sin(theta) * cos(phi)
    y = R * sin(theta) * sin(phi)
    z = R * cos(theta)
    return (x, y, z)
end
Out[0]:
spherical_to_cartesian (generic function with 1 method)

Применение функции для смены системы координат:

In [ ]: 
for i in 1:length(phi)
    for j in 1:length(theta)
        result = spherical_to_cartesian(phi[i], theta[j])
        x_earth[i,j] = result[1]
        y_earth[i,j] = result[2]
        z_earth[i,j] = result[3]
    end
end

Визуализация полученных результатов:

Визуализация траектории космического аппарата двигающегося вокруг Земли:

In [ ]: 
plotlyjs() # Подключение бэкенда - метода отображения графики

# Визуализация траектории спутника
plot(x, y, z, zcolor=z, legend=false, labels="v0 = 7660 м/с", linewidth=5,
    color=:red, title="Траектории движения спутника<br>при разной начальной скорости",
    titlefont=font(11)) # Начальная скорость КА 7660 м/с
plot!(x1, y1, z1, zcolor=z, legend=false, labels="v0 = 8660 м/с", linewidth=5,
    color=:green) # Начальная скорость КА 8660 м/с
plot!(x2, y2, z2, zcolor=z, legend=false, labels="v0 = 9660 м/с", linewidth=5,
    color=:purple) # Начальная скорость КА 9660 м/с

# Визуализация поверхности Земли
surface!(x_earth, y_earth, z_earth, color=:blue, labels="Земля", alpha=0.9)
Out[0]:

Визуализация результатов расчёта скоростей:

In [ ]: 
plot(times/60, absolute_velocity/1000, linewidth=3, color=:red, labels="v0 = 7660 м/с",
    xlabel="Время, мин", ylabel="Скорость, км/с", title="Изменение скорости спутника во времени")
plot!(times/60, absolute_velocity1/1000, linewidth=3, color=:green, labels="v0 = 8660 м/с")
plot!(times/60, absolute_velocity2/1000, linewidth=3, color=:purple, labels="v0 = 9660 м/с",
    legend=:bottomleft )
Out[0]:

Визуализация результатов расчёта высоты полёта:

In [ ]: 
plot(times/60, absolute_altitude/1000000, linewidth=3, color=:red, labels="v0 = 7660 м/с",
    xlabel="Время, мин", ylabel="Высота, тыс. км", title="Изменение высоты полёта во времени")
plot!(times/60, absolute_altitude1/1000000, linewidth=3, color=:green, labels="v0 = 8660 м/с")
plot!(times/60, absolute_altitude2/1000000, linewidth=3, color=:purple, labels="v0 = 9660 м/с",
    legend=:left )
Out[0]:

Вывод:

В данном примере был рассмотрен упрощённый расчёт параметров движения космического аппарата вокруг Земли при разных начальных скоростях полёта, продемонстрирован пример использования трёхмерного графика библиотеки Plots для визуализации результатов расчёта. На графике были получены траектории, а также были выведены скорости и расстояния от поверхности Земли.