// 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

// 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)]

// 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