WaBis

walter.bislins.ch

WGS84 JavaScript Module

This is the source code of the WGS84 JavaScript module, used in some of my Calculators and Apps. The source of the constants and equations used is given as the URL's in the code. Functions without source links are trivial and derived by me.

The World Geodetic System 1984 (WGS84) is a standard used in cartography, geodesy, and satellite navigation including GPS. The WGS 84 defines an Earth-centered, Earth-fixed (ECEF) coordinate system and a geodetic datum, and also describes the associated Earth Gravitational Model (EGM) and World Magnetic Model (WMM). The standard is published and maintained by the United States National Geospatial-Intelligence Agency.

/////////////////////////////////////////
// WGS84 Model
//
// Dependencies:
//
// * x.js: some basic browser independent DOM functions
// * JsgVectMat3.inc: V3 = JsgVect3: 3D vector functions
// * NewtonSolver.js
//
// Abbreviations:
//
// lat   = latitude in degrees
// lng   = longitude in degrees
// h     = height above surface
// azim  = azimuth angle in degrees
// elev  = elevation angle in degrees
// dir   = direction vector
// accel = acceleration
// dist  = distance

// Some Alias Definitions and Utility Functions

var V3 = JsgVect3;
var Pi    = Math.PI;     // 180°
var Pi2   = Math.PI/2;   //  90°
var sqrt  = Math.sqrt;
var sin   = Math.sin;
var asin  = Math.asin;
var cos   = Math.cos;
var acos  = Math.acos;
var tan   = Math.tan;
var atan  = Math.atan;
var atan2 = Math.atan2;
var abs   = Math.abs;
function rad(x) { return x * Math.PI / 180; }
function deg(x) { return x * 180 / Math.PI; }
function sqr(x) { return x * x; }
function limit1(x) { return x < -1 ? -1 : x > 1 ? 1 : x; }
V3.Add3 = function( v1, v2, v3 ) { return V3.Add( V3.Add( v1, v2 ), v3 ); }

// The WGS84 model -----------------------------------------

