How to calculate the very special orbit of 2020 SO

Question

From a NASA JPL page:

Earth May Have Captured a 1960s-Era Rocket Booster

Source

(Sorry, I got an error when I tried to include an animated gif.)

2020 SO was captured by Earth's gravity on Nov. 8, 2020. It will escape in March 2021

This very special orbit could not be calculated as a two body problem.

At least we need to calculate the influence of the Sun and the Earth.

With the analysis of more than 170 detailed measurements of 2020 SO's position over the last three months, including observations made by the NASA-funded Catalina Sky Survey in Arizona and ESA's (European Space Agency's) Optical Ground Station in Tenerife, Spain, the impact of solar radiation pressure became evident and confirmed 2020 SO's low-density nature. The next step was to figure out where the suspected rocket booster could have come from.

So a three body problem is not enough, the solar radiation pressure should be calculated too.

Is there a simple solution using Python and Skyfield or another package?

Answer 1

I'll post this and then wait to see if the OP want's to take it from here and post an updated improved script, or would like improvements here.

---

I'm no expert at this but I will give it a try with an approximate calculation and then mention ways to make it better.

I'll start with positions of the Sun, Earth, Moon and 2020SO on 01-Nov-2020 from JPL horizons in barycentric coordinates and then propagate their motion.

I am lazy so for the Sun I will just let it drift relative to the solar system barycenter. We could add the four big planets if we wanted to get the motion of the Sun better.

I haven't included photon pressure on the 2020 SO. Since it varies as $1/r^2$ just like gravity, you could add that by reducing GMsun for just the line when the Sun's gravitational acceleration on 2020 SO is calculated, but I don't think it will make much difference over six months,

Force due to photon pressure will be something like this (pointed away from the Sun)

$$F \approx f \ A \ \frac{\text{1361 W/m}^2}{c^2} \frac{\text{1 AU}^2}{r^2}$$

where $f$ is some fudge factor of order unity involving reflectivity and diffuseness of scattering and $A$ is some average cross-sectional area.

Ways to improve:

1. Include relativistic corrections to the acceleration from answers to How to calculate the planets and moons beyond Newtons's gravitational force?
2. Add more planets to account for their effects on the motion of these bodies
3. Use actual ephemeris positions from Skyfield or Horizons or SPICE, either by calculating once and interpolating, or for Skyfield calling it repeatedly in the integration loop
4. Update to the new SciPy initial value problem ODE integrator scipy.integrate.solve_ivp which handles shaped arrays (so they don't have to be flattened like this) and possibly other goodies
5. Include some simple model for photon pressure

---

above: Earth in the center and the frame is rotating with the Earth's orbit around the Sun. below: Earth in the center and the frame doesn't rotate

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.integrate import odeint as ODEint

def rotz(vecs, th):
    x, y, z = vecs
    cth, sth = np.cos(th), np.sin(th)
    xrot = x * cth - y * sth
    yrot = y * cth + x * sth
    return np.vstack((xrot, yrot, z))

def deriv(X0, t):
    x, v = X0.reshape(2, -1)
    xsun, xearth, xmoon, xso = x.reshape(4, -1)
    
    asun = np.zeros(3) # just let the Sun drift (approximation)
    
    aearth = -GMs * (xearth-xsun) * (((xearth-xsun)**2).sum())**-1.5
    aearth += -GMm * (xearth-xmoon) * (((xearth-xmoon)**2).sum())**-1.5
    
    amoon = -GMs * (xmoon-xsun) * (((xmoon-xsun)**2).sum())**-1.5
    amoon += -GMe * (xmoon-xearth) * (((xmoon-xearth)**2).sum())**-1.5

    aso = -GMs * (xso-xsun) * (((xso-xsun)**2).sum())**-1.5
    aso += -GMe * (xso-xearth) * (((xso-xearth)**2).sum())**-1.5
    aso += -GMm * (xso-xmoon) * (((xso-xmoon)**2).sum())**-1.5

    acc = np.hstack((asun, aearth, amoon, aso))

    return np.hstack((v, acc))
    
# initial state vectors from JPL Horizons JD=2459154.500000000 (A.D. 2020-Nov-01 00:00:00.0000)

x0 = [-9.310337316714592E+05,  9.481788919374773E+05,  1.382575687749073E+04,
      1.148762587112595E+08,  9.386327656117186E+07,  9.411282390393317E+03,
      1.151732738477067E+08,  9.413906042360909E+07, -1.276638399431854E+04,
      1.164364627893544E+08,  9.471022406041500E+07, -2.793366484529004E+05]

v0 = [-1.242704706611964E-02, -9.038092487749897E-03,  3.763127777363607E-04,
      -1.913151322927344E+01,  2.310482932662013E+01, -1.501188465050873E-03,
      -1.979573309672475E+01,  2.381249918077239E+01,  6.837768022510993E-02,
      -1.954880640404515E+01,  2.295580439951134E+01,  6.207646654113397E-02]

X0 = np.array(x0 + v0) * 1000. # convert kilometers to meters

GMs = 1.32712440018E+20
GMe = 3.986004418E+14
GMm = 4.9048695E+12

# sample at 0.1 day intervals (integrator timesteps are variable and internal)

days = np.arange(0, 180.1, 0.1) # 181
times = days * 3600 * 24

# learn to use the newer scipy.integrate.solve_ivp later

answer, info = ODEint(deriv, X0, times, full_output=True)

positions, velocities = answer.T.reshape(2, 4, 3, -1)
sunpos, earthpos, moonpos, SOpos = positions

r_SO_moon = np.sqrt(((SOpos - moonpos)**2).sum(axis=0))
r_SO_earth = np.sqrt(((SOpos - earthpos)**2).sum(axis=0))
r_moon_earth = np.sqrt(((moonpos - earthpos)**2).sum(axis=0))

plt.figure()

#plt.plot(days, r_SO_moon/1000.)
#plt.plot(days, r_SO_earth/1000.)
#plt.plot(days, r_moon_earth/1000.)
## modified for legend
plt.plot(days, r_SO_moon/1000., label='2020 SO to Moon')  
plt.plot(days, r_SO_earth/1000., label='2020 SO to Earth')
plt.plot(days, r_moon_earth/1000., label='Moon to Earth')

plt.xlabel('days since 01-Nov-2020')
plt.ylabel('distance (km)')
plt.title('Separation between 2020 SO, Moon & Earth')

plt.legend(title='distance')   ## added for legend title

plt.show()

x, y, z = earthpos
theta = np.arctan2(y, x)

sunpos_rot, earthpos_rot, moonpos_rot, SOpos_rot = [rotz(thing, -theta)
                                                    for thing in positions]

plt.figure()
plt.plot([0], [0], 'ob')
x, y, z = moonpos_rot - earthpos_rot
plt.plot(x, y, '-g')
x, y, z = SOpos_rot - earthpos_rot
plt.plot(x, y, '-r')
plt.gca().set_aspect('equal')
plt.show()

plt.figure()
plt.plot([0], [0], 'ob')
x, y, z = moonpos - earthpos
plt.plot(x, y, '-g')
x, y, z = SOpos - earthpos
plt.plot(x, y, '-r')
plt.gca().set_aspect('equal')
plt.show()

Setup for Horizons to obtain initial state vectors in the script:

above: solar system barycenter. below: Type-2 state vectors, km, km/sec units, ecliptic and mean equinox, etc.

Answer 2

This is only a fragmentary answer. Using the data from this page

I calculated the gravitational acceleration to 2020 SO by the Sun and some planets:

body           mass    distance    distance  acceleration
units           kg        AU          km       m/s^2
perihelion  1.988e+30  6.189e-01  9.259e+07  1.548e-02
aphelion    1.988e+30  1.534e+00  2.295e+08  2.520e-03
mercury     3.300e+23  2.925e-01  4.375e+07  1.150e-08
venus       4.870e+24  1.907e-02  2.853e+06  3.994e-05
earth       5.970e+24  4.580e-02  6.852e+06  8.488e-06
mars        6.420e+23  1.236e-01  1.849e+07  1.253e-07
jupiter     1.898e+27  3.601e+00  5.387e+08  4.365e-07
saturn      5.680e+26  8.359e+00  1.250e+09  2.425e-08

At this point of the orbit the influence of the Sun but also of Venus and the Earth.

I modified the Python script by Uhoh to get this plot:

Here is the log scaled y axis :

d_SO_sun = SOpos - sunpos
d_earth_sun = earthpos - sunpos
a = -GMs * d_SO_sun * ((d_SO_sun**2).sum(axis=0))**-1.5
b = -GMs * d_earth_sun * ((d_earth_sun**2).sum(axis=0))**-1.5
diff = np.sqrt(((a-b)**2).sum(axis=0))

Now I would like to see the influence of some close planets and some very heavy planets.