Determining the first future intersection possible between ships and a planet

Question

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

Answer 1

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)² = (p - cos(t))² + sin²(t)
        = p² + 1 - 2p cos(t)

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

T(t)² = (s(t - t_0))²
        =s²t² - 2s²t*t_0 + s²t_0²

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

s²t² - 2s²t*t_0 + s²t_0² = p² + 1 - 2p cos(t)

Which we can rephrase as:

cos(t) = (1/2p)(-s²t² + 2s²t*t_0 - s²t_0² + p² + 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.

Answer 2

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.