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