Engee 文档
Notebook

航天器环绕地球运动的轨迹计算和可视化

本示例将展示航天器环绕地球飞行的轨迹计算，以及在不同初始飞行速度下建立航天器运动三维可视化的可能性。

连接必要的库：

In [ ]: 
Pkg.add(["LinearAlgebra"])
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 库中的三维图将计算结果可视化的示例。在图上获得了轨迹，并得出了速度和与地球表面的距离。