Trajectory calculation and visualisation of spacecraft motion around the Earth

This example will show the calculation of the trajectory of a spacecraft in flight around the Earth, as well as the possibility of building a three-dimensional visualisation of its motion at different initial flight speeds.

Connection of necessary libraries:

Pkg.add(["LinearAlgebra"])

using Plots
using LinearAlgebra

Spacecraft trajectory calculation

Definition of constants and calculation step:

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

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

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)

Definition of initial conditions:

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

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

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

Calculating points in a spherical coordinate system:

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

The parameterisation of spherical coordinates is used to determine the points that characterise the Earth's surface.

The parameter phi represents the angle between the x-axis and the projection of the radius-vector of the point on the x-y plane. It varies from 0 to 2π, which corresponds to a complete revolution around the z-axis.

The parameter theta represents the angle between the radius-vector of the point and the positive direction of the z-axis.

For each combination of phi and theta values, the coordinates of the 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.

Definition of the function for transition from the spherical coordinate system to Cartesian coordinate system for displaying on the graph:

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)

Applying a function to change coordinate system:

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

Visualising the results obtained:

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

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)

Visualisation of the results of velocity calculations:

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 )

Visualisation of flight altitude calculation results:

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 )

Conclusion:

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