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 17476 days, 20 hours, 8 minutes ago.