1

I made an animation of an orbit that plots an orbit using an initial altitude and velocity.

Here is the code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
import matplotlib.animation as animation
plt.style.use('dark_background')
mu_earth = 398600
# Planet radius, km
R = 6371
fig = plt.figure()
ax1 = fig.add_subplot(111)
orbit_alt = 500
v_mag = 7.7 # km/s
# time step (delta-time)
dt = 10
satellite_size = 100
r_mag = R + orbit_alt
v_esc = np.sqrt(2*mu_earth*1e9/(R*1000 + orbit_alt*1000))
earth = plt.Circle((0, 0), 6371, color='blue')
ax1.add_artist(earth)

def diff_eq(t, y): 
    rx, ry, vx, vy = y # state vectors
    r = np.array([rx, ry])
    radial_vector = np.linalg.norm(r)
    ax, ay = -r * mu_earth / radial_vector ** 3
    return [vx, vy, ax, ay]

# Initial states
r0 = [r_mag, 0]
v0 = [0, v_mag]
specific_energy = (v_mag * 1000 * v_mag * 1000) / 2 - ((mu_earth * 1e9) / ((R + orbit_alt) * 1000))
semi_major = -mu_earth * 1e9 / (2 * specific_energy)
orbit_period = 2 * np.pi * np.sqrt(semi_major ** 3 / (mu_earth * 1e9))
if v_mag * 1000 < v_esc:
    tspan = 5 * orbit_period
else:
    tspan = 1500000
eccen = ((v_mag * 1000) ** 2 * (r_mag)*1000)/(mu_earth*1e9) - 1
perigee = (semi_major * (1 - np.abs(eccen)))/1000
apogee = (semi_major * (1 + np.abs(eccen)))/1000
n_steps = int(np.ceil(tspan/dt))
ts = np.zeros((n_steps, 1))
ys = np.zeros((n_steps, 4))
ys0 = r0 + v0
ts[0] = 0
ys[0] = ys0
solver = ode(diff_eq)
solver.set_integrator('lsoda')
solver.set_initial_value(ys0, 0)
rs = ys[:, :2]

def animate(i):
    solver.integrate(solver.t + dt)
    ts[i] = solver.t
    ys[i] = solver.y
    point = plt.Circle((ys[i][0], ys[i][1]), satellite_size, facecolor=(1, 1, 1))
    ax1.add_artist(point)
    satellite, = ax1.plot(rs[:, 0], rs[:, 1], color='r', alpha = 1)
    return satellite, point

ax1.set_xlim([-apogee, apogee])
ax1.set_ylim([-apogee, apogee])
ax1.set_xlabel('X (km)')
ax1.set_ylabel('Y (km)')
ax1.set_title("ORBIT SIMULATION")
plt.gca().set_aspect('equal', adjustable='box')
ani = animation.FuncAnimation(fig, animate, interval=0.01, blit=True)
plt.show()

And here is a screen shot of the output:

enter image description here

I want to omit the line that goes to the origin (0, 0). That line is a part of the satellite line object. What I believe is happening is that the origin is considered the first point, and the satellite is the second point, and matplotlib is connecting them. So how could that line leading to the middle be omitted?

2
  • Use zorder to specify which is on top? Commented Dec 22, 2020 at 18:14
  • zorder would cause the blue circle to cover the red line, but I belive there would likely still be some of the red line exposed. I believe the correct approach is to just omit that first point when plotting the satellite line object by plotting e.g. [1:] Commented Dec 22, 2020 at 18:51

1 Answer 1

2

The main problem is that your data is a full vector of zeros and your code plots the full dataset. Change from rs[:, 0] to rs[:, 0][:i] to only plot the points that have been simulated so far.

Second consideration is that FuncAnimation is usually used to update data already plotted. See small change below in animate(i).

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
import matplotlib.animation as animation
plt.style.use('dark_background')
mu_earth = 398600
# Planet radius, km
R = 6371
fig = plt.figure()
ax1 = fig.add_subplot(111)
orbit_alt = 500
v_mag = 7.7 # km/s
# time step (delta-time)
dt = 10
satellite_size = 100
r_mag = R + orbit_alt
v_esc = np.sqrt(2*mu_earth*1e9/(R*1000 + orbit_alt*1000))
earth = plt.Circle((0, 0), 6371, color='blue')
ax1.add_artist(earth)

def diff_eq(t, y):
    rx, ry, vx, vy = y # state vectors
    r = np.array([rx, ry])
    radial_vector = np.linalg.norm(r)
    ax, ay = -r * mu_earth / radial_vector ** 3
    return [vx, vy, ax, ay]

# Initial states
r0 = [r_mag, 0]
v0 = [0, v_mag]
specific_energy = (v_mag * 1000 * v_mag * 1000) / 2 - ((mu_earth * 1e9) / ((R + orbit_alt) * 1000))
semi_major = -mu_earth * 1e9 / (2 * specific_energy)
orbit_period = 2 * np.pi * np.sqrt(semi_major ** 3 / (mu_earth * 1e9))
if v_mag * 1000 < v_esc:
    tspan = 5 * orbit_period
else:
    tspan = 1500000
eccen = ((v_mag * 1000) ** 2 * (r_mag)*1000)/(mu_earth*1e9) - 1
perigee = (semi_major * (1 - np.abs(eccen)))/1000
apogee = (semi_major * (1 + np.abs(eccen)))/1000
n_steps = int(np.ceil(tspan/dt))
ts = np.zeros((n_steps, 1))
ys = np.zeros((n_steps, 4))
ys0 = r0 + v0
ts[0] = 0
ys[0] = ys0
solver = ode(diff_eq)
solver.set_integrator('lsoda')
solver.set_initial_value(ys0, 0)
rs = ys[:, :2]

orbit, = ax1.plot(rs[:, 0], rs[:, 1], color='r', alpha=1)
point = plt.Circle((ys[0][0], ys[0][1]), satellite_size, facecolor=(1, 1, 1))
ax1.add_artist(point)

def animate(i):
    solver.integrate(solver.t + dt)
    ts[i] = solver.t
    ys[i] = solver.y
    orbit.set_data(rs[:, 0][:i], rs[:, 1][:i])
    point.set_center((ys[i][0], ys[i][1]))
    return orbit, point

ax1.set_xlim([-apogee, apogee])
ax1.set_ylim([-apogee, apogee])
ax1.set_xlabel('X (km)')
ax1.set_ylabel('Y (km)')
ax1.set_title("ORBIT SIMULATION")
plt.gca().set_aspect('equal', adjustable='box')
ani = animation.FuncAnimation(fig, animate, interval=10, blit=True)
plt.show()
Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.