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    # The gravitational constant
M = 5.972e24       # The mass of the Earth, kg
R = 6371000.0      # The radius of the Earth, m
dt = 50;           # Time step, with

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

In [ ]: 
function calculate_trajectory(initial_velocity, initial_position, time)
    # Initializing arrays for storing results
    x = [initial_position[1]]
    y = [initial_position[2]]
    z = [initial_position[3]]
    
    # Initializing an array to store the absolute speed
    absolute_velocity = []

    # Initializing an array to store the distance from the center of the Earth
    absolute_altitude = []

    # Initializing arrays for storing speeds and accelerations
    velocities = []
    accelerations = []

    # Calculating new speeds using the speed formula
    v = initial_velocity

    for t in 0:dt:time
        # Calculation of accelerations according to the law of gravity
        r = norm([x[end], y[end], z[end]])
        a = -G * M / r^3 * [x[end], y[end], z[end]]
        
        # Calculating new speeds using the speed formula
        v += a * dt
        
        # Calculating the absolute speed
        absolute_v = norm(v)
        
        # Calculating the distance from the Earth's surface
        absolute_a = norm(r) - R
        
        # Calculating new coordinates using formulas
        x_new = x[end] + v[1] * dt
        y_new = y[end] + v[2] * dt
        z_new = z[end] + v[3] * dt
        
        # Adding new coordinates to arrays
        push!(x, x_new)
        push!(y, y_new)
        push!(z, z_new)
        
        # Adding a new absolute speed to the array
        push!(absolute_velocity, absolute_v)
        
        # Adding a new distance from the center of the Earth to the array
        push!(absolute_altitude, absolute_a)

        # Adding new speeds and accelerations to arrays
        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 [ ]: 
# Setting the initial velocity of the satellite
initial_velocity = [0, 7660.0, 0]  # m/s
initial_velocity1 = [0, 8660.0, 0]  # m/s
initial_velocity2 = [0, 9660.0, 0]  # m/s

# Setting the initial coordinates of the satellite
initial_position = [R + 408000.0, 0, 0]  # m

# Setting the satellite's travel time
total_time = 11200.0;  # sec

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 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.

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() # Enabling a backend graphics display method

# Visualization of the satellite trajectory
plot(x, y, z, zcolor=z, legend=false, labels="v0 = 7660 m/s", linewidth=5,
    color=:red, title="Satellite trajectories<br>at different initial speeds",
    titlefont=font(11)) # The initial velocity of the spacecraft is 7660 m/s
plot!(x1, y1, z1, zcolor=z, legend=false, labels="v0 = 8660 m/s", linewidth=5,
    color=:green) # The initial velocity of the spacecraft is 8660 m/s
plot!(x2, y2, z2, zcolor=z, legend=false, labels="v0 = 9660 m/s", linewidth=5,
    color=:purple) # The initial velocity of the spacecraft is 9660 m/s

# Visualization of the Earth's surface
surface!(x_earth, y_earth, z_earth, color=:blue, labels="Earth", alpha=0.9)
Out[0]:

Visualization of speed calculation results:

In [ ]: 
plot(times/60, absolute_velocity/1000, linewidth=3, color=:red, labels="v0 = 7660 m/s",
    xlabel="Time, min", ylabel="Speed, km/s", title="Changing the satellite's speed over time")
plot!(times/60, absolute_velocity1/1000, linewidth=3, color=:green, labels="v0 = 8660 m/s")
plot!(times/60, absolute_velocity2/1000, linewidth=3, color=:purple, labels="v0 = 9660 m/s",
    legend=:bottomleft )
Out[0]:

Visualization of flight altitude calculation results:

In [ ]: 
plot(times/60, absolute_altitude/1000000, linewidth=3, color=:red, labels="v0 = 7660 m/s",
    xlabel="Time, min", ylabel="Altitude, thousand km", title="Change in flight altitude over time")
plot!(times/60, absolute_altitude1/1000000, linewidth=3, color=:green, labels="v0 = 8660 m/s")
plot!(times/60, absolute_altitude2/1000000, linewidth=3, color=:purple, labels="v0 = 9660 m/s",
    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.