Engee documentation
Notebook

Calculation of the trajectory and visualization of the spacecraft's motion around the Earth

This example will show the calculation of the spacecraft's trajectory during flight around the Earth, as well as the possibility of constructing a three-dimensional visualization of its movement at different initial flight speeds.

Connecting the necessary libraries:

In [ ]: 
Pkg.add(["LinearAlgebra"])
In [ ]: 
using Plots
using LinearAlgebra

Calculation of the spacecraft trajectory

Definition of constants and calculation step:

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

Defining a function that calculates the acceleration, velocity, and coordinates of a satellite:

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)

Determining the initial conditions:

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;  # сек

Calculation of spacecraft motion parameters:

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

Getting coordinates to define the Earth's surface on a graph:

Calculation of points in a spherical coordinate system:

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

Parameterization of spherical coordinates is used to determine the points characterizing the Earth's surface.

Parameter phi It represents the angle between the x-axis and the projection of the point's radius vector onto the x-y plane. It varies from 0 to 2π, which corresponds to a full rotation around the z axis.

Parameter theta represents the angle between the radius vector of the point and the positive direction of the z axis.

For each combination of values phi and theta The coordinates of a point are calculated using the transformation from spherical to Cartesian coordinates.

Each element of the resulting array x_earth-y_earth-z_earth corresponds to one point on the Earth's surface.

Defining a function for moving from a spherical coordinate system to a Cartesian coordinate system for displaying on a graph:

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)

Applying a function to change the coordinate system:

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

Visualization of the results obtained:

Visualization of the trajectory of a spacecraft moving around the Earth:

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]:

Visualization of the speed calculation results:

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]:

Visualization of flight altitude calculation results:

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]:

Conclusion:

In this example, a simplified calculation of the parameters of spacecraft movement around the Earth at different initial flight speeds was considered, and an example of using a three-dimensional graph of the Plots library to visualize the calculation results was demonstrated. Trajectories were obtained on the graph, as well as speeds and distances from the Earth's surface.