# Intersection with object on ellipsoid path

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)

• Are ship and planet in a plane? Is ship follows the planet or moves in straight line? – Nosrati Dec 25 '16 at 18:23
• Ah yes it's a 2D plane, and ship is traveling in a straight line. – Betalord Dec 25 '16 at 18:44
• Is the full velocity vector given (length and direction)? – mvw Dec 25 '16 at 18:58
• No, only length is given... direction is already a solution to the problem! – Betalord Dec 25 '16 at 21:42
• I don’t entirely understand the conditions of the problem. First, is the origin at the center of the ellipse or at one of the foci? Second, is $\omega$ really a constant angular velocity? That’s rather different from its being, say a coefficient in the parameterization $t\mapsto(a\cot\omega t,b\sin\omega t)$. I would find it helpful if you included an explicit formula for the position of the planet. – amd Dec 26 '16 at 3:39

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.)

• @mvw Once you find one of the (finite number of) solutions of the equation $f(t)=0$, that gives you the location where the planet will be at the time of interception. Set $\phi$ to aim the spaceship at that point. – David K Dec 27 '16 at 18:19
• @mvw An alternative approach is to solve for $\phi$ in $g(\phi)=0$ where $g(\phi)$ is some function that is zero if interception occurs; for example, we might let $g(\phi)$ be the distance of closest approach between the spaceship and planet. But it seemed to me that's a harder equation to set up, let alone to solve. – David K Dec 27 '16 at 18:21
• Thanks. I think I got your idea now. – mvw Dec 27 '16 at 18:37
• Nice work, you actually explained the algorithm, which I think will be useful for me when I get around to it again (I'm pausing right now, have spent 5 days implementing my own algorithms that I thought were really intelligent but in the end none of them works well). I've read your post a couple of times and liked the idea about doing the derivative and then searching for a single minimum since all others are just multiples of $2\pi/\omega$ (didn't think of that). – Betalord Dec 29 '16 at 9:49

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.

• And how do we obtain positions intersection with these equations? – Betalord Dec 25 '16 at 21:40
• @Betalord With some calculations and delete $t$ between these equations, you will get an equation like $A\sin\omega t+B\cos\omega t=C$ and then take $\tan\frac{\omega t}{2}$ to solve it! – Nosrati Dec 25 '16 at 21:52
• Well, that is my problem :) I don't know how to find an equation for t. My math knowledge is limited. I know how to set up the system of equations as well, but not how to solve it, unfortunately. If you can help me with that, I'd be most thankful. – Betalord Dec 25 '16 at 22:21
• You might rewrite your second set of equations. They look like the factor of $t$ is part of the argument to the trigonometric functions, which it isn’t. – amd Dec 26 '16 at 8:48
• Combining these equations produces $\tan{\omega t}={avt\sin\theta+ay_0\over bvt\cos\theta+bx_0}$. I don’t know of an analytic solution of this equation for $t$. If you do, please edit your answer to include it. – amd Dec 26 '16 at 8:50

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.

• Yes I came there too, but I am searching for a solution, not a system of equations. Can you solve this system for t? – Betalord Dec 25 '16 at 21:40
• OK perhaps I was not clear... I can't use some external program to find numerical solution to the problem. I'm writing a game and I need analytical solution so that I can code it in my game. For example, I need an equation like t=..., or something similar. My problem is that I don't know how to find an analytical (symbolic) solution to your equations. Can you? – Betalord Dec 25 '16 at 22:18
• @Betalord I hope you fiddled a bit with the online interactive model. It should give you an idea how the parameters lead to problems with solutions or not. – mvw Dec 26 '16 at 16:44
• I've implemented the Newton method today and it didn't work... it always returned INF. The problem is that this function is not a good one for Newton's method. It goes up and down (since it has trigonometric functions) and hence tangent doesn't hit y=0 and so ends up in infinity. One should already have a good initial value in order for the Newton to work. Currently I am thinking about some other methods... but I guess the way to go is a numeric way. I think an analytic solution doesn't exist (or it's just very hard to derive). – Betalord Dec 26 '16 at 17:29
• I documented my version (code and sample runs) and it gives a root if one is in a proper neighbourhood. An analytic solution ought to be complicated because there are not always solutions. Alternatives are bisection and fixed point iteration. – mvw Dec 27 '16 at 10:54