Faster routines for predicting occultations

[rodluger]

If we're interested in finding the next PPO among two planets, we need to integrate the system forward in time from the current state in order to plot its light curve. In other words, if I use system.next_occultation() and find that the next PPO is 100 days from now, I can't just run planetplanet in the vicinity of this occultation, since that won't capture the planet-planet interactions that will occur in the next 100 days. Instead, I would need to integrate the system forward in time, starting from a point where I have the complete orbital solution. Currently, there's no way to do it without computing the full light curve over those 100 days. Ideally, the user should be able to do the N-body integration up to the desired point quickly, and only evaluate the light curve where they actually need it.

This shouldn't be too hard to implement, but since we're still a year (at least) away from being able to do this, I'm deferring this to v2.