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).
x=742453.224
y=-6345418.953
z=-3394102.502
vx=-5142.381
vy=-4487.449
vz=-7264.009
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.
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.figure()
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')
plt.show()
