Intersection with object on ellipsoid path

Question

In my game I need to calculate the time when a ship arrives at the planet. What is known is ship's starting position and velocity (which is constant), and planet’s position at a given time (it follows an elliptic path). To be more specific:

Regarding the ship:

$x_0, y_0$ - ship's initial position
$v$ - ship's constant speed

Regarding the planet:

The planet’s position is given by $$\begin{align}x(t)&=a\cos{\omega_0t}\\y(t)&=b\sin{\omega_0t}\end{align}$$ with $a\ge b\gt0$ and $\omega_0\gt0$. I.e., the planet’s path follows an ellipse in standard position with a phase of $0$.

Now, what I need to know, is at what point does the two positions intersect (ship's position and planet's position). So I am either looking for the intersection time $t_i$, or an intersection point $p_i$, or angle $\alpha$ at which the ships should be fired, or distance $d_i$ of intersection point from starting point $x_0, y_0$. Any of these 4 things should do (figuring one from the other is trivial, of course). I am looking for the closest such intersection point of course (the smallest such $d_i$ or $t_i$ of all possible ones).

In other words, I need to know at what angle should I send the ship from $x_0, y_0$ so that it will arrive at the planet in smallest possible time.

UPDATE 29.12.2016

I am pausing my work on this problem, since I've already spent 5 days on it and am really tired of it. I've tried all kinds of approaches, but the code is buggy and this problem is much more complex than I first thought. I will finish it at some point in the future, but right now there are other aspects of the game I have to implement. So I'd like to help everyone who's contributed, I'll upvote/accept your answers once I get to finish the implementation (I'll present my final algorithm at that point and post it for anyone to use it). Thanks again to everybody for now!

