Determining the first future intersection possible between ships and a planet

Somewhat of an extension of this.

The Problem

I have a ship being launched from a planet, and I would like it to travel at a constant speed and direction such that it intersects another planet the soonest that it can.

(Beautiful diagram courtesy of Ausa)

I am ignoring initial velocity, end velocity, and gravity. The angular velocity of the planets and the velocity of the ship is constant. The physics that I am trying to achieve need not be realistic.

Using some quick maths I created a function! The function determines the smallest fly-by distance from the ship to the planet over time if the ship were to take the optimal route (a linear route since there's no gravity). So by determining when the fly-by distance is 0 (when the paths of the ship and planet intersect) I can work backwards to easily determine the angle of which I must launch the ship at to achieve this.

I created an implementation of the Bisection Method to find the first root of the function but depending on the variables of the function, the bisection method isn't the best solution to this problem. So now I'm trying to solve it mathematically! This is essentially what I need to solve:

Any solutions or ideas are greatly appreciated, thank you!

TL;DR:

This means I want to find the x-intersection of the function nearest to the origin.

Desmos

If you're interested in playing around with this function, I graphed it in Desmos https://www.desmos.com/calculator/tdtzt1gqkl

r_0 radius from the planet to the sun ("magic" units)

a_0 current angle of the planet in relation to the run (clockwise radians)

s_0 angular speed of the planet (clockwise radians per second

x_1 current x-position of the ship ("magic" units)

y_1 current y-position of the ship ("magic" units)

s_1 speed of the ship ("magic" units per second)

x_p planet x-position at x time ("magic" units)

y_p planet y-position at x time ("magic" units)

d_x delta x-position from ship to planet at x time ("magic" units)

d_y delta y-position from ship to planet at x time ("magic" units)

d_net net delta position from ship to planet at x time ("magic" units)

v_x optimal x-component of ship velocity ("magic" units per second)

v_y optimal y-component of ship velocity ("magic" units per second)

x_s ship x-position at x time ("magic" units)

y_s ship y-position at x time ("magic" units)

f(x) net distance from planet to ship at x time

I don't think a purely analytical solution to this is possible.

Any way we slice it, we have part of the problem in polar coordinates and another part in cartesian coordinates. So we end up trying to find an intersection between a line or curve and a sinusoid/periodic trig function. If the line were a constant we could use the arcsine/arccosine/arctangent functions. But unfortunately we're dealing with curves or lines with slope so no standard function gets us what we need.

In a previous answer I showed how we can solve this to within the needed accuracy by binary search. But since you've specified you don't want to solve it that way, let's see if we can reframe our problem differently.

First, to simplify our math, let's apply some (reversible) changes of coordinate system and units.

• Define the center of the planet's orbit as the origin

• Define the radius of the orbit as 1 (effectively our distances are now measured in AU)

• Define our unit of time as the duration it takes the planet to travel 1 unit of space (so now 1 planetary year = 2 pi time units)

• Define time t=0 as the most recent moment that the planet passed between the spaceship's start position and the center of the orbit, and time t_0 is the "start" time when our ship begins moving from that location

This lets us standardize the situation like so:

So our parameters are:

• t_0, the time forward from which we want to find solutions:

t_0 = angleDifference(planetStartAngle, atan2(shipStartPosition - orbitCenter))
if(t_0 < 0)
t_0 += 2 * pi


This is just the clockwise (orbit-wise) angle in radians between the ship's start position and the planet's, as viewed from the center of the orbit. Starting from perfect alignment with the ship at t = 0, the planet needs to orbit through angle t_0 before it arrives at its start position relative to the ship.

• p, the scalar distance of the ship's start position from the center of the orbit:

p = length(shipStartPosition - orbitCenter)/orbitRadius

• s, the speed of the ship in our new spacetime coordinate system's scale:

s = shipSpeed * 2 * pi / (planetAngularSpeed * orbitRadius)


We're interested in the function L(t), the distance from the ship's start position to the planet.

L(t)2 = (p - cos(t))2 + sin2(t)
= p2 + 1 - 2p cos(t)

And we want to set this equal to the distance the ship has traveled by time t > 0,

T(t)2 = (s(t - t_0))2
=s2t2 - 2s2t*t_0 + s2t_02

At the times of intercept, T(t) = L(t) so...

s2t2 - 2s2t*t_0 + s2t_02 = p2 + 1 - 2p cos(t)

Which we can rephrase as:

cos(t) = (1/2p)(-s2t2 + 2s2t*t_0 - s2t_02 + p2 + 1)

Which is "just" an intersection between a parabola and a garden-variety cosine wave (unscaled, unshifted).

Take the parabola side of that equation and call it P(t). We know we can't possibly have an intersection before P(t) drops to less than or equal to 1, so we can use the garden variety quadratic formula and find the greater of the two solutions P(t_low) = 1 to get a lower bound on t.

We know the cosine wave rises back up to 1 by t_high = ceil(t_low/2pi)*2pi, but the parabola is strictly decreasing, so our earliest intersection must lie between t_low and t_high, at most one wavelength apart.

We can also find the greatest solution for P(t) = -1 to get a potentially tighter bound on t_high

Unfortunately this seems to be about as far as a pure analytical approach can take us, and we're again forced to use some type of iterative root finding method in this span, but at least we've reduced it to about the simplest case we could hope for. There could still be multiple roots in this span, but our bounds are tight enough that this should be rare.

If your planet orbits much faster than your ship moves, the parabola will be nearly flat in this span and you can closely approximate it as a constant, using acos to find a very close first guess between t_low and t_high. Or if your ship travels much faster than your planet, the parabola will be nearly vertical the solution will be quite close to t_low. For intermediate cases like the example above, you might be able to approximate each curve as a straight line and iteratively refine from there.

• Couldn't you find the intersection point of the ship's path and the planet's orbit circle? Then you could check if the time the ship needs to get to those positions is the same as the planet needs. I might've misunderstood something though. – Bálint Nov 22 '17 at 7:28
• An acceptable start DM, but actually a bit too complex. One could move to isometric space, ie. making the parabola a straight line and omitting cos() completely (as sin and cos work hand in hand). But one always ends up with sin(mT) = nT, which is the annoying part. – Stormwind Nov 22 '17 at 10:02
• @Bálint no, the ship doesn't yet have an absolute direction. We need to choose its direction such that, moving at its speed, it will reach a point in the orbit at the same time the planet does (effectively finding an intersection between the ship's future cone and the planet's future helix). Guessing a direction then checking whether the planet is ahead of or behind us in the orbit when we get there could be used for another style of iterative search though. – DMGregory Nov 22 '17 at 11:34
• I wonder if there are cases where the most effective way to get to the planet the fastest is to just sit still for a certain amount of time, and then travel to the location? Or to travel to a point where the planet will be, then stop and wait for the planet? – Tim Holt Nov 22 '17 at 18:46
• @TimHolt in any such case, we could instead aim a little "upstream" in the orbit and travel a little longer to intercept it earlier in its path, instead of twiddling our thumbs. Think of the ship's potential positions as a circular ripple expanding continuously over time, or a widening cone in spacetime. There is always a point where the spacetime helix traced by the planet intersects this cone for the first time, and if we steer to that point we'll intercept it at the earliest possible moment, with no waiting required. – DMGregory Nov 22 '17 at 18:50

You can solve this by expressing the position of the ship using parametric equations, expressing the position of the planet using the equation for one side of a circle, and solving a system of equations.

The ship's x position can be expressed as x = vcos(theta)t, where theta is the absolute value of the ship's angle below the horizontal.

The ship's y position can be expressed as y = vsin(theta)t.

The planet's x-position can be expressed as sqrt(1 - (y - k)^2) + h, where k is the translation of the circular orbit from the origin of the coordinate system.

Combining these equations yields: vcos(theta)t = sqrt(1 - (vsin(theta)t - k)^2) + h

To make this a system of equations, the arc length traveled by the planet from a_0 to a_0 + theta is equal to s_0t(2 * pi * r). Due to the same-side angles theorem, this theta is the same as the ship's theta below the horizontal.

The arc length can be given by the equation: Assume o is theta for the rest of this answer.

integral(from 0 to (a_0 + theta))(sqrt(1 - (dr/do)^2)dx - integral(from 0 to
a_0)(sqrt(1 - (dr/do)^2)dx.


The equation for the circle is r^2 - 2rr_0cos(o - o_2) + r_0^2 = a^2.

(r_0, o_2) is the location of the center of the circle

dr/do = (2r_0rcos(o - o_2))/(2r + 2r_0sin(o))


Finally,

integral(from 0 to (a_0 + theta))(sqrt(1 - (dr/do)^2)dx - integral(from 0 to
a_0)(sqrt(1 - (dr/do)^2)dx = s_0t(2 * pi * r)


The above equation most likely will require the trapezoidal rectangular approximation method or Simpson's approximation method in order to solve the integrals.

After solving for t, substitute t into the two ship-position equations to get the x and y positions of the collision.

• That won't work, i'm afraid. You'd need to solve 2 things (theta and t) with one equation system. But if you can prove it works, i'll change my mind of course :-). – Stormwind Nov 22 '17 at 9:48
• @Stormwind You can find theta just by finding the inverse tangent of the ship's y distance passed during the first second over the ship's x distance passed during the first second, but by storing the coordinates in those two different times, and not by using the equation. – Cpp plus 1 Nov 22 '17 at 13:54
• The ship's direction is unknown - that's one of the variables we're solving for. Depending on where the planet is in its orbit, we might want to launch in a different direction to get the earliest possible interception. – DMGregory Nov 22 '17 at 14:36
• @DMGregory I have another idea which eliminates the need for an angle, but requires that you know how many radians around the center of orbit the planet is in. – Cpp plus 1 Nov 22 '17 at 16:09
• This gives the correct solution only if the earliest interception point is on the nearest or furthest point of the orbit. It does not handle cases where we intercept the planet at other parts of its orbit. – DMGregory Nov 22 '17 at 16:13