tesla=sim.particles[-1]earth=sim.particles[3]r=np.linalg.norm(np.array(tesla.xyz)-np.array(earth.xyz))v=np.linalg.norm(np.array(tesla.vxyz)-np.array(earth.vxyz))energy=0.5*v*v-earth.m/rc3=2.*energy*887.40652# from units where G=1, length=1AU to km and sprint("c3 = %f (km^2/s^2)"%c3)

c3 = 11.943788 (km^2/s^2)

That seems about right! So let's look at the orbit. It starts at Earth's orbit, crosses that of Mars and then enters the asteroid belt.

In [45]:

rebound.OrbitPlot(sim,lim=1.8,slices=True,color=True);

And then integrate it forward in time. Here, we use the hybrid integrator MERCURIUS. You can experiment with other integrators which might be faster, but since this is an eccentric orbit, you might see many close encounters, so you either need a non-symplectic integrator such as IAS15 or a hybrid integrator such as MERCURIUS.

To check the sensitivity of the integrations, let us perturb the initial orbit by a small factor equal to the confidence interval posted by Bill Gray. Instead of just integrating one particle at a time, we here add 10 test particles. We also switch to the high precision IAS15 integrator to get the most reliable result.

sim.dt=sim.particles[1].P/60.# small fraction of Mercury's periodsim.integrator="ias15"N=1000times=np.linspace(0.,2.*np.pi*1e3,N)a_log=np.zeros((N,Ntesla))e_log=np.zeros((N,Ntesla))fori,tinenumerate(times):sim.integrate(t,exact_finish_time=0)forjinrange(Ntesla):orbit=sim.particles[9+j].calculate_orbit(primary=sim.particles[0])a_log[i][j]=orbit.ae_log[i][j]=orbit.e

When plotting the semi-major axis and eccentricity of all orbits, note that their kicks are correlated. This is because they are all due to close encounters with the Earth. This fast divergence means that we cannot predict the trajectory for more than a hundred years without knowing the precise initial conditions and all the non-gravitational effects that might be acting on a car in space.