I actually ran into this about 6 months ago when I was fooling around with SA/Aurora-style (i.e. max velocity, instant acceleration) intercepts. It (the orbits of orbits of orbits) is a REALLY nasty problem. A couple of suggestions:
1) Hierarchically throw away detail. By this I mean you want to predict the intercept to the outermost orbit (i.e. the 2nd star) 1st. That will give you an intercept time. If the intercept time is long compared to the next period, stop there and pretend you're aiming for the central body. If it's not recursively go to the next body. As an example in the Earth/Luna system, first you'd calculate an intercept to Sol. If the time was more than a year, you'd stop and simply aim at Sol. If the time was less than a year but more than a month, you'd calculate the intercept to Earth and ignore Luna. If less than a month, you'd simplify the orbit of Earth (see below) and intercept Luna. Note that this assumes that an orbiter has a significantly shorter period than an orbitee - if that's not the case you've got a big problem 
2) Use a "constant acceleration" approximation (i.e. a Taylor expansion) to approximate orbits whose period is much longer than the intercept time. In the example below (assuming Sol is not moving), when aiming at Luna you would not model Earth's motion as a circle. Instead you'd calculate Earth's velocity and acceleration. (I would calculate them relative to Sol, then add in Sol's velocity and accel to get a total. This is because, relative to Sol, they'll always have the same magnitude and simply have different directions.) You'd then approximate Luna's motion as a circle orbiting a object moving along a parabola. The reason for doing this is that a) it makes the math easier and b) it's recursive - no matter how many nestings you've got, you're still always solving the same problem. This is simply a more sophisticated version of what you're doing right now, i.e. simply using the 0'th derivative (position) to approximate the correct course. When you get really close to the target (say an intercept time of less than Period/4), then you can approximate the targets motion too, at which point you're calculating an intercept to something with a simple parabolic path (i.e. intercepting a constant acceleration body).
3) Looking at the above, if the "circle around a parabola" problem is too hard, then you could just ignore the orbiting body until intercept time was less than Period/4, then pop over directly to approximating its motion with constant acceleration. So at that point the two changes to your existing algorithm would be a) ignore nested orbits until the intercept time drops to some fraction of a period and b) approximate all orbital motion as constant velocity (no accel) or constant accel, which allows you to use the same intercept code you'd use for ships.
4) Another thing I just thought of: when calculating the intercept time to see if you should igore a sub-orbit, I really meant "intercept time to the closest point on the orbital circle". In other words, you should use the motion of the central body, but subtract off the orbital radius from the distance to the central body (you should probably use this as the course target as well). This will keep you from getting falsely small time values if you happen to be inside the orbit of whatever you're chasing. So the real suggestion when ignoring inner orbits is to calculate the intercept to the motion (due to central body motion) of the closest point on the orbital circle to your current position.
The orbit calculations in Aurora are fairly easy as I cheat slightly
. Rather than work out objects moving in parabolas, I store the current bearing of each body from its parent body. Then, based on how much time has passed since the last orbital movement phase, I work out how far around its orbit the planet has moved (time since last move in seconds / orbital period in seconds). I multiply that result by 360 so I know the new bearing then just calculate the position of the body based on a line the length of the orbital radius drawn from the location of the parent body on the new bearing. Using this method, it is very fast to work out exactly where a moon orbiting a planet orbiting a star orbiting another star will be in a specified period of time (assuming you do all the bodies in the right order
I am now loading all system bodies into memory and orbital movement takes place during every sub-pulse. It's faster than I expected, much more accurate and looks cooler too
. A body won't move if the length of the time period means that its orbit changes by less than 1/10,000th of a degree.
For fleets, I now estimate the time taken to reach the target body, based on time required for acceleration from current speed to max fleet speed (if one is set), the time required for deceleration from there to the escape velocity of the planet and any time at top speed (which is sometimes zero because max speed cannot be achieved before deceleration must start). Then I work out where the planet will be after that amount of time and recalculate the fleet time of arrival using that new estimated planetary position as the destination. I run that cycle 10x and by the tenth time the estimated destination is changing by a very small amount.
However, a new problem arose despite the above. In some cases, the fleet was still chasing the planet and for the life of me I couldn't figure out why. After a few days of frustration it finally occured to me to check the speed at which the target planet (Earth in this case) was moving along its orbital path. It turns out to be 30 km/s. As the fleet was slowing to Earth's escape velocity of 11.2 km/s, Earth was moving away from it. The fleet would then speed up to catch Earth and repeat the cycle.
Which creates an entirely new problem, especially as the fleet could be approaching a body on a orbital path from any direction. The fleet would therefore be closing on the planet at different velocities depending on direction and that closing velocity would change during the approach as the planet moved on its curved path. Not to mention that a moon, its parent planet and the parent star could all be moving in different relative directions to the approaching fleet, which severely complicates any calculations regarding suitable approach speeds. Aaagh!
Therefore, instead of using a planet's escape velocity as the maximum speed at which a fleet can successfully enter close orbit, I think I will use the maximum speed of the planet through space plus its escape velocity. For example, that means that if you reach Earth at a speed of less than 41.2 km/s, you can enter orbit. In this case of approaching from behind the planet, this is perfectly reasonable in reality. In the case of approaching from ahead then I am breaking the laws of physics a little
. Of course, for moons this means including both the orbital speed of the moon and its parent planet, because for part of the moon's orbit its motion through space (from the perspective of a stationary observer - well stationary relative to the system anyway) it would be moving at a combined velocity of its own orbital speed and that of the parent planet. As you can imagine, this is more complicated for planets and moons orbiting the outer stars in a multiple system. Even given all that, I still think this is still a reasonable compromise given the horrendous maths involved in the alternative.
Because of this solution, I will compensate by also adding the motion of a moon/planet to the momentum of a ship that leaves orbit. This will require me to optionally display the current bearing and speed of all system bodies (including any motion from parent planets/stars) so that the player can see the tactical implications of leaving orbit at any given time. A combination of entering and leaving orbit could be used for slingshot maneuvers, which while unrealistic in some circumstances will at least provide the flavour of such manoeuvers and create an extra tactical dimension.
Steve
