I'm trying to validate an collision detector, and for that I'm starting from a pair of initial conditions for a satellite. After that I will use an RK method to integrate. How can I determine a pair of initial conditions for a second satellite that would lead to a collision? Let's say into an interval of one day (86400 sec).

  • 1
    $\begingroup$ For two satellites to collide their positions (x y z) need to match at the same time t. Unless you constrain the velocity some way, there are an infinite number of orbits that go through a position (x y z) at time t. Your problem is too general for an answer. $\endgroup$
    – cms
    Commented Aug 16, 2018 at 1:35
  • $\begingroup$ @cms OP asks for "a pair" not "the pair". You're right that there are an infinite number, but that doesn't stop the OP from just choosing a few. $\endgroup$
    – uhoh
    Commented Aug 16, 2018 at 3:47
  • 1
    $\begingroup$ Slightly related: Algorithmic methods or techniques to find conjunctions in high N Keplerian element ensembles? and also it's companion question. $\endgroup$
    – uhoh
    Commented Aug 16, 2018 at 4:08

1 Answer 1


How can I determine a pair of initial conditions for a second satellite that would lead to a collision?

Propagating an orbit is solving a set of differential equations by numerically integrating over time. The main effects are gravity and drag. As long as the spacecraft has a constant drag coefficient (it's not tumbling) these equations can be run backwards in time just as easily as they can be run forwards. As noted again below, when running backwards in time don't forget to change the sign of the drag force and if you are using any spherical harmonics other than the zonal harmonics (J2, J3, J4, ...), don't forget to rotate the Earth backwards as well!

The procedure would be to integrate the motion of the first spacecraft from your initial state $\mathbf{v_0}, \mathbf{x_0}$ for one day. Call the final position and velocity $\mathbf{x_1}, \mathbf{v_1}$.

Double check your integration by starting the spacecraft at $\mathbf{x_1}, -\mathbf{v_1}$ (the final position, but with the opposite velocity) and propagating for one day to make sure it arrives at $\mathbf{x_0}, -\mathbf{v_0}$.

Now start a second spacecraft at the collision point $\mathbf{x_1}$ but choose a velocity different than $-\mathbf{v_1}$, and propagate it for one day. Call the result $\mathbf{x_2}, -\mathbf{v_2}$ (note the negative sign in front of velocity).

When you choose a different velocity, don't forget to make sure it doesn't also intersect the Earth's atmosphere. Once you fix your initial orbit (see below) a good way to do this is to keep the radial velocity the same and just rotate the tangential velocity around the radius vector. That will result in a nearly identically shaped orbit but with a different inclination.

You are done!

If you start the first spacecraft at $\mathbf{v_0}, \mathbf{x_0}$ and the second spacecraft at $\mathbf{x_2}, \mathbf{v_2}$ and integrate both for one day, they should both arrive at $\mathbf{x_1}, \mathbf{v_1}$ within the precision limits of your RK method integrator and other numerical artifacts.

This should work for a central gravitational force with spherical harmonics (J2, etc.) but no Sun or Moon, plus a drag term dependent on velocity squared and on altitude using some atmospheric model. (See Atmospheric drag effect for some ideas.) However, if you use drag, you need to change the sign of your drag force when propagating backwards in time. Reminder that this only works for a constant drag coefficient and no pathological issues like atmospheric reentry.

Also, if you are using any spherical harmonics for gravity besides the zonal harmonics (J2, J3, J4, ...), you'll need to spin the Earth backwards as well.

You also need to pay attention to where each orbit goes. Your initial conditions seem to have a periapsis of about 3660 km, which is an altitude of about -2710 km!! which is below the surface of the Earth.

Below I've integrated using only the central force and no drag with a simple Python script. You can see that the spacecraft dips below the Earth's surface. You can add J2 and higher harmonics and drag after you get this going. (See this answer for some ideas about adding J2 and relativistic correction, but you'll need a more formal coordinate system to do it correctly. )

I got a specific energy of -5.42E+06 Joules/kg.

an orbit that hits the Earth

def deriv(X, t):

    x, v = X.reshape(2, -1)
    a    = -GMe * x * ((x**2).sum())**-1.5
    return np.hstack((v, a))

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

# https://space.stackexchange.com/questions/30125/how-to-determine-an-orbit-of-a-satellite-for-a-collision-detection

GMe = 3.986E+14
X0  = np.array([742453.224, -6345418.953, -3394102.502,
                -5142.381, -4487.449, -7264.009])

r0 = np.sqrt((X0[:3]**2).sum())
v0 = np.sqrt((X0[3:]**2).sum())

KE, PE = 0.5 * v0**2,  -GMe/r0
E      = KE + PE
print "E (Joules/kg): ", E

minutes = np.arange(2000.)
time    = 60. * minutes

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

xx, vv = answer.T.reshape(2, 3, -1)
rr     = np.sqrt((xx**2).sum(axis=0))

if True:
    plt.subplot(2, 1, 1)
    for thing in xx:
        plt.plot(minutes,  thing/1000.)
    plt.title('x, y, z (km) vs. time (minutes)')
    plt.subplot(2, 1, 2)
    plt.title('radius (km) vs. time (minutes)')
    plt.plot(minutes, rr/1000.)
    plt.plot(minutes, 6378.137 * np.ones_like(rr), '-k')
  • 4
    $\begingroup$ If only the zonal harmonics are used (J2, J3, J4, ...) the gravitational acceleration varies only with the latitude, so you don't need to rotate the Earth. While if also the tesseral and/or sectorial harmonics are used, you need to rotate the Earth. $\endgroup$
    – Cristiano
    Commented Aug 16, 2018 at 9:36
  • 1
    $\begingroup$ Probably my comment is not clear. If you write: "if you are using spherical harmonics above J2, don't forget to rotate the Earth backwards as well!" it seems to me that you are referring to the zonal harmonics only (J2, J3, J4, ..., Jn). If your statements "above J2" and "beyond J2" mean that you also include tesseral (J2,1, J3,1, J3,2, ...) and sectorial (J2,2, J3,3, ...) then I agree with your statements. $\endgroup$
    – Cristiano
    Commented Aug 16, 2018 at 14:55
  • 1
    $\begingroup$ @Cristiano I couldn't stop thinking about this so I did some further reading, and decided that there certainly could be reasons to cherry-pick just the zonal harmonics above J2, one example might be checking analytical expressions for perturbation against a numerical test, but I'm sure there are others. So I've taken your comment to heart and edited in two places. How does it look now? $\endgroup$
    – uhoh
    Commented Aug 18, 2018 at 4:51
  • 2
    $\begingroup$ Seems much better to me. It's common practice to include only the lower degree zonal harmonics (J2+J3 or J2 to J5) in simulations. When I don't need to use a full gravitational model (zonal, tesseral and sectorial harmonics), good results are also obtained with J2+J3 only (which is also much faster). $\endgroup$
    – Cristiano
    Commented Aug 18, 2018 at 10:20
  • 1
    $\begingroup$ Did you delete the answer? $\endgroup$
    – Cristiano
    Commented Sep 27, 2018 at 8:57

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.