Subject: Re: Formula for dis. between two points on the Globe
From: Jan van Puffelen (puffelej@XXX.XXX)
Date: Mon Oct 21 1996 - 18:08:37 EDT
Vincent (a dutch name?):
At 11:33 20-10-96 -0400, you wrote:
>Hello,
>
>Is there somebody that has the Formula to calculate the distance
>between two points on the globe.
>
>eg. 52 29.66N and 52 32.86N
> 005 02.45E 004 40.67E
A very Dutch position!
>
>My GPS tells me that it is 13.6nm CTS106
>
>
Yes there are:
1. The straight line course (sailing a constant course from source to
destination) als called "loxodrome"
2. The shortest distance between two points (but this requires a constantly
changing course and can lead to unacceptably high latitudes).
3. The composite track (a good compromise between the two).
I give the loxodrome course in the form of a small BASIC program, as
developed for the CASIO FX7602P, a handheld computer with a small LCD
screen. This machine could store up to 10 BASIC programs and has served me
well over the years:
____________________________________________________________________________
____________________
LOXODROME 218
____________________________________________________________________________
____________________
1 INP "SR",R,"SS",S,"V",V,"HH",U
SET F1
2 PRC 1,D-B Normalize difference in longitude
RPC X,Y
Z=180/pi*LN ((TAN C+1/COS C)/(TAN A+1/COS A)) Calculate integral of
difference in longitude on the spheroid; an ellipsoid is a bit overdone
3 RPC Z,Y
PRT "VH=";X/Z*(C-A)*60 distance
K=Y-ASN (S/V*SIN (R-Y)) corrected for current
4 PRC V,K course vector
Z=8*U/V/V leeway angle
K=K-Z correct for leeway
I=X
J=Y
PRC S,R current vector
K=Y-360*INT (K/360) normalise course
5 RPC X+I,Y+J vector addition
PRT "BV=";X resulting speed including current
SET F0
PRT "DR=";Z,"WK=";K leeway and course
____________________________________________________________________________
_____________________
Distance and course are calculated between any two points on earth with
coordinated (Lat, Long) in (A,B) and (B,C) following a straight line on the
map (Mercator Projection) i.e. with a constant course (Loxodrome). The
course is corrected for a given current vector and a calculated leeway
angle. Leeway is calculated from the angle of heel and the speed.
The CASIO has standard 26 predefined variables A-Z which are used as
follows: The first column contains the variable, the second column the
external name and the 3rd column an explanation
____________________________________________________________________________
_____________________
Var Disp Contains
____________________________________________________________________________
_____________________
A Latitude point of departure
B Longitude point of departure
C Latitude destination
D Longitude of destination
I
J
K Course
R SR Direction of current
S SS Speed of current
V V Speed of the ship
U HH Angle of heel +/-
X
Y
Z Leeway angle
____________________________________________________________________________
_____________________
When the program is started the following information is started/displayed:
____________________________________________________________________________
_____________________
Display/Accept Description
____________________________________________________________________________
_____________________
SR Current direction
SS Current speed
V Speed of ship relative to water
HH Angle of heel (negative to larboard, positive to starboard)
VH= m.m Distance in NM
BV= k.k Speed relative to the ground
DR= g Leeway
WK= True course to steer (to be corrected for variation and
deviation on a magnetic compass)
____________________________________________________________________________
_____________________
N.B. the leeway angle is approximated with a formula from Sparkman & Stephens.
Most BASIC instructions are self explanatory but the following:
____________________________________________________________________________
_____________________
Function Description
____________________________________________________________________________
_____________________
N.B. all trig functions work by default in degrees
SET F1 Set display to 1 decimal position
PRC (a),(b) Translates the Polar coordinates with length (a) and direction
(b) into Rectangular Coordinates, stored in X and Y.
The given expression is equivalent to:
X=(a)COS(b)
Y=(a)SIN(b)
RPC (a),(b) Translates the Rectangular coordinates (a) and (b) into Polar
Coordinates, stored in X (length) and Y(direction).
The given expression is equivalent to:
X=SQRT((a)^2+(b)^2) [plain old Pythagoras]
Y=ATAN2((a),(b))
DMS a display a in format dd.mmss, i.e. DMS 3.5 displays 3° 30' 00.00"
pi the famous constant, 3.14159265358979 etc.
____________________________________________________________________________
_____________________
The other program is slightly more complex. It calculates not only the
distance and course according to the loxodrome, but also distance and
initial heading (since the course changes all the time) of a great circle
course, The coordinates of the vertex (the highest latitude reached, this is
the reason why transatlantic flights fly over the north pole) and many more
calculations:
____________________________________________________________________________
_____________________
Courses 628
____________________________________________________________________________
_____________________
1 PRC 1,D-B Normalize difference in longitude
RPC X,Y
K=Y Remember sign of difference in longitude
Z=180/pi*LN ((TAN C+1/COS C)/(TAN A+1/COS A)) meridional parts on sphere
2 RPC Z,Y
X=X/Z*(C-A) distance
PRT "LX"; print "LX" to indentify loxodrome
GSB 100 print distance and course
3 X=ACS (COS (D-B)*COS C*COS A+SIN A*SIN C) Great Circle distance
Y=ACS ((SIN C-SIN A*COS X)/COS A/SIN X Initial heading
IF SIN (D-B)<0; normalize initial heading
Y=360-Y
5 PRT "GC"; print GC to identify Great Circle
GSB 100 print distance and course
Z=ATN (1/(TAN Y*SIN A)) difference in longitude with vertex
I=ATN (TAN A/COS Z) latitude of vertex
6 PRC 1,B+Z normalize longitude vertex
RPC X,Y
J=Y
PRT "VTX: LAT="; print to identify vertex
DMS I print latitude vertex
PRT "LON=";
DMS J print longitude vertex
7 INP "GC CALC: LATX",R input latitude X
IF ABS R>=ABS I; if this is greater than the latitude of the vertex
PRT "NO LONX" there are no points on the great circle
GOTO 10 skip remainder of calculation
8 S=ACS (TAN R/TAN I) difference in longitude with vertex
PRC 1,J-S normalize western point on GC
RPC X,Y
PRT "LONX1="; print first longitude X
DMS Y
9 PRC 1,J+S normalize eastern point on GC
RPC X,Y
PRT "LONX2="; print second longitude X
DMS Y
10 INP "LONY",U input longitude Y
T=ATN (COS (J-U)*TAN I) compute latitude Y
PRT "LATY="; print latitude Y
DMS T
SET F0 set display to 0 decimals
11 M=ACS (SIN (J-U)*SIN I) course in point LatY, LonY
N=ASN (SIN M*SIN (R-T)) minimum distance from point LatX, LonY to Great Circle
IF K<0; direction along Great Circle
M=M+180 course is opposite
12 O=ATN (1/(COS (R-T)*TAN M)) bearing from point LatX, LonY to nearest
point on Great Circle
IF N>0 N of the Great Circle
O=180-O bearing is complement
GOTO 14
13 O=-O S of the Great Circle
IF O<0; normalize bearing
O=O+360
14 PRT "COURSE(YY)=";M,"BEARING(XY)=";O print course in YY and bearing in XY
SET F1 set display to 1 decimal
PRT "MIN(XY)=";N*60,"LATM="; print minimum distance in NM and coordinates
of this point of the great circle
15 PRC ABS N,O sail this distance in the direction of the bearing
DMS R+X print the latitude of LatM
PRC 1,U+Y/COS (R+X/2) normalize the longitude of LonM
RPC X,Y
PRT "LONM="; print LonM
DMS Y
GOTO 7 return for more calculations
100 IF Y<0; subroutine to normalize a given course
Y=Y+360
101 SET F1 set the display to one decimal
PRT ": DIST=";X*60;" COURSE=";###;Y print the distance in NM and Course
RET
____________________________________________________________________________
_____________________
This program calculates a loxodrome course, a great circle course and
several things around great circle courses Between the point of departure
(latitude in A, longitude in B) and the destination (latitude in C,
longitude in D), course and distance according to the loxodrome, distance
and initial heading of the Great Circle and the Latitude and Longitude of
the Vertex (point with the highest latitude on this great circle) are
calculated.
With a given LatX the corresponding LonX1 and LonX2 are calculated (for a
given latitude there are either 2, 1 or no points on the great circle).
With the given LonY the corresponding LatY is calculated as well as the
course of the ship in point LatY, LonY.
>From a point with given LatX, LonY, the minimum distance to the great circle
is calculated as well as the (radio/radar/visual) bearing to the ship. The
coordinates LatM, LonM of the point on the great circle with minimum
distance to the given point LatX, LonY are calculated as well.
____________________________________________________________________________
_____________________
Variable Display Contains
____________________________________________________________________________
_____________________
A Latitude point of departure (N/S=+/-)
B Longitude point of departure (E/W=+/-)
C Latitude destination (N/S=+/-)
D Longitude destination (E/W=+/-)
I Latitude vertex (N/S=+/-)
J Longitude vertex (E/W=+/-)
K Difference in Longitude (E/W=+/-)
M Course in LatY, LonY
N Minimum distance from LatX, LonY to great circle
O Bearing from LatX, LonY to point with minimum distance on great circle
R LatX (N/S=+/-)
S Diff Longitude LonX (E/W=+/-)
T LatY (N/S=+/-)
U LonY (E/W=+/-)
X
Y
Z
____________________________________________________________________________
_____________________
The conversation goes as follows:
____________________________________________________________________________
_____________________
Display/Accept Description
____________________________________________________________________________
_____________________
LX: DIST= m.m COURSE= g Distance and Course loxodrome
GC: DIST= m.m COURSE= g Distance and init. Course great circle
VTX: LAT= gg mm ss.s Latitude Vertex (N/S=+/-)
LON= gg mm ss.s Longitude Vertex (E/W=+/-)
_____________
Is repeated:
_____________
GC CALC: LATX Input LatX (N/S=+/-)
NO LONX Latitude is too high
LONX1= gg mm ss.s Longitude of first point on great circle
LONX2= gg mm ss.s Longitude of 2nd point on great circle
LONY Input LonY (E/W=+/-)
LATY= Display corresponding LatY
COURSE(YY)= g Course in point LatY, LonY
BEARING(XY)= g Bearing from LatX, LonY to point on GC with
min dist.
MIN(XY)= m.m Minimal distance great circle to point LatX, LonY
LATM= gg mm ss.s Latitude point on GC with min dist to LatX, LonY
LONM= gg mm ss.s Longitude point on GC with min dist to LatX, LonY
____________________________________________________________________________
_____________________
Example:
San Francisco (37° 25'N, 122° 30'W) to Yokohama (35° 30'N, 139° 40'E)
Loxodrome 4722.1 NM, course 269°
Great Circle 4479.0 NM, initial course 303°
Vertex 48° 22'5N,
169° 40.2'W
For LatX, LonY we enter the coordinates of Dutch Harbour (53° N, 166° W).
This gives the following results:
· There are no points with the given latitude LatX=53° N on the great circle
(53° has a higher latitude than the vertex), normally there are two points
on the great circle LonX1 and LonX2.
· For the given longitude LonY=166° W the latitude LatY=48° 19'N
· The course in point LatY, LonY is 273°
· The bearing from Dutch Harbour to the nearest point on the great circle is
183°
· The minimum distance the ship on the great circle passes Dutch Harbour is
280.7 NM
· The point of minimum distance on the great circle to Dutch Harbour is
LatM=48° 19.6'N, LonM=166° 21.3'W.
N.B. Both programs work accurate to any distance as the loxodrome is
calculated based on the meridional parts formula (Dutch: Vergrotende Breedte).
According to this last program, the loxodrome distance between your given
two points is 13,6 NM and a course of 284 true. According to the great
circle, the distance is the same (the distance is very short) and the
initial heading is 284 too.
N.B. If you are interested I can send the Composite track program as well.
Such tracks are used for long courses between relatively high northern or
southern latitudes , for instance from cape town to melbourne, from which
the great circle would lead to unacceptable high latitudes (through the
polar region).
>
>THANKS!
You are welcome,
Jan van Puffelen <puffelej@XXX.XXX>
52d 24.5N
4d 55.0E
>
>
>
>*** Jippie Kajee KaJoo ***
>
>
>
|