Download or view navigation.frink in plain text format
// This file contains high-accuracy navigation calculations for the earth.
// Sign conventions
// It is *highly* recommended that you use these constants instead of
// assuming a sign convention.
North := +1
South := -1
West := -1
East := +1
// Turn a latitude measurement into the appropriate letter or word. (e.g.
// "N" or "S"
// Pass in the latitude and an optional flag to produce the letter
// or word. Default is letter only.
latitudeName[lat, short=true] :=
{
if lat >= 0 degrees
short ? "N" : "North"
else
short ? "S" : "South"
}
// Turn a longitude measurement into the appropriate letter or word. (e.g.
// "N" or "S"
// Pass in the longitude and an optional flag to produce the letter
// or word. Default is letter only.
longitudeName[long, short=true] :=
{
if long <= 0 degrees
short ? "W" : "West"
else
short ? "E" : "East"
}
// Calculates the initial bearing for getting from point 1
// to point 2. This corrects for the ellipsoidal shape of the earth to very
// high accuracy.
// This is based on a paper by T. Vincenty:
// "Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application
// of Nested Equations", Survey Review XXII, 176, April 1975.
// http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
// This implements the "inverse formula."
// Returns:
// [distance, initialBearing, finalBearing]
//
// TODO: Allow this to pass in a datum.
//
// TODO: Implement the version described in:
// http://link.springer.com/article/10.1007%2Fs00190-012-0578-z
// Algorithms for geodesics
/// Charles F. F. Karney // Journal of Geodesy
/// January 2013, Volume 87, Issue 1, pp 43-55,
earthDistanceAndBearing[lat1, long1, lat2, long2] :=
{
L = long2 - long1
if L < -180 degrees
L = L + 360 degrees
if L > 180 degrees
L = L - 360 degrees
lambda = L
if (lat1 == lat2) and (long1 == long2)
return [0 m, 0 degrees, 0 degrees] // Points are identical.
a = earthradius_equatorial
b = earthradius_polar
f = earth_flattening
// Calculate "reduced" latitudes
U1 = arctan[(1-f) tan[lat1]]
U2 = arctan[(1-f) tan[lat2]]
cU1 = cos[U1]
cU2 = cos[U2]
sU1 = sin[U1]
sU2 = sin[U2]
var slambda
var clambda
var oldlambda
var calpha2
var cos2sigmam
var sinSigma
var sinalpha
do
{
oldlambda = lambda
slambda = sin[lambda]
clambda = cos[lambda]
sinSigma = sqrt[(cU2 slambda)^2 + (cU1 sU2 - sU1 cU2 clambda)^2] // Eq.14
cossigma = sU1 sU2 + cU1 cU2 clambda // Eq. 15
tansigma = sinSigma / cossigma // Eq. 16
sinalpha = cU1 cU2 slambda / sinSigma // Eq. 17
calpha2 = 1 - sinalpha^2
if (calpha2 == 0) // Equatorial points
cos2sigmam = 0
else
cos2sigmam = cossigma - 2 sU1 sU2 / calpha2 // Eq. 18
C = f/16 calpha2 (4 + f (4-3 calpha2)) // Eq. 10
sigma = arctan[sinSigma, cossigma]
lambda = L + (1-C) f sinalpha (sigma +
C sinSigma ( cos2sigmam + C cossigma (-1 + 2 cos2sigmam^2)))
// Eq. 11
} while (abs[oldlambda - lambda] > 1e-6 arcsec and abs[lambda] < pi radians)
slambda = sin[lambda]
clambda = cos[lambda]
alpha1 = arctan[cU2 slambda, cU1 sU2 - sU1 cU2 clambda]
alpha2 = arctan[cU1 slambda, -sU1 cU2 + cU1 sU2 clambda]
u2 = calpha2 (a^2-b^2)/b^2
A = 1 + u2/16384 (4096 + u2 (-768 + u2 (320 - 175 u2))) // Eq. 3
B = u2/1024 (256 + u2 ( -128 + u2 (74 - 47 u2))) // Eq. 4
deltaSigma = B sinSigma ( cos2sigmam +
1/4 B (cos[sigma](-1 + 2 cos2sigmam^2) -
1/6 B cos2sigmam (-3 + 4 sinSigma^2)(-3+4 cos2sigmam^2)))
// Eq. 6
dist = b A (sigma - deltaSigma)
return [dist, alpha1 mod circle, alpha2 mod circle]
}
// Convenience method to just return the distance.
earthDistance[lat1, long1, lat2, long2] := earthDistanceAndBearing[lat1,long1,lat2,long2]@0
// Convenience method to just return the bearing.
earthBearing[lat1, long1, lat2, long2] := earthDistanceAndBearing[lat1,long1,lat2,long2]@1
// Given the lat/long of starting point, and traveling a specified distance,
// at an initial bearing, calculates the lat/long of the resulting location.
// This corrects for the ellipsoidal shape of the earth to very high accuracy.
// This is based on a paper by T. Vincenty:
// "Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application
// of Nested Equations", Survey Review XXII, 176, April 1975.
// http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
// This implements the "direct formula."
//
// Returns:
// [lat, long]
resultantLatLong[lat1, lon1, dist, bearing] :=
{
f = earth_flattening
// Calculate "reduced" latitude
U1 = arctan[(1-f) tan[lat1]]
cU1 = cos[U1]
sU1 = sin[U1]
a = earthradius_equatorial
b = earthradius_polar
cosalpha1 = cos[bearing]
sinalpha1 = sin[bearing]
sigma1 = arctan[tan[U1],cosalpha1] // Eq. 1
sinalpha = cU1 sinalpha1 // Eq. 2
calpha2 = 1 - sinalpha^2
u2 = calpha2 (a^2-b^2)/b^2
A = 1 + u2/16384 (4096 + u2 (-768 + u2 (320 - 175 u2))) // Eq. 3
B = u2/1024 (256 + u2 ( -128 + u2 (74 - 47 u2))) // Eq. 4
baseS = dist/(b A)
sigma = baseS
do
{
lastsigma = sigma
twoSigmam = 2 sigma1 + sigma // Eq. 5
cos2sigmam = cos[twoSigmam]
sinSigma = sin[sigma]
deltaSigma = B sinSigma ( cos2sigmam +
1/4 B (cos[sigma](-1 + 2 cos2sigmam^2) -
1/6 B cos2sigmam (-3 + 4 sinSigma^2)(-3+4 cos2sigmam^2)))
// Eq. 6
sigma = baseS + deltaSigma // Eq. 7
} while (abs[lastsigma - sigma] > 1e-6 arcsec)
sinSigma = sin[sigma]
cosSigma = cos[sigma]
twoSigmam = 2 sigma1 + sigma // Eq. 5
cos2sigmam = cos[twoSigmam]
lat2 = arctan[sU1 cosSigma + cU1 sinSigma cosalpha1,
(1-f)sqrt[sinalpha^2 + (sU1 sinSigma - cU1 cosSigma cosalpha1)^2]]
// Eq. 8
lambda = arctan[sinSigma sinalpha1, cU1 cosSigma - sU1 sinSigma cosalpha1]
//Eq.9
C = f/16 calpha2 (4 + f (4-3 calpha2)) // Eq. 10
L = lambda - (1-C) f sinalpha (sigma +
C sinSigma ( cos2sigmam + C cosSigma (-1 + 2 cos2sigmam^2)))
// TODO: Calculate final azimuth? (eq. 12)
return [lat2, lon1+L]
}
// Calculates the perimeter of a polygon on the earth. This uses great
// circle distances. The input value is an array containing
// 3 or more pairs of [lat, long] values.
earthPerimeter[polygon] :=
{
perimeter = 0 m
size = length[polygon]
for i = 0 to size-1
{
[lat1, long1] = polygon@i
[lat2, long2] = polygon@((i+1) mod size)
perimeter = perimeter + earthDistance[lat1, long1, lat2, long2]
}
return perimeter
}
// Calculates the area of a polygon on the earth. This uses great
// circle distances. The input value is an array containing
// 3 or more pairs of [lat, long] values.
earthArea[polygon, radius=earthradius] :=
{
anglesum = 0 radians
size = length[polygon]
for i = 0 to size-1
{
[lat1, long1] = polygon@i
[lat2, long2] = polygon@((i+1) mod size)
[lat3, long3] = polygon@((i+2) mod size)
angle = abs[earthBearing[lat2, long2, lat3, long3] -
earthBearing[lat2, long2, lat1, long1]]
if (angle > 180 degrees)
angle = (360 degrees) - angle
anglesum = anglesum + angle
}
return (anglesum/radians - (size-2) pi) radius^2
}
// Distance between two points on the earth. This is from the high-accuracy
// system of equations in Meeus, chapter 11. These equations are not numbered,
// but begin *after* equation 11.2.
// These correct for the non-sphericity of the earth.
// This function is deprecated and uses the wrong sign convention.
lowAccuracyEarthDistance[lat1, long1, lat2, long2] :=
{
F = (lat1 + lat2) / 2
//println["F: " + (F -> degrees)]
G = (lat1 - lat2) / 2
//println["G: " + (G -> degrees)]
lambda = (long1-long2) / 2
//println["lambda: " + (lambda -> degrees)]
S = sin[G]^2 cos[lambda]^2 + cos[F]^2 sin[lambda]^2
C = cos[G]^2 cos[lambda]^2 + sin[F]^2 sin[lambda]^2
//println["S: $S"]
//println["C: $C"]
omega = arctan[sqrt[S/C]]
//println["omega: " + (omega -> degrees)]
R = sqrt[S*C] / omega
//println["R: $R"]
D = 2 omega earthradius_equatorial
//println["D: $D"]
H1 = (3 R - 1)/(2 C)
H2 = (3 R + 1)/(2 S)
//earth_flattening = 1/298.257222
s = D (1 + earth_flattening (H1 sin[F]^2 cos[G]^2 -
H2 cos[F]^2 sin[G]^2))
return s
}
// Calculates the initial great circle bearing for getting from point 1
// to point 2. This is only uses spherical (and not ellipsoidal) geometry and
// does *not* correct for the fact that the earth is not quite a perfect
// sphere. This uses the wrong sign convention for longitude. This function
// is deprecated.
lowAccuracyEarthBearing[lat1, long1, lat2, long2] :=
{
arctan[sin[long1-long2]*cos[lat2],
cos[lat1]*sin[lat2]-sin[lat1]*cos[lat2]*cos[long1-long2]] mod circle
}
// Given the lat/long of starting point, and traveling a specified distance,
// at an initial bearing, calculates the lat/long of the resulting location.
// Equation taken from:
// http://williams.best.vwh.net/avform.htm#LL
// This does *not* correct for the fact that the earth is not quite a
// perfect sphere.
// This method is deprecated.
lowAccuracyResultantLatLong[lat1, lon1, dist, bearing, radius=earthradius] :=
{
d = dist/radius // Convert distance to radians
lat =arcsin[sin[lat1]*cos[d]+cos[lat1]*sin[d]*cos[bearing]]
dlon=arctan[sin[bearing]*sin[d]*cos[lat1], cos[d]-sin[lat1]*sin[lat]]
lon= ((lon1-dlon + pi) mod (2*pi))-pi
return [lat, lon]
}
"Ok"
Download or view navigation.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 20203 days, 10 hours, 32 minutes ago.