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

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

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

calculate_trajectory (generic function with 1 method)

# Задание начальной скорости спутника
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;  # сек

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

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

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

spherical_to_cartesian (generic function with 1 method)

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

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)

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 )

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 )