The W3FB10 routine computes the bearing and great circle distance.
from point (1) to point (2) assuming a spherical earth. The
North and South poles are special cases. If latitude of point
(1) is within 1e-10 degrees of the North pole, bearing is the
negative longitude of point (2) by convention. If latitude of
point (1) is within 1e-10 degrees of the south pole, bearing is
the longitude of point (2) by convention. If point (2) is within
1e-6 radians of the antipode of point (1), the bearing will be
set to zero. If point (1) and point (2) are within 1e-10 radians
of each other, both bearing and distance will be set to zero.
USAGE: CALL W3FB10(DLAT1, DLON1, DLAT2, DLON2, BEARD, GCDKM)
Input argument list:
DLAT1 - REAL LATITUDE OF POINT (1) IN DEGREES NORTH.
DLON1 - REAL LONGITUDE OF POINT (1) IN DEGREES EAST.
DLAT2 - REAL LATITUDE OF POINT (2) IN DEGREES NORTH.
DLON2 - REAL LONGITUDE OF POINT (2) IN DEGREES EAST.
Output argument list:
BEARD - REAL BEARING OF POINT (2) FROM POINT (1) IN
COMPASS DEGREES WITH NORTH = 0.0, VALUES FROM
-180.0 TO +180.0 DEGREES.
GCDKM - REAL GREAT CIRCLE DISTANCE FROM POINT (1) TO
POINT (2) IN KILOMETERS.
According to the NMC Handbook, the earth's radius is
6371.2 kilometers. This is what we use, even though the value
recommended by the Smithsonian Meteorological Handbook is
6371.221 KM. (I wouldn't want you to think that i didn't know
what the correct value was.)
METHOD: The poles are special cases, and handled separately.
Otherwise, from spherical trigonometry, the law of cosines is used
to calculate the third side of the spherical triangle having
sides from the pole to points (1) and (2) (the colatitudes).
Then the law of sines is used to calculate the angle at point
(1). A test is applied to see whether the arcsine result may be
be used as such, giving an acute angle as the bearing, or whether
The arcsine result should be subtracted from pi, giving an obtuse
angle as the bearing. This test is derived by constructing a
right spherical triangle using the pole, point (2), and the
meridian through point(1). The latitude of the right-angled
vertex then provides a test--if latitude (1) is greater than this
latitude, the bearing angle must be obtuse, otherwise acute.
If the two points are within 1e-6 radians of each other
a flat earth is assumed, and the four-quadrant arctangent
function is used to find the bearing. The y-displacement is
the difference in latitude and the x-displacement is the
difference in longitude times cosine latitude, both in radians.
distance is then the diagonal.
Fundamental trigonometric identities are used freely, such
as that cos(x) = sin(pi/2 - x), etc. See almost any mathematical
Handbook, such as the C.R.C. standard math tables under 'relations
in any spherical triangle', or the National Bureau of Standards
'Handbook of mathematical functions' under section 4.3.149,
formulas for solution of spherical triangles.
Double precision is used internally because of the wide
range of geographic values that may be used.
||Library contains Fortran 90 decoder/encoder
routines for GRIB edition 1. (Fortran90)
Date posted: 2/22/2007