Shortest path on a sphere direction descision

Question

I'm trying to write an algorithm that does the following.

Given a current position (in Azimuth and Inclination) and a target position (again in A, I) in what direction should I travel to travel over the shortest path. The return value could be something like a vector A = -1, I = +0.5, that I can then scale for step size/time.

The shortest path can be found by using a great circle, this is easy to visualize, but it's hard to implement like above because my coordinate system isn't continuous.

My coordinate system is as followed (imagine standing in front of the sphere)

The azimuth is 0 ~ pi when traveling along the equator along the front side, it's 0 ~ -pi when traveling along the equator along the rear side.

The inclination is 0~+pi when traveling from the top to the bottom of the sphere.

So given this non-continuous coordinate system, how do I create a decision function that says 'increase A' to travel over the shortest path?

Answer 1

You have a couple of alternatives. The first is to use a Haversine formulation. There is some Javascript source code here. It requires using more traditional lat / lon where the equator is at 0 latitude and the poles are at +/- π or +/- 90° latitude (depending on your units) and longitude is in the range [-180°, 180°) or [-π, π) again depending on your units. You can repeatedly find the midpoint until you have an approximate path that suites your needs. The azimuth / inclination vector would just be the difference in lat / lon between two adjacent points, though over time this will likely induce an error if you repeatedly apply those lat / lon deltas to the location of your agent.

Another approach that may work well for you is to transform your spherical coordinates of your starting and ending location to cartesian coordinates, call them points u_b and u_e for beginning and end points. The normal vector v of the great circle connecting the two points is the cross product of the two (i.e. v = u_b x u_e) and the angle θ is just the arccosine of the normalized inner product (ie. θ = cos^-1( (u_e ∙ u_e) / (|u_b||u_e)). You can then use quaternion rotation and iterate from 0 to θ about the vector v to actually navigate the path. With this approach, the actual instantaneous vector at some point p along the path is just the p x v, or you can just approximate this by using the cartesian difference between two adjacent points along the path.