Cruise control ability to navigate around the planet

[Alastor-git]

If the designated target is partially or completely around the planet, the cruise control would navigate the ship through the gravity field or even through the planet.
It would also try to cover large distances laterally.
Instead, it should point and navigate along a great circle around the gravity field and only point towards the planet for final approach.

For now, ignore the suggested algorithm. I think I made a mistake, it's not always a great circle.
I know a better way to find the right algorithm, but I have to formalize it.
I'll come back and edit with the correct one.

To navigate along the great circle, I suggest the following algorithm :

• Let's call up = (xp, yp, zp) the planet coordinates. u1 = (x1, y1, z1) the ship initial coordinates, u2 = (x2, y2, z2) the final coordinates.
• u1 = u1 - up, same for u2 to go from absolute to relative coordinates.
• In the (x, y, z) reference frame, a circle of radius r has the parametric equations :

x = r cos(theta)
y = r sin(theta)
z = 0

• The great circle is a rotation of this circle around the center of the planet to align with the plane that contains the vectors pointing to the starting and end point. It's along this circle that we find the shortest path on a sphere. To find it, let's find the reference frame of said plane to then rotate the circle in the desired plane.
• One special case where the rest of the computation would fail is if u1/|u1| = - u2/|u2|. But in that case, the equator plane is the right choice, so we can skip most of the following.
• We decide that the x axis should be defined by the starting vector u1, so that the starting parameter is 0.
• The coordinate vectors of the new plane are :

ex = u1/|u1| = (exx, exy, exz)
ey = ez ^ ex = (eyx, eyy, eyz)
ez = u1/|u1| ^ u2/|u2| = (ezx, ezy, ezz)

To rotate objects in the new plane, we must use the transformation matrix :

 (exx, eyx, ezx
  exy, eyy, ezy
  exz, eyz, ezz)

• Using this, the transformation of the circle is :

x = r ( exx cos(theta) + eyx sin(theta) )
y = r ( exy cos(theta) + eyy sin(theta) )
z = r ( exz cos(theta) + eyz sin(theta) )

• We now know the coordinates of the circle we will follow. We know we start at theta1 = 0, but we don't know where to stop. To find theta2, we must stop for x = x2, y = y2 and z = z2. One of the obvious elements of the solution is r = r2. The rest is a linear system, yielding : theta2 = Atan2(y2 exx - x2 exy, x2 eyy - y2 eyx) it should be in the [-Pi;Pi] range. .
• Since theta1 = 0, we know that whatever theta2 is, going directly from theta1 to theta2 is the shortest path.
• We parametrize the path by s in [0 ; 1] with :

theta(s) = s theta2 # theoretically, theta2 - theta1, but theta1 is 0
r(s) = r1 + s (r2-r1)
x(s) = r(s) ( exx cos( theta(s) ) + eyx sin( theta(s) ) )
y(s) = r(s) ( exy cos( theta(s) ) + eyy sin( theta(s) ) )
z(s) = r(s) ( exz cos( theta(s) ) + eyz sin( theta(s) ) )
u(s) = ( x(s), y(s), z(s) ) + up

• u(s) gives us the ideal trajectory to follow.

I may have a look at implementing this later, the main challenge is to have a guidance algorithm that allows us to follow a trajectory. But at least we have the algorithm written down for said trajectory.