Get heading between two points

[RisingFury]

I'm trying to calculate the heading between two points on Earth. I have both their equatorial positions, however, either my concept, my formula or my implementation is wrong. The angle is off by a few degrees...

I know Orbiter doesn't have a problem calculating heading - if you set a target base, the little triangle at the top points to it.

My question is: How does Orbiter do that accurately? Here's my formula:

I designated x as longitude, y as latitude.

My guess is that the sine and cosine functions aren't 100% accurate.

I also get two singularities in the formula - that is to be expected, since I'm dividing with 0 when I'm directly north or south to the target. What's odd is that the singularity doesn't appear at one particular angle, but ranges over a degree in width or so.

 

[agentgonzo]

I doubt that it's a problem with the lack of accuracy of the cos and sin functions as they should be accurate to many many decimal places.

You seem to be overcomplicating the derivation though, and possbily making a mistake somewhere.

Using the same notation as in your diagram, we know the values of a c (from latitude of target position), b (from latitude of ship position) and if we label the angle A as being the angle SNT (opposite side a on the triangle), then we know that this angle A is the difference in longitudes of the ship and target.

Therefore, using the law of spherical triangles, you can get:

cos(β) = (cos(c) - cos(a)cos(b)) / (sin(a)sin(b))

as you have correctly determined. The only unknown in this is cos(a).

This can be determined by the same spherical triangles rule:

cos(a) = cos(b)cos(c) + sin(b)sin(c)cos(A)

You can just plug this directly into the first equation to get an equation for β.

I think that this is what is used to derive the Haversine function, which you could just use directly.

I don't know how orbiter does this, but this assumes that the Earth is spherical (which it isn't). You can modify this in a simple case to take account of the oblateness of the Earth. For more accurate results (which you most probably don't need), look up a paper written by Vincenty in 1975 that gives a numerical solution of very high accuracy to the problem, but I would imagine that the spherical approximation works fine.

---------- Post added at 13:56 ---------- Previous post was at 13:39 ----------

Vincenty formula if you fancy being overcomplicated:http://www.movable-type.co.uk/scripts/latlong-vincenty.html

 

[RisingFury]

Escapetomsfate answered me in PM, but I'm putting this here as a reference.

Try this:
tc1=mod(atan2(sin(lon2-lon1)*cos(lat2),
cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon2-lon1)),
2*pi)

variables:
lon1 = longitude of point A
lat1 = latitude of point A
lon2 = longitude of point B
lat2 = latitude of point B
tc1 = direction of point B from point A (angle east of north)

Thanks!

 

[Kozak]

My question is: How does Orbiter do that accurately? Here's my formula:

Orbiter AFAIK does that in Mercator projection, so it gives out loxodrome instead of geodesic curve. That gives an error around the mentioned one-two degrees.
EDIT: escapetomsfate actually does the same, so your method will be more precise. If you get rid of singularities, of course.