navigation.frink

View or download 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"


View or download 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 17419 days, 16 hours, 42 minutes ago.