(P.S., if interested, this is how the thing works when it works: https://youtu.be/KjQCOkWVIvg)

Answer 1

I think that finding $t_i$ may be the simplest way to approach this problem, although no approach is really simple.

The distance from your spaceship's starting position to the position of the planet at time $t$ is $$ \sqrt{(a\cos(\omega_0 t) - x_0)^2 + (b\sin(\omega_0 t) - y_0)^2}. $$ Assuming the spaceship starts at time $t_0,$ the distance from the the spaceship's starting position to the spaceship's position at time $t$ is $v(t - t_0).$ (I introduced the parameter $t_0$ because it was not clear that you wanted the spaceship to start at the instant when the planet passed the point $(a,0).$ If you did want the spaceship to start at that exact instant, just set $t_0=0$ in all the equations; it will simplify them a bit.)

In order for the spaceship to intercept the planet at time $t_i,$ the spaceship and the planet must be at the same distance from the spaceship's starting point at that instant. that is, $$ \sqrt{(a\cos(\omega_0 t_i) - x_0)^2 + (b\sin(\omega_0 t_i) - y_0)^2} = v(t_i - t_0). $$ So if we define a function $f$ by $$ f(t) = \sqrt{(a\cos(\omega_0 t_i) - x_0)^2 + (b\sin(\omega_0 t_i) - y_0)^2} - v(t_i - t_0), $$ one way to find a time when the spaceship can intercept the planet is to solve for $t$ in the equation $f(t)= 0.$

Unfortunately, I'm fairly sure there is no closed-form solution for this equation, at least not using the functions you would have available in a typical programming environment. So the only way to solve the equation is by numeric methods--basically, making guesses and refining the guesses until you get "close enough" to the exact solution.

The distance of the planet from from the spaceship's initial position, $(x_0,y_0),$ periodically increases and decreases. The rate of increase has some maximum value $u_\max$. If $v$ is greater than or equal to $u_\max$ then there is exactly one solution to the equation $f(t)= 0.$ The planet can never increase its distance from $(x_0,y_0)$ faster than it is traveling, and it never travels faster than the speed $a\omega_0,$ so if $v \geq a\omega_0$ the problem is slightly simpler than it might be. If $v < a\omega_0$ then you either have to figure out the value of $u_\max$ so that you can determine whether $v \geq u_\max$, or you can decide to solve the problem without making the assumption that the solution to $f(t) = 0$ is unique.

Let's consider first what happens when we know the solution is unique. We know that $f(t)>0$ when $t=t_0$; now find a time $t_1$ such that $f(t_1)< 0.$ Setting $r_0=\max\{a,b\},$ every part of the ellipse is inside the circle of radius $r$ around $(0,0),$ so no part of the ellipse can be further than $r_0 + \sqrt{x_0^2+y_0^2}$ from the starting position of the spaceship. So any value of $t_1$ such that $$ t_1 > t_0 + \frac{r_0 + \sqrt{x_0^2+y_0^2}}{v} $$ will be sufficient. Now we can use one of several ways to find an approximate value of $t_i,$ but the simplest may be the "bisection" method: find the midpoint of your interval of time, compute $f(t)$ at the midpoint, and then change your interval to the interval bounded by the midpoint and one of the two previous endpoints of the interval so that $f(t)$ changes sign between the new endpoints. In other words, split your interval of time in half and replace the old interval with the half-interval in which $f(t)$ changes from positive to negative. Repeat until the interval is so small that it no longer makes a difference where in the interval the solution is (that is, no matter which time you choose, the spaceship will get close enough to the planet at that time that you consider it a "hit"). If a separation of $\delta$ units of distance is "close enough" in space, then a time difference of $\delta/v$ will be "close enough" (possibly even better than "close enough") in time.

In the simple case (unique solution), at this point you're done. The rest of this answer concerns the more complicated "possibly non-unique solution" version of the problem.

If you do not know that the solution is unique, you can still use the bisection method, but if you use it as described above it may not find the earliest solution. It is possible that $f(t)$ has a local minimum that is less than zero, but that $f(t)$ is positive at some time after that. There would then be three (or five or seven or more) solutions before $f(t)$ goes negative for the last time. To be sure that you don't miss the first solution of $f(t)=0,$ you have to find out at what times the local minimums occur.

To find when the local minimums occur, take the derivative of $f(t)$ with respect to $t$ (or have Wolfram Alpha do it for you), and solve for $t$ in the equation in which that function is zero. You only need to find one solution $t=t_2,$ it doesn't matter whether it is the "earliest", and then all other solutions are just $t_2$ plus or minus some multiple of $2\pi/\omega.$

Find the earliest local minimum; if $f(t)<0$ at that time then the solution is between $t_0$ and that time; otherwise try the next local minimum, and the next, and so forth until you find two local minimums such that $f(t)>0$ at the first one and $f(t)<0$ at the second, and then look for a solution between those two times; it will be the earliest solution. (Of course if $f(t)=0$ at any of the local minimums then that's your solution. You could even accept the local minimum as a solution if $-\delta < f(t) < \delta$ where $\delta$ is a distance that you consider "close enough" to count as an interception.)

Answer 2

The parametric equations for ellipse are $$x=a\cos\omega t~~~;~~~y=b\sin\omega t$$ and the parametric equations for line are $$x=x_0+v_0\cos\theta t~~~;~~~y=y_0+v_0\sin\theta t$$ that $\theta$ is path angle with respect to $x$ axis. Then we can obtain positions intersect with these equations.

‎I investigate some cases for this problem.

First Stage: With $a$, $b$, $x_0$, $y_0$ and with prescribed $\theta$ and $\omega$ which are constant here, we want to determine $v$, we have: $$t=\frac{a\cos\omega t-x_0}{v\cos\theta}=\frac{b\sin\omega t-y_0}{v\sin\theta}$$ so $$(b\cos\theta)\sin\omega t+(a\sin\theta)\cos\omega t=x_0\sin\theta-y_0\cos\theta$$ by substituating $\displaystyle\sin\omega t=\frac{2\tan\frac{\omega t}{2}}{1+\tan^2\frac{\omega t}{2}}$ and $\displaystyle\cos\omega t=\frac{1-\tan^2\frac{\omega t}{2}}{1+\tan^2\frac{\omega t}{2}}$ and simplyfing $\tan\frac{\omega t}{2}=k$ we conclude that $$(b\cos\theta)\frac{2k}{1+k^2}+(a\sin\theta)\frac{1-k^2}{1+k^2}=x_0\sin\theta-y_0\cos\theta$$ so $$\tan\frac{\omega t}{2}=\frac{b\cos\theta\pm\sqrt{(b^2-y_0^2)\cos^2\theta+(a^2-x_0^2)\sin^2\theta+x_0y_0\sin2\theta}}{(a+x_0)\sin\theta-y_0\cos\theta}~~~~~~~~(1)$$ or $$t=\frac{2}{\omega}\tan^{-1}\frac{b\cos\theta\pm\sqrt{(b^2-y_0^2)\cos^2\theta+(a^2-x_0^2)\sin^2\theta+x_0y_0\sin2\theta}}{(a+x_0)\sin\theta-y_0\cos\theta}$$ this is the time witch the planet arrives to straight line and for accident, the ship must travel a distance in $t$-time to goes there, or $$v=\frac{x-x_0}{t\cos\theta}$$ if the prescribed $v$ is not equal to this value, the ship and the planet won't collide.

Second Stage: With $a$, $b$, $x_0$, $y_0$ and with prescribed $\theta$ and $v$ which are constant here, we want to determine $\omega$: $$t=\frac{x-x_0}{v\cos\theta}$$ with (1): $$\omega=\frac{2v\cos\theta}{x-x_0}\tan^{-1}\frac{b\cos\theta\pm\sqrt{(b^2-y_0^2)\cos^2\theta+(a^2-x_0^2)\sin^2\theta+x_0y_0\sin2\theta}}{(a+x_0)\sin\theta-y_0\cos\theta}~~~~~~~~(2)$$ like previous stage, the planet must travel a distance on it's path in $t$-time to goes there, otherwise if the prescribed $\omega$ is not equal to this value, the ship and the planet won't collide.

Third Stage: This case is complicated. With $a$, $b$, $x_0$, $y_0$ and with prescribed $\omega$ and $v$ which are constant here, we want to determine $\theta$ which for collide occurs. For this purpose we delete time between equations $$\left\lbrace\begin{array}{c l}\cos\theta=\frac{x-x_0}{vt}=\frac{a\cos\omega t-x_0}{vt},\\\sin\theta=\frac{y-y_0}{vt}=\frac{b\sin\omega t-x_0}{vt}.\end{array}\right.$$ Since I don't know these cases are interest of @Betalord, I finish my notes.

Answer 3

You have the ellipse of the planets trajectory $$ \left( \frac{x}{a} \right)^2 + \left( \frac{y}{b} \right)^2= 1 $$ (Note: I use a coordinate system, where the origin is not on the ellipse).

Further you have a start point $u_0 = (x_0, y_0)$ where the straight line of the ship's trajectory $$ u = u_0 + v t \quad (*) $$ will originate from.

A possible trajectory of the planet on the ellipse is $$ (x, y) = (a \cos \omega t, b \sin \omega t) \quad (**) $$ An intersection at some time $t \ge 0$ will lead to the equation: $$ (a \cos \omega t, b \sin \omega t) = (x_0 + v_x t, y_0 + v_y t) $$ or $$ (f(t), g(t)) = (x_0 + v_x t - a \cos \omega t, y_0 + v_y t - b \sin \omega t) = (0, 0) = 0 \quad (***) $$

Update:

Knowing now that the direction of $v$ is a degree of freedom, we can proceed. $$ v = (v_x,v_y) = \lVert v\rVert (\cos(\phi),\sin(\phi)) $$ We have everything together to model the intersection problem in GeoGebra. The given information which determines the problem instance is $$ x_0, y_0, a, b, \omega, \lVert v\rVert \quad (P) $$ the variables are $$ \phi, t $$ In the image below point $A$ is the start location of the ship. $B$ is the ship position at $t=1$. The green ray from $A$ through $B$ is the trajectory $(*)$ of the ship. The red ellipse is the trajectory of the planet $(**)$. The intersection points of $(*)$ and $(**)$ are $C$ and $D$. Point $S$ is the position of the ship at time $t$, $P$ is the position of the planet at time $t$.

In the interactive version (see link below), one can twiddle with the $t$ slider and watch the ship and planet move on their trajectories.

The goal is to find a value for $\phi$ which will lead to a common position for ship and planet at some time $t\ge 0$.

The key are the two functions from $(***)$, shown as $f(x)$ (purple line) and $g(x)$ (blue line). In the interactive version one notes that these graphs change depending on the value of $\phi$. To solve the problem one has to find a value for $\phi$ which causes the graphs of $f$ and $g$ to both have a root (zero value) at some $x\ge 0$. In the image you will notice two such arguments around $x=1.5$ and $x=4$. They correspond to the intersections $C$ and $D$. The present time $t$ of the scene is visualized by the vertical black line around $x=3.6$ (see $t$ slider as well), so ship and planet are shortly before impact at $D$, $t=4$.

In general, depending on the problem instance $(P)$, there might be no, one or two intersections.

To implement this procedure on a machine, we have to replace the human fiddling with the slider for $\phi$, looking for simultaneous roots of $f$ and $g$ by some numerics, solving equation $(***)$, thus $$ F(\phi, t) = \begin{pmatrix} x_0 + \lVert v\rVert \cos(\phi) t - a \cos(\omega t) \\ y_0 + \lVert v\rVert \sin(\phi) t - b \sin(\omega t) \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} $$

Common methods are root finding by Newton-Raphson iteration or simple bisection, or reformulation as fixed point iteration.

If we look at the square of the distance between ship and planet we get this to $$ q(\phi, t) = (x_0 + \lVert v\rVert \cos(\phi)t - a \cos(\omega t))^2 + (y_0 + \lVert v\rVert \sin(\phi)t - b \sin(\omega t))^2 $$

Update:

Here is a link to the interactive version: GeoGebra

Update:

I explained the Newton-Raphson iteration to my friend Ruby (see here), and she performed it a couple of times for me.

irb> newton($u0, 1e-10, 20)
x0 = 1.0
y0 = 2.0
a  = 2.0
b  = 1.0
v  = 1.0
w  = 1.0
0: u=Vector[1.0, 2.0] f(u)=Vector[2.9128982848305642, 2.7736445427901115]
1: u=Vector[0.382829472675263, 0.3248271541942971] f(u)=Vector[-0.5940984322993181, 1.8021930994105337] d=1.8975913455808797
 .
 .
20: u=Vector[-94.18834695323662, 1959.9371619837884] f(u)=Vector[1955.6477996333003, 118.82040799198684] d=1959.2540941812401
maxstep = 20 reached, fail
=> nil

So it did not converge, when we started to look at $u_0 = (x_0, y_0) = (1,2)$ for a root. On the other hand $u=(\phi, t)$, so it was not a good idea in the first place.

irb> newton(Vector[-1,4], 1e-10, 20)
 .
 .
20: u=Vector[138.25481797929191, 454.49648517971787] f(u)=Vector[456.3793855184854, 12.384069736222514] d=456.5473783841737
maxstep = 20 reached, fail

Another miss.

irb> newton(Vector[-1,1], 1e-10, 20)
0: u=Vector[-1, 1] f(u)=Vector[0.45969769413186023, 0.317058030384207]
1: u=Vector[-1.5668891620557184, 1.0077918056648623] f(u)=Vector[-0.06352118626769587, 0.14656055765459897] d=0.15973396058722875
2: u=Vector[-1.664180141632932, 1.103100791097277] f(u)=Vector[-0.004523100381582101, 0.009095949944068482] d=0.010158481306123844
3: u=Vector[-1.67043452556281, 1.1098339180715397] f(u)=Vector[-1.951140039613275e-05, 4.57909054815131e-05] d=4.9774509241529196e-05
4: u=Vector[-1.670469710132332, 1.1098684183520295] f(u)=Vector[-6.100806526632141e-10, 1.337344013307984e-09] d=1.4699276896108456e-09
5: u=Vector[-1.6704697111186804, 1.1098684193565171] f(u)=Vector[-2.220446049250313e-16, 2.220446049250313e-16] d=3.1401849173675503e-16
d=3.1401849173675503e-16 < eps=1.0e-10, success

Yeah, this one worked!

irb> newton(Vector[0,0], 1e-10, 20)
0: u=Vector[0, 0] f(u)=Vector[-1.0, 2.0]
Jacobian Matrix[[-0.0, 1.0], [0.0, -1.0]] is not regular at u=Vector[0, 0]!
=> nil

Ouch, another bad start. This reminds, that one has to come up with a good start value, close to the suspected root, for this method to work.