Mercator Sailing: Course and Distance

[Dieter]

(08-24-2018 07:21 PM)Gene222 Wrote:  I modified Eddie's program where South latitudes and West longitudes must be entered as negative numbers in degrees, but I ran into a problem when the two points have the same latitude, such as

Lat 40 N, Long 25 W
Lat 40 N, Long 30 W

In this case, the course and distance equations give division by zero error messages. How do you calculate the course and distance on a parallel?

I wrote a Mercator Saling program for the 35s and encountered the same problem. However, a bit of calculus yields the following result:

Code:

                                           1 - e²·sin²(lat)
distance = 60·|long2 - long1| · cos(lat) · ----------------
                                                1 - e²

The course c then is 0° (exactly East or West). Or 90° / 270° true course.

Here "e" is the eccentricity of the assumed ellipsoid (around 0,082). This is not explicitely considered in the approximate formula given in the book (as used in Namir's program).

But we can do better: calculate the meridional parts (M) as follows using hyperbolic functions:

Code:

M = 60·180/pi · [artanh(sin(lat)) – e·artanh(e·sin(lat))]

You can also do it without hyperbolics:

Code:

M = 60·180/pi · ln[tan(45° + lat/2) · ((1–e·sin(lat))/(1+e·sin(lat))^(e/2)]

BTW, the constant 7915,7045 in the approximate formula is 60·180/pi · ln10 due to the use of the base 10 logarithm instead of the natural log.

Both formulas are more accurate and allow implementing different ellipsoids. I prefer the one with the two artanhs since it is more compact and ...elegant. ;-)

For "e" you may use these values:
Clarke 1880 ellipsoid: e = 0,082483399
WGS84 / GRS80 ellipsoid: e = 0,081819191

For your example the distance is 230,7 nmi (for both ellipsoids the results differ by less than 0,02 nmi).

Check: use your program "as is" and make the latitudes slightly different, e.g. 40° and 40,0001°. This should give (nearly) the same result.

Remark: instead of calculating the arctan you can determine the angle with a R→P conversion (automatically selects the correct quadrant) followed by mod 360° to get a positive angle. This way the four IF cases are not required (BTW, are you sure the last two are both correct?).

Dieter