# How to calculate the very special orbit of 2020 SO

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?

• Don't forget that (according to the animation) the capture is triggered by the Moon as is its ejection. First it passes in front of the Moon and slows down, then it follows the Moon and speeds up. So it's a 4-body problem plus solar wind. Nov 29, 2020 at 10:25
• @asdfex You are right, 2020 SO gets very close to the Moon two times.
– Uwe
Nov 29, 2020 at 10:31
• – uhoh
Nov 29, 2020 at 14:09
• Here are the distances to the 8 planets, I would like to calculate and compare the forces, may be it is more than a 4-body problem.
– Uwe
Nov 29, 2020 at 14:40
• @uhoh the short Wikipedia article is not written in German but in Netherlands language. I only understand some very few words. But there is a German version too. I am sorry but I don't know the special English words for mechanical engineering
– Uwe
Jan 7, 2021 at 2:17

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(, , '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(, , '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. • Thanks a lot for the script. I added a legend to the last plot. There was no easy example for a legend at mathplotlib. So I tried to add label='2020 SO to Moon'. Only one additional line was neccessary for the title of the legend. I add the modified plot and the modified script lines.
– Uwe
Nov 30, 2020 at 8:38
• @Uwe sweet! :-)
– uhoh
Nov 30, 2020 at 9:51
• @uhoh Could you please tell me how you got those initial state vectors (ephemerides?!). I tried searching for them on JPL HORIZONS but couldn't get these values. Maybe I am searching it wrong. Thanks! Dec 2, 2020 at 19:13
• @uhoh Ok got it for all the celestial objects. That second image was really helpful where I needed to change the table settings to Type 2 and output units to km & km/s. Thank you so much! Also, I always look forward to your Python scripts - they make the answers much more interesting :) Dec 3, 2020 at 15:18
• @OrangeDurito Python rocks! :-)
– uhoh
Dec 3, 2020 at 23:36 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.

• This is interesting, here's some comments: 1) for the plot it would be interesting to see the accelerations on a log scale. Small effects may not seem so important but they would be easier to see. It did pass pretty close to Earth so even small effects might have had some impact. 2) In terms of effects near time of cis-lunar interaction by large-but-distant bodies, remember that they accelerate the Earth, Moon and SO all about the same amount, so they have only a tiny differential (tidal) effect on SO relative to the Earth-Moon barycenter.
– uhoh
Nov 30, 2020 at 22:31
• So if you for example plotted the difference between Jupiter's acceleration on the Earth-Moon barycenter and on SO, it might actually be very tiny and Venus' effect might be much larger if it happened to be close by. Tidal effects I think scale as $1/r^3$ so Venus' differential effect could be 20x larger than Jupiter's even though Jupiter's mass is 400x larger. (if at inferior conjunction) Since Venus's perturbations are more "pulse-like" rather than steady, it can be a real trouble-maker for short duration phenomenon like a temporary mini-moon.
– uhoh
Nov 30, 2020 at 22:38
• I added the log scale, looks much better. We see a lot more of the Moons influence.
– Uwe
Nov 30, 2020 at 22:52
• It would be interesting to add one more to the plot; the difference between Sun's acceleration on SO and on Earth (since Earth is pretty close to the Earth-Moon barycenter it serves as a good proxy). I think that will drop way down. This is why we don't include the Sun's gravitational effect on satellites in LEO (short term) unless they are far from Earth. How do SDP4's “Deep space” corrections to SGP4 account for the Sun's and Moon's gravity?
– uhoh
Nov 30, 2020 at 22:58
• For the difference between Sun's acceleration on SO and on Earth, I calculate the distance Sun to Earth and Sun to SO, get both accelerations from those distances and then difference of both accelerations? But should we calculate a three dimensional difference of those forces? The angle between the forces is not constant.
– Uwe
Nov 30, 2020 at 23:24