var WGS84 = {

  error: '',  // contains error from NewtonSolver if EcefToGeo can not be calculated
  nLoops: 0,

  precision: 1e-13,
  maxLoops: 100,

  // constants ----------------------------

  // GM = G * Mass of earth including mass of atmosphere:
  GM: 3.986004418e14,  // m^3/s^2, http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf

  // gravitational constant
  G: 6.67408e-11,

  // mass of the earth. Note WGS84 uses a slightly different value because they use G = 6.673e-11 m^3/kg/s^2
  // M WGS84 = 5.9733328e24 kg -> http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
  M: 3.986004418e14 / 6.67408e-11, // GM / G

  // sidereal earth rotation period
  T: 86164.098903691,  // s,     https://en.wikipedia.org/wiki/Earth

  // semi-major axis
  a: 6378137,          // m,     http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf

  // semi-minor axis
  b: 6356752.314245,   // m,     http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf

  // mean earth radius
  R: 6371008.7714,     // m,     http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf

  // gravitational acceleration at equator
  ge: 9.7803253359,    // m/s^2, http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf

  // gravitational acceleration at poles
  gp: 9.8321849378,    // m/s^2, http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf

  // flattening of the ellipsoid
  f: 1 / 298.257223563,          // https://en.wikipedia.org/wiki/Geodetic_datum

  // Reciprocal of flattening
  f_inv: 298.257223563,          // https://en.wikipedia.org/wiki/Geodetic_datum

  // First eccentricity squared
  e_sqr: 6.69437999014e-3,       // https://en.wikipedia.org/wiki/Geodetic_datum

  // Second eccentricity squared
  e2_sqr: 6.73949674228e-3,      // https://en.wikipedia.org/wiki/Geodetic_datum

  // Earth rotation rate (sidereal day) = 2 pi / T
  // Note: WGS84 uses a less accurate value
  // Omega WGS84 = 7.292115.0e-5 rad/s -> http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
  Omega: 7.2921151467069e-5,
  OmegaVect: [ 0, 0, 7.2921151467069e-5 ],

  // functions ----------------------------

  GeoToEcef: function( lat, lng, h ) {
    // or GeoToEcef( { lat, lng, h } )
    //
    // Converts lat, lng, h to [ x, y, z ]
    //
    // https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#Coordinate_system_conversion

    if (xObj(lat)) {
      // lat == { lat, lng, h }
      return this.GeoToEcef( lat.lat, lat.lng, lat.h );
    }

    var lat = rad( lat );
    var lng = rad( lng );
    var N = this.a / sqrt( 1 - this.e_sqr * sqr(sin(lat)) );
    var x = (N + h) * cos(lat) * cos(lng);
    var y = (N + h) * cos(lat) * sin(lng);
    var z = (N * sqr( this.b / this.a ) + h) * sin(lat);

    return [ x, y, z ];

  },

  EcefToGeo: function( pos, lngDefault ) {
    // Calculates latitude, longitude and height of pos.
    //
    // If pos is on the axis of rotation, longitude is not defined and can be set by lngDefault.
    // lngDefault is set to 0 if not defined.
    //
    // return { lat, lng, h }
    //
    // If algorithmus cannot find a solution, this.error is set
    // and returns { lat: NaN, lng: NaN, h: NaN }
    // All angles are in degrees.
    //
    // https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_geodetic_coordinates

    this.error = '';
    lngDefault = xDefNum( lngDefault, 0 );
    var x = pos[0];
    var y = pos[1];
    var z = pos[2];

    // compute longitude
    var lat, lng;
    var vXY = [ x, y, 0 ];
    if (V3.Length(vXY) == 0) {
      // lng is undefined -> set lng = 0
      lng = rad( lngDefault );
    } else {
      // assert V3.Length(vXY) > 0, so Norm returns no null vector
      lng = acos( limit1( V3.ScalarProd( [1,0,0], V3.Norm(vXY) ) ) );
      if (vXY[1] < 0) lng *= -1;
    }

    // compute latitude
    var p_sqr = sqr(x) + sqr(y);
    var p = sqrt( p_sqr );
    if (p < this.precision) {

      lat = z < 0 ? -Pi2 : Pi2;

    } else if (abs(z) < this.precision) {

      lat = 0;

    } else {

      // assert p != 0 and z != 0
      // store some constants in solverData because this.kFunc() has no access to 'this'
      var solverData = { MaxLoops: this.maxLoops, p_sqr: p_sqr, z: z };
      var k0 = 1 / (1 - this.e_sqr);
      var k = SolveWithNewton( this.kFunc, 0, k0, this.precision, solverData );
      if (solverData.Status) {
        this.error = solverData.Status;
        this.nLoops = solverData.NLoops;
        return { lat: NaN, lng: NaN, h: NaN };
      }
      lat = atan( k * z / p );

    }

    // compute height
    var h;
    var cos_lat = cos(lat);
    if (abs(cos_lat) < this.precision) {
      h = abs(z) - this.b;
    } else {
      var N = this.a / sqrt( 1 - this.e_sqr * sqr(sin(lat)) );
      h = p / cos_lat - N;
    }

    return { lat: deg(lat), lng: deg(lng), h: h };

  },

  kFunc: function( k, data ) {
    // Private function used in NewtonSolver of EcefToGeo()
    // https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_geodetic_coordinates

    var nom = WGS84.e_sqr * WGS84.a * k;
    var denom = sqrt( data.p_sqr + (1 - WGS84.e_sqr) * sqr(data.z * k) );

    return k - 1 - nom / denom;

  },

  Perpendicular: function( lat, lng ) {
    // or Perpendicular( { lat, lng, [h] } )
    // or Perpendicular( pos: array[3] )
    //
    // Returns a unit vector perpendicular to the surface through position at lat, lng.
    // Note: h is not needed and ignored
    //
    // Sets this.error if pos can not be converted into Geo coordinates
    // and returns [ NaN, NaN, NaN ] in this case

    // handle alternative argument formats
    if (xObj(lat)) {
      // lat == { lat, lng, [h] }
      return this.Perpendicular( lat.lat, lat.lng );
    }
    if (xArray(lat)) {
      // lat == pos
      var geo = this.EcefToGeo( lat );
      if (this.error) { return [ NaN, NaN, NaN ]; }
      lat = geo.lat;
      lng = geo.lng;
    }

    lat = rad(lat);
    lng = rad(lng);
    var r = cos(lat);
    var x = r * cos(lng);
    var y = r * sin(lng);
    var z = sin(lat);
    // assert V3.Length( [x,y,z] ) == 1

    return [ x, y, z ];

  },

  SurfacePoint: function( lat, lng, hDefault ) {
    // or SurfacePoint( { lat, lng, [h] }, hDefault )
    // or SurfacePoint( pos, hDefault )
    //
    // hDefault is optional offset of returned point above reference ellipsoid, default = 0.
    //
    // Returns a point with same lat/lng as pos but at height hDefault.
    // Returns [ NaN, NaN, NaN ] and sets this.error if pos can not be transformed to Geo coordinate system.

    if (xObj(lat)) {
      // lat == { lat, lng, [h] }, lng == hDefault
      return this.SurfacePoint( lat.lat, lat.lng, lng );
    }
    if (xArray(lat)) {
      // lat == pos, lng == hDefault
      var geo = this.EcefToGeo( lat );
      if (this.error) { return [ NaN, NaN, NaN ]; }
      hDefault = lng;
      lat = geo.lat;
      lng = geo.lng;
    }
    hDefault = xDefNum( hDefault, 0 );

    return this.GeoToEcef( lat, lng, hDefault );

  },

  GeocentricRadius: function( lat ) {
    // or GeocentricRadius( { lat, [lng, h] } )
    // or GeocentricRadius( pos: array[3] )
    //
    // The distance from the Earth's center to a point on the ellipsoid surface at latitude lat.
    // lat in degrees.
    // https://en.wikipedia.org/wiki/Earth_radius#Geocentric_radius
    //
    // Sets this.error if pos can not be converted into Geo coordinates
    // and returns NaN in this case

    // handle alternative argument formats
    if (xObj(lat)) {
      // lat == { lat, [lng, h] }
      return this.GeocentricRadius( lat.lat );
    }
    if (xArray(lat)) {
      // lat == pos
      var geo = this.EcefToGeo( lat );
      if (this.error) { return NaN; }
      lat = geo.lat;
    }

    lat = rad( lat );
    var cos_lat = cos( lat );
    var sin_lat = sin( lat );
    var A = sqr( sqr( this.a ) * cos_lat ) + sqr( sqr( this.b ) * sin_lat );
    var B = sqr( this.a * cos_lat ) + sqr( this.b * sin_lat );

    return sqrt( A / B );

  },

  MeridialRadius: function( lat ) {
    // or MeridialRadius( { lat, [lng, h] } )
    // or MeridialRadius( pos: array[3] )
    //
    // Returns the Earth's radius of curvature in the (north–south) meridian at lat.
    // lat in degrees.
    // https://en.wikipedia.org/wiki/Earth_radius#Meridional
    //
    // Sets this.error if pos can not be converted into Geo coordinates
    // and returns NaN in this case

    // handle alternative argument formats
    if (xObj(lat)) {
      // lat == { lat, [lng, h] }
      return this.MeridialRadius( lat.lat );
    }
    if (xArray(lat)) {
      // lat == pos
      var geo = this.EcefToGeo( lat );
      if (this.error) { return NaN; }
      lat = geo.lat;
    }

    lat = rad( lat );
    var M = sqr( this.a * this.b ) / Math.pow( sqr( this.a * cos(lat) ) + sqr( this.b * sin(lat) ), 1.5 );

    return M;

  },

  PrimeVerticalRadius: function( lat ) {
    // or MeridialRadius( { lat, [lng, h] } )
    // or MeridialRadius( pos: array[3] )
    //
    // Returns the prime vertical curvature which is perpendicular to MeridialRadius()
    // lat in degrees.
    // https://en.wikipedia.org/wiki/Earth_radius#Prime_vertical
    //
    // Sets this.error if pos can not be converted into Geo coordinates
    // and returns NaN in this case

    // handle alternative argument formats
    if (xObj(lat)) {
      // lat == { lat, [lng, h] }
      return this.PrimeVerticalRadius( lat.lat );
    }
    if (xArray(lat)) {
      // lat == pos
      var geo = this.EcefToGeo( lat );
      if (this.error) { return NaN; }
      lat = geo.lat;
    }

    lat = rad( lat );
    var N = sqr( this.a ) / sqrt( sqr( this.a * cos( lat ) ) + sqr( this.b * sin( lat ) ) );

    return N;

  },

  DirectionalRadius: function( azim, lat ) {
    // or DirectionalRadius( azim, { lat, [lng, h] } )
    // or DirectionalRadius( azim, pos: array[3] )
    //
    // Returns ellipsoid radius along a course azimuth (measured clockwise from north) at lat
    // Angles in degrees
    // https://en.wikipedia.org/wiki/Earth_radius#Directional
    //
    // Sets this.error if pos can not be converted into Geo coordinates
    // and returns NaN in this case

    // handle alternative argument formats
    if (xObj(lat)) {
      // lat == { lat, [lng, h] }
      return this.DirectionalRadius( azim, lat.lat );
    }
    if (xArray(lat)) {
      // lat == pos
      var geo = this.EcefToGeo( lat );
      if (this.error) { return NaN; }
      lat = geo.lat;
    }

    var a = rad( azim );
    var M = this.MeridialRadius( lat );
    var N = this.PrimeVerticalRadius( lat );
    var R = 1 / ( (sqr(cos(a)) / M) + (sqr(sin(a)) / N) );

    return R;

  },

  EnuCoordinateSystem: function( lat, lng ) {
    // or EnuCoordinateSystem( { lat, lng, [h] } )
    // or EnuCoordinateSystem( pos: array[3], lngDefault )
    //
    // pos in ECEF coordinates
    // lngDefault (optional, default = 0) is used if pos is at the poles to define
    // the direction of the ENU sysem.
    //
    // Returns local coordinate system as a matrix [ east[], north[], up[] ] at position pos.
    // All direction vectors are normalized to length 1. lngDefault in degrees.
    // Requires pos is not zero vector.
    //
    // Sets this.error if pos can not be converted into Geo coordinates
    // and returns [ [NaN,NaN,NaN], [NaN,NaN,NaN], [NaN,NaN,NaN] ] in this case

    // handle alternative argument formats
    var pos = lat;
    if (xObj(lat)) {
      // lat == { lat, lng, [h] }
      return this.EnuCoordinateSystem( lat.lat, lat.lng );
    }

    lng = xDefNum( lng, 0 );
    var up = this.Perpendicular( lat, lng ); // also handles variant (pos, lngDefault)
    if (this.error) { return [ [NaN,NaN,NaN], [NaN,NaN,NaN], [NaN,NaN,NaN] ]; }

    // if up is parallel to earth axis [0,0,1] then longitude gives the direction
    var east = V3.Mult( [0,0,1], up );
    if (V3.Length(east) < this.precision) {
      lng = rad( lng );
      if (up[2] > 0) {
        up = [0,0,1];
      } else {
        up = [0,0,-1];
      }
      east = [ -sin(lng), cos(lng), 0 ];
    } else {
      east = V3.Norm( east );
    }
    var north = V3.Mult( up, east );
    // ensures north is normed because up and east are normed perpendicular vectors

    return [ east, north, up ];

  },

  DirectionEcefToEnu: function( dir, enu ) {
    // Transforms dir from ECEF coordinate system into the
    // enu = [ east[], north[], up[] ] local coordinate system enu
    // (see EnuCoordinateSystem()).
    // The length of the returned vector is the same as dir.
    //
    // Returns [ east, north, up ] components of dir.

    var e = V3.ScalarProd( dir, enu[0] );
    var n = V3.ScalarProd( dir, enu[1] );
    var u = V3.ScalarProd( dir, enu[2] );

    return [ e, n, u ];

  },

  DirectionEnuToEcef: function( dirEnu, enu ) {
    // Transforms dirENU from enu = [ east[], north[], up[] ] coordinate system
    // into the ECEF coordinate system (see EnuCoordinateSystem()).
    // The length of the returned vector is the same as dir.
    //
    // Returns [ x, y, z ] in ECEF coordinates

    var e = V3.Scale( enu[0], dirEnu[0] );
    var n = V3.Scale( enu[1], dirEnu[1] );
    var u = V3.Scale( enu[2], dirEnu[2] );
    var dir = V3.Add3( e, n, u );

    return dir;

  },

  DirectionEnuFromAzimuthElevation: function( azim, elev ) {
    // or DirectionEnuFromAzimuthElevation( { azim, elev } )
    //
    // Calculates the direction vector in the local coordinate system ENU.
    //
    // Returns [ e, n, u ] unit vector.

    // handle alternative argument formats
    if (xObj(azim)) {
      // azim == { azim, elev }
      return this.DirectionEnuFromAzimuthElevation( azim.azim, azim.elev );
    }

    var az = rad( azim );
    var el = rad( elev );
    var r = cos( el );
    var n = r * cos( az );
    var e = r * sin( az );
    var u = sin( el );

    return [ e, n, u ];

  },

  AzimuthElevationFromDirectionEnu: function( dirEnu ) {
    // dirEnu as a vector in a local east/north/up coordinate system
    // dirEnu has not to be a unit vector, but must not be of 0 length
    //
    // returns { azim, elev } of dirEnu in degrees
    // azim = 0..360, elev = -90..90

    var dirEnu = V3.Norm( dirEnu );
    var elev = asin( dirENU[2] );
    var r = sqrt( sqr(dirENU[1]) + sqr(dirENU[0]) );
    var azim = 0;
    if (r >= this.precision) {
      azim = acos( limit1( dirENU[1] / r ) );
      if (dirENU[0] < 0) {
        azim = 2 * Pi - azim;
      }
    }

    return { azim: deg(azim), elev: deg(elev) };

  },

  DirectionFromAzimuthElevation: function( enu, azim, elev ) {
    // or DirectionFromAzimuthElevation: function( enu, { azim, elev } )
    //
    // enu = [ east[], north[], up[] ], local coordinate system in ECEF coordinates
    // azim, elev = angles in degree
    // returns normalized direction vector [ x, y, z ] ECEF coordinates
    //
    // Use this.EnuCoordinateSystem() to calculate enu.

    var dirENU = this.DirectionEnuFromAzimuthElevation( azim, elev );
    var dirECEF = this.DirectionEnuToEcef( dirENU, enu );
    // ensures dirECEF is unit vector

    return dirECEF;

  },

  DirectionFromPosAzimuthElevation: function( pos, azim, elev, lngDefault ) {
    // or DirectionFromPosAzimuthElevation( pos, { azim, elev }, lngDefault )
    //
    // pos is a location vector in ECEF coordinates.
    // azim, elev = angles in degree.
    // Returns normalized direction vector [ x, y, z ] in ECEF coordinates.
    //
    // Use GeoToEcef() to convert from lat, lng, h to ECEF coordinates:
    // this.DirectionFromPosAzimuthElevation( this.GeoToEcef(lat,lng,h), azim, elev, lng );
    //
    // lngDefault (optional, default = 0) is used if posECEF is at the poles to
    // define the direction of the ENU sysem.

    // handle alternative argument formats
    if (xObj(azim)) {
      // azim == { azim, elev }, elev = lngDefault
      return this.DirectionFromPosAzimuthElevation( pos, azim.azim, azim.elev, elev );
    }

    var enu = this.EnuCoordinateSystem( pos, lngDefault );

    return this.DirectionFromAzimuthElevation( enu, azim, elev );

  },

  AzimuthElevationFromPosDirection: function( lat, lng, dir, lngDefault ) {
    // or AzimuthElevationFromPosDirection( { lat, lng, [h] }, dir, lngDefault } )
    // or AzimuthElevationFromPosDirection( pos: array[3], dir, lngDefault )
    //
    // Calculates azimuth and elevation angles at location pos given in lat/lng or ECEF
    // for the direction dir given in ECEF by first calculating the local
    // coordinate system at pos.
    //
    // Returns { azim, elev } in degrees
    // If computation failes, this.error is set and { azim: NaN, elev: NaN } is returned.
    //
    // (lat/lng) or pos can be any location in space except the center of the earth.
    // If pos is on the poles, the longitude is lngDefault (default = 0)
    // is used to calculte the local coordinate system, from which azimuth
    // depends on.

    // handle alternative argument formats
    if (xObj(lat)) {
      // lat == { lat, lng, [h] }, lng == dir, dir == lngDefault
      return this.AzimuthElevationFromPosDirection( lat.lat, lat.lng, lng, dir );
    }
    var pos;
    if (xArray(lat)) {
      // lat == pos, lng == dir, dir == lngDefault
      lngDefault = dir;
      dir = lng;
      pos = lat;
    } else {
      pos = this.GeoToEcef( lat, lng, 0 );
    }

    // compute local coordinate system at point pos
    var enu = this.EnuCoordinateSystem( pos, lngDefault );
    if (this.error) { return { azim: NaN, elev: NaN }; }

    var dirENU = this.DirectionEcefToEnu( dir, enu );

    return this.AzimuthElevationFromDirectionEnu( dirENU );

  },

  Distance: function( lat1, lng1, lat2, lng2 ) {
    // or Distance( { lat: lat1, lng: lng1, [h: h1] }, { lat: lat2, lng: lng2, [h: h2] } )
    // or Distance( pos1, pos2 )
    //
    // Calculates the shortest distance along the surface of the reference ellipsoid
    // between 2 locations. This is approximately a great circle distance.
    // Returns NaN and sets this.error if distance can not be calculated.
    //
    // If you need the azimuth angles too, use the function DistanceAzimuthFromLocations().
    // All angles in degrees
    //
    // Note: this function uses DistanceAzimuthFromLocations() to calculate the great circle
    // distance.

    // handle alternative argument formats
    if (xObj(lat1)) {
      // lat1 == { lat1, lng1, [h1] }, lng1 == { lat2, lng2, [h2] }
      return this.Distance( lat1.lat, lat1.lng, lng1.lat, lng1.lng );
    }
    var distAzim;
    if (xArray(lat1)) {
      // lat1 == pos1, Lng1 == pos2
      distAzim = this.DistanceAzimuthFromLocations( lat1, lng1, false );
    } else {
      distAzim = this.DistanceAzimuthFromLocations( lat1, lng1, lat2, lng2, false );
    }
    // ensures this.error == '' and distAzim.dist != NaN or this.error and distAzim.dist != NaN

    return distAzim.dist;

  },

  // Vincenty's formulae to calculate distances and azimuths on the WG84 ellipsoid
  // Source: https://en.wikipedia.org/wiki/Vincenty%27s_formulae
  //
  // Notation:
  //
  // this.a          : length of semi-major axis of the ellipsoid (radius at equator)
  // this.b          : length of semi-minor axis
  // this.f          : flattening of the ellipsoid
  // lat1, lat2      : latitudes of the points
  // lng1, lng2      : longitude of the points
  // U1, U2          : reduced latitude (latitude on the auxiliary sphere)
  // L = lng2-lng1   : difference in longitude of two points
  // lambda          : Difference in longitude of the points on the auxiliary sphere
  // alpha1, alpha2  : forward azimuths at the points
  // alpha           : forward azimuth of the geodesic at the equator, if it were extended that far
  // s               : ellipsoidal distance between the two points
  // sigma           : angular separation between points
  // sigma1          : angular separation between the point and the equator
  // sigma_m         : angular separation between the midpoint of the line and the equator

  DistanceAzimuthFromLocations: function( lat1, lng1, lat2, lng2, compAzimuths ) {
    // or DistanceAzimuthFromLocations( { lat: lat1, lng: lng1, [h: h1] }, { lat: lat2, lng: lng2, [h: h2] }, compAzimuths )
    // or DistanceAzimuthFromLocations( pos1, pos2, compAzimuths )
    //
    // Calculates the shortest distance along the surface of the reference ellipsoid
    // and azimuth angles from 2 locations, using Vincenty's formula.
    // Source: https://en.wikipedia.org/wiki/Vincenty%27s_formulae
    //
    // If compAzimuths = false (default = true) then the azimuth calculation is skipped.
    // Returns { dist, azim1, azim2 }. dist is approximately a great circle distance.
    //
    // If points are antipodal the algorithm does not converge.
    // In this case this.error is set to 'Distance does not converge to a solution' and
    // { dist: NaN, azim1: NaN, azim2: NaN } is returned.
    // this.nLoops is set to the number of loops the algorithm takes.
    // All angles in degrees, azimuts 0...360

    // handle alternative argument formats
    if (xObj(lat1)) {
      // lat1 == { lat1, lng1, [h1] }, lng1 == { lat2, lng2, [h2] }, lat2 == compAzimuths
      return this.DistanceAzimuthFromLocations( lat1.lat, lat1.lng, lng1.lat, lng1.lng, lat2 );
    }
    if (xArray(lat1)) {
      // lat1 == pos1, lng1 == pos2, lat2 == compAzimuths
      compAzimuths = lat2;
      var geo2 = this.EcefToGeo( lng1 );
      if (this.error) { return { dist: NaN, azim1: NaN, azim2: NaN }; }
      lng2 = geo2.lng;
      lat2 = geo2.lat;
      var geo1 = this.EcefToGeo( lat1 );
      if (this.error) { return { dist: NaN, azim1: NaN, azim2: NaN }; }
      lng1 = geo1.lng;
      lat1 = geo1.lat;
    }

    var U1, U2, L, lambda, lambda_old, sin_U1, cos_U1, sin_U2, cos_U2, sin_lambda, cos_lambda, sin_sigma, cos_sigma, sigma, sin_alpha, sqr_cos_alpha, cos_2sigma_m, A, B, C, D, E, sqr_u, v, k1, sqr_k1, delta_sigma, s, alpha1, alpha2;

    compAzimuths = xDefBool( compAzimuths, true );
    lat1  = rad( lat1 );
    lat2  = rad( lat2 );
    lng1  = rad( lng1 );
    lng2  = rad( lng2 );
    U1    = atan( (1 - this.f) * tan( lat1 ) );
    U2    = atan( (1 - this.f) * tan( lat2 ) );
    L     = lng2 - lng1;

    sin_U1 = sin(U1);
    cos_U1 = cos(U1);
    sin_U2 = sin(U2);
    cos_U2 = cos(U2);

    lambda = L;  // init
    this.nLoops = 0;
    this.error = '';

    do {

      if (this.nLoops >= this.maxLoops) {
        // algorithm does not converge
        this.error = 'Distance does not converge to a solution';
        return { dist: NaN, azim1: NaN, azim2: NaN };
      }
      this.nLoops++;

      sin_lambda = sin(lambda);
      cos_lambda = cos(lambda);
      sin_sigma  = sqrt( sqr(cos_U2 * sin_lambda) + sqr(cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda) );
      if (sin_sigma == 0) {
        // points are identical
        return { dist: 0, azim1: 0, azim2: 0 };
      }

      cos_sigma     = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda;
      sigma         = atan2( sin_sigma, cos_sigma );
      sin_alpha     = cos_U1 * cos_U2 * sin_lambda / sin_sigma;
      sqr_cos_alpha = 1 - sqr(sin_alpha);
      if (sqr_cos_alpha != 0) {
        cos_2sigma_m  = cos_sigma - 2 * sin_U1 * sin_U2 / sqr_cos_alpha;
        C = (this.f / 16) * sqr_cos_alpha * ( 4 + this.f * (4 - 3 * sqr_cos_alpha) );
      } else {
        cos_2sigma_m = -1;
        C = 0;
      }
      D = 2 * sqr(cos_2sigma_m) - 1;  // reused below the loop
      E = sigma + C * sin(sigma) * (cos_2sigma_m + C * cos_sigma * D);

      lambda_old = lambda;
      lambda     = L + (1 - C) * this.f * sin_alpha * E;

    } while( abs(lambda_old - lambda) > this.precision );

    sqr_u  = sqr_cos_alpha * ( (sqr(this.a) - sqr(this.b)) / sqr(this.b) );
    var v  = sqrt(1 + sqr_u);
    k1     = (v - 1) / (v + 1);
    sqr_k1 = sqr(k1);

    A = (1 + 0.25 * sqr_k1) / (1 - k1);
    B = k1 * (1 - 0.375 * sqr_k1);
    C = cos_sigma * D - (B/6) * cos_2sigma_m * (4 * sqr(sin_sigma) - 3) * (4 * sqr(cos_2sigma_m) - 3);

    sin_lambda  = sin(lambda);
    cos_lambda  = cos(lambda);
    delta_sigma = B * sin_sigma * (cos_2sigma_m + 0.25 * B * C);

    s = this.b * A * (sigma - delta_sigma);

    if (compAzimuths) {

      alpha1 = atan2( cos_U2 * sin_lambda, cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda );
      alpha2 = atan2( cos_U1 * sin_lambda, -sin_U1 * cos_U2 + cos_U1 * sin_U2 * cos_lambda );

      alpha1 = deg( alpha1 );
      if (alpha1 < 0) alpha1 = 360 + alpha1;
      alpha2 = deg( alpha2 );
      if (alpha2 < 0) alpha2 = 360 + alpha2;

    } else {

      alpha1 = 0;
      alpha2 = 0;

    }

    return { dist: s, azim1: alpha1, azim2: alpha2 };

  },


  PositionAzimuthFromPositionAzimuthDistance: function( lat, lng, azim, dist ) {
    // or PositionAzimuthFromPositionAzimuthDistance( { lat, lng, [h] }, azim, dist )
    // or PositionAzimuthFromPositionAzimuthDistance( pos, azim, dist )
    //
    // Returns { lat: lat2, lng: lng2, azim: azim2 }.
    // In case of an error, this.error is set and { NaN, NaN, NaN } is returned.
    // All angles in degrees, azimuts 0...360
    //
    // Calculates the position and azimut from a geodesics starting at lat,lng in direction
    // azim a dist apart. The new point is returned as lat, lng and the direction of the
    // geodesics is returned as azim in the returned object.
    //
    // Using Vincenty's formula on the WGS84 ellipoid.
    // Source: https://www.movable-type.co.uk/scripts/latlong-vincenty.html
    // see also: https://en.wikipedia.org/wiki/Vincenty%27s_formulae

    // handle alternative argument formats
    if (xObj(lat)) {
      // lat == { lat, lng, [h] }, lng == azim, azim == dist
      return this.PositionAzimuthFromPositionAzimuthDistance( lat.lat, lat.lng, lng, azim );
    }
    if (xArray(lat)) {
      // lat == pos, lng == azim, azim == dist
      dist = azim;
      azim = lng;
      var geo = this.EcefToGeo( lat );
      lng = geo.lng;
      lat = geo.lat;
    }

    this.error  = '';
    var phi1    = rad( lat );
    var lambda1 = rad( lng );
    var alpha1  = rad( azim );

    var sin_alpha1 = sin( alpha1 );
    var cos_alpha1 = cos( alpha1 );

    var tan_U1 = (1-this.f) * tan(phi1);
    var cos_U1 = 1 / sqrt((1 + tan_U1*tan_U1));
    var sin_U1 = tan_U1 * cos_U1;

    var sigma1        = atan2( tan_U1, cos_alpha1 ); // sigma1 = angular distance on the sphere from the equator to P1
    var sin_alpha     = cos_U1 * sin_alpha1;         // alpha = azimuth of the geodesic at the equator
    var cos_sqr_alpha = 1 - sin_alpha*sin_alpha;
    var u_sqr         = cos_sqr_alpha * (sqr(this.a) - sqr(this.b)) / sqr(this.b);

    var A = 1 + u_sqr / 16384 * (4096 + u_sqr * (-768 + u_sqr * (320 - 175 * u_sqr)));
    var B = u_sqr / 1024 * (256 + u_sqr * (-128 + u_sqr * (74 - 47 * u_sqr)));

    var sin_sigma, cos_sigma, delta_sigma; // sigma = angular distance P1 P2 on the sphere
    var cos2_sigma_m; // sigma_m = angular distance on the sphere from the equator to the midpoint of the line
    var sigma_old;

    var sigma = dist / (this.b * A);
    this.nLoops = 0;

    do {

        cos2_sigma_m = cos( 2 * sigma1 + sigma );
        sin_sigma    = sin( sigma );
        cos_sigma    = cos( sigma );
        delta_sigma  = B * sin_sigma * (cos2_sigma_m + B/4 * (cos_sigma  *(-1 + 2 * cos2_sigma_m*cos2_sigma_m) -
            B/6 * cos2_sigma_m * (-3 + 4 * sin_sigma*sin_sigma) * (-3 + 4 * cos2_sigma_m*cos2_sigma_m)));
        sigma_old    = sigma;
        sigma        = dist / (this.b * A) + delta_sigma;

    } while (abs(sigma-sigma_old) > this.precision && ++this.nLoops < this.maxLoops);

    if (this.nLoops >= 100) {
      this.error = 'Vincenty formula failed to converge';
      return { lat: NaN, lng: NaN, azim: NaN };
    }

    var x      = sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_alpha1;
    var phi2   = atan2( sin_U1 * cos_sigma + cos_U1 * sin_sigma * cos_alpha1, (1-this.f) * sqrt( sin_alpha*sin_alpha + x*x ));
    var lambda = atan2( sin_sigma * sin_alpha1, cos_U1 * cos_sigma - sin_U1 * sin_sigma * cos_alpha1 );

    var C = this.f/16 * cos_sqr_alpha * (4 + this.f * (4 - 3 * cos_sqr_alpha));
    var L = lambda - (1-C) * this.f * sin_alpha * (sigma + C * sin_sigma * (cos2_sigma_m + C * cos_sigma * (-1 + 2 * cos2_sigma_m*cos2_sigma_m)));

    var lambda2 = lambda1 + L;

    var alpha2 = deg( atan2( sin_alpha, -x ) );
    if (alpha2 < 0) alpha2 = 360 + alpha2;

    return { lat: deg(phi2), lng: deg(lambda2), azim: alpha2 };

  },

  Centrifugal: function( lat, lng, h ) {
    // or Centrifugal( { lat, lng, h } )
    // or Centrifugal( pos: array[3], lngDefault )
    //
    // Calculates the centrifugal acceleration at position (lat,lng,h) or pos.
    // If pos is on the axis of rotation, longitude is not defined and can be set by lngDefault.
    // lngDefault is set to 0 if not defined.
    //
    // returns { accel: full vector, magn: magnitude, dir: normed vector }

    // handle alternative argument formats
    if (xObj(lat)) {
      // lat == { lat, lng, h }
      return this.Centrifugal( lat.lat, lat.lng, lat.h );
    }
    var pos;
    if (xArray(lat)) {
      // lat == pos, lng == lngDefault
      pos = lat;
      lng = xDefNum( lng, 0 );
    } else {
      pos = this.GeoToEcef( lat, lng, h );
    }

    var r = sqrt( sqr(pos[0]) + sqr(pos[1]) );
    if (r < this.precision) {
      lng = rad(lng);
      var dir = [ cos(lng), sin(lng), 0 ];
      return { accel: [0,0,0], magn: 0, dir: dir };
    }
    var dir = V3.Norm( [ pos[0], pos[1], 0 ] );
    var a = sqr( this.Omega ) * r;

    return { accel: V3.Scale( dir, a ), magn: a, dir: dir };

  },

  Coriolis: function( velocityEcef ) {
    // returns coriolis acceleration vector in ECEF coordinates.
    //
    return V3.Scale( V3.Mult( this.OmegaVect, velocityEcef ), -2 );
  },

  Gravity: function( lat, lng, h ) {
    // or Gravity( { lat, lng, h } )
    // or Gravity( pos: array[3] )
    //
    // Effective (normal) gravitational acceleration at pos inclusive centrifugal correction.
    // This gravity acts perpendicular to the ellipsiodal surface through pos.
    //
    // returns { accel: full vector, magn: magnitude, dir: normed vector }
    //
    // If algorithmus cannot find a solution, this.error is set and
    // { accel: [NaN,NaN,NaN], magn: NaN, dir: [NaN,NaN,NaN] } is returned
    //
    // see https://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
    //     Normal Gravity Above the Ellipsoid
    //
    // Note: this calculation is only accurate for low altitudes. A more accurate version
    // for moderate and high altitudes can be found in the link above.
    // If you want the result to an accuracy of eps (e.g. eps = 0.01 for 1% accuracy), h
    // must not be greater than h = eps * a / 2, where a is the semi-major axis.

    // handle alternative argument formats
    if (xObj(lat)) {
      // lat == { lat, lng, h }
      return this.Gravity( lat.lat, lat.lng, lat.h );
    }
    if (xArray(lat)) {
      // lat == pos
      var geo = this.EcefToGeo( lat );
      if (this.error) { return { accel: [NaN,NaN,NaN], magn: NaN, dir: [NaN,NaN,NaN] }; }
      lat = geo.lat;
      lng = geo.lng;
      h   = geo.h;
    }

    var up = this.Perpendicular( lat, lng );

    lat = rad( lat );
    lng = rad( lng );

    // normal gravity at the ellipsoidal surface
    var k = (this.b * this.gp) / (this.a * this.ge) - 1;
    var g = this.ge * ( 1 + k * sqr(sin(lat)) ) / sqrt( 1 - this.e_sqr * sqr(sin(lat)) );

    if (abs(h) >= this.precision) {

      // normal gravity above the ellipsoidal surface

      var m = sqr(this.Omega * this.a) * this.b / this.GM;
      var gh = g * ( 1 - 2/this.a * ( 1 + this.f + m - 2 * this.f * sqr(sin(lat)) ) * h + 3/sqr(this.a) * sqr(h) );
      g = gh;

    }

    return { accel: V3.Scale( up, -g ), magn: g, dir: V3.Scale( up, -1 ) };

  },

  GravityMass: function( lat, lng, h ) {
    // or GravityMass( { lat, lng, h } )
    // or GravityMass( pos: array[3] )
    //
    // Gravitational acceleration at pos due to earth's mass alone (no centrifugal correction).
    //
    // returns { accel: full vector, magn: magnitude, dir: normed vector }
    //
    // If algorithmus cannot find a solution, this.error is set and
    // { accel: [NaN,NaN,NaN], magn: NaN, dir: [NaN,NaN,NaN] } is returned
    //
    // This function calculates first the effective gravitational acceleration
    // and then subtracts the centrifugal acceleration.
    //
    // Note: this calculation is only accurate for low altitudes. A more accurate version
    // for moderate and high altitudes can be found in the link above.
    // If you want the result to an accuracy of eps (e.g. eps = 0.01 for 1% accuracy), h
    // must not be greater than h = eps * a / 2, where a is the semi-major axis.

    // handle alternative argument formats
    if (xObj(lat)) {
      // lat == { lat, lng, h }
      return this.GravityMass( lat.lat, lat.lng, lat.h );
    }
    var pos;
    if (xArray(lat)) {
      // lat == pos
      pos = lat;
    } else {
      pos = this.GeoToEcef( lat, lng, h );
    }

    var c = this.Centrifugal( pos );
    var g = this.Gravity( pos );
    if (this.error) { return { accel: [NaN,NaN,NaN], magn: NaN, dir: [NaN,NaN,NaN] }; }

    var gmVect = V3.Sub( g.accel, c.accel );
    var gm = V3.Length( gmVect );
    var dir = V3.Norm( gmVect );

    return { accel: V3.Scale( dir, gm ), magn: gm, dir: dir };

  },

  TransformToCS: function( cs, vect ) {
    // Transforms vect [x,y,z] from ECEF to the coordinate system,
    // given as a matrix of the basis vectors cs = [ e1, e2, e3 ].
    // An example for cs is the ENU coordinate system as optained
    // by the function EnuCoordinateSystem(). But other local
    // coordinat systems can be used, given we know their basis
    // vectors in ECEF coordinatas.
    //
    // returns the vector [x,y,z] in the new coordinate system.
    //
    var x = V3.ScalarProd( cs[0], vect );
    var y = V3.ScalarProd( cs[1], vect );
    var z = V3.ScalarProd( cs[2], vect );
    return [ x, y, z ];
  },

};

More Page Infos / Sitemap
Created Friday, June 9, 2023
Scroll to Top of Page
Changed Wednesday, August 16, 2023