WaBis

walter.bislins.ch

Code: Globe and Flat Earth Transformations and Mappings

JavaScript source code for the page Globe and Flat Earth Transformations and Mappings.

#INCLUDE NewtonSolver.inc
#INCLUDE JsGraphX3D.inc
#INCLUDE ControlPanel.inc
#INCLUDE EarthMap.inc
#INCLUDE WGS84.inc
#INCLUDE FlatEarthMath.inc

<jscript>

var Pi = Math.PI;
var Pi2 = Math.PI/2;
function ToRad( d ) { return d * Pi / 180; }
function ToDeg( d ) { return d * 180 / Pi; }
function Limit1( x ) { return x < -1 ? -1 : x > 1 ? 1 : x; }
function Zero( x ) { return abs(x) < 0.0001 ? 0 : x; }
function sqr( x ) { return x * x; }
var abs = Math.abs;
var sin = Math.sin;
var cos = Math.cos;
var asin = Math.asin;
var acos = Math.acos;
var sqrt = Math.sqrt;
var V3 = JsgVect3;

///////////////////////////////////////////////////////////////////

var Model = {

  Lat: -20,
  Long: 47,
  Alt: 3000000,

  Xglobe: 0,
  Yglobe: 0,
  Zglobe: 0,

  Xfe: 0,
  Yfe: 0,
  Zfe: 0,

  VlenGlobe: 0,
  VlenFe: 0,

  HAng: 30,
  VAng: 25,
  GlobeZoom: 0.7,
  FEZoom: 1.0,
  Alpha: 0.25,

  ShowContinents: true,
  ShowGrid: true,
  ShowCoordSystem: true,
  ShowVector: true,
  ShowCartCoord: true,
  ShowGeoCoord: false,
  ShowData: true,

  CoordSysSize: 1.3,
  CameraDist: 1e10, // distance of virtual camera from earths center
  Scale: 2.5,
  R: 6371000,
  Rfe: 0,

  Graph: null,

  Create: function() {
    this.Rfe = Pi * this.R;

    this.Graph = NewGraphX3D( {
      Id: 'ModelCanvas',
      Width: '100%',
      Height: '50%',
      DrawFunc: function() { Model.Draw(); },
      AutoReset: true,
      AutoClear: true,
      AutoScalePix: true
    } );

    EarthMap.SetLandColor( 'Antarctica', '#ddd', '#bbb' );
    this.Update();
  },

  Update: function( field ) {

    if (ControlPanels.MatchesField( field, 'Lat,Long,Alt' )) {

      var ecef = WGS84.GeoToEcef( this.Lat, this.Long, this.Alt );
      this.Xglobe = ecef[0];
      this.Yglobe = ecef[1];
      this.Zglobe = ecef[2];
      var ecef = FlatEarthMath.GeoToEcef( this.Lat, this.Long, this.Alt );
      this.Xfe = ecef[0];
      this.Yfe = ecef[1];
      this.Zfe = ecef[2];

    } else if (ControlPanels.MatchesField( field, 'Xglobe,Yglobe,Zglobe' )) {

      var geo = WGS84.EcefToGeo( [ this.Xglobe, this.Yglobe, this.Zglobe ] );
      this.Lat = geo.lat;
      this.Long = geo.lng;
      this.Alt = geo.h;
      var ecef = FlatEarthMath.GeoToEcef( this.Lat, this.Long, this.Alt );
      this.Xfe = ecef[0];
      this.Yfe = ecef[1];
      this.Zfe = ecef[2];

    } else if (ControlPanels.MatchesField( field, 'Xfe,Yfe,Zfe' )) {

      var geo = FlatEarthMath.EcefToGeo( [ this.Xfe, this.Yfe, this.Zfe ] );
      this.Lat = geo.lat;
      this.Long = geo.lng;
      this.Alt = geo.h;
      var ecef = WGS84.GeoToEcef( this.Lat, this.Long, this.Alt );
      this.Xglobe = ecef[0];
      this.Yglobe = ecef[1];
      this.Zglobe = ecef[2];

    }

    this.VlenGlobe = V3.Length( [ this.Xglobe, this.Yglobe, this.Zglobe ] );
    this.VlenFe = V3.Length( [ this.Xfe, this.Yfe, this.Zfe ] );

    ControlPanels.Update();
    this.Draw();
  },

  Draw: function() {
    var g = this.Graph;
    g.Reset3D();
    g.SetAngleMeasure( 'deg' );

    this.DrawGlobe( g );
    this.DrawFE( g );
  },

  DrawGlobe: function( g ) {
    EarthMap.Radius = this.R;

    // calculate scene size and setup virtual 3D camera
    var alpha = acos( EarthMap.Radius / this.CameraDist );
    var scale = this.Scale * EarthMap.Radius / cos( Pi2 - alpha );
    g.SetCameraScale( scale );
    g.SetCameraView( [0,0,-2500000], this.HAng, this.VAng, this.CameraDist );
    g.SetCameraZoom( this.GlobeZoom );

    // set left viewport and map camera to this viewport
    // camera clipping clips graphic that extends behind the camera sensor plane
    g.SetViewport( 0, 0, g.CanvasWidth/2, g.CanvasHeight, false, true );
    g.SetViewportRel( 0, 0, 2, 0 );
    g.SetWindowToCameraScreen();
    g.SetCameraClipping( 0.001 );

    if (this.ShowContinents) {
      g.SetAlpha( this.Alpha );
      EarthMap.DrawGlobe( g );
      g.SetAlpha( 1 );
    } else {
      // compute clip plane so earth grid gets clipped correctly
      EarthMap.CompClipPlane( g );
    }

    if (this.ShowGrid) {
      // requires globe is drawn or clipping plane is computed, see EarthMap.CompClipPlane()
      g.SetLineAttr( 'black', 1 );
      g.SetAlpha( 0.15 );
      EarthMap.DrawGlobeGrid( g );

      g.SetAlpha( 0.5 );
      EarthMap.DrawGlobeEquator( g, 0 );
      EarthMap.DrawGlobeMeridian( g, 0 );
    }

    if (this.ShowCoordSystem) {
      this.DrawGlobeCoordSys( g, this.CoordSysSize * this.R );
    }

    if (this.ShowVector) {
      this.DrawVectorGlobe( g );
    }

    if (this.ShowCartCoord) {
      this.DrawGlobeCartCoord( g );
    }

    if (this.ShowGeoCoord) {
      this.DrawGlobeGeoCoord( g );
    }

    if (this.ShowData) {
      this.DrawGlobeData( g );
    }

  },

  DrawFE( g ) {
    EarthMap.Radius = this.Rfe;

    // calculate scene size and setup virtual 3D camera
    var alpha = acos( EarthMap.Radius / this.CameraDist );
    var scale = this.Scale * EarthMap.Radius / cos( Pi2 - alpha );
    g.SetCameraScale( scale );
    g.SetCameraView( [0,0,-2500000], this.HAng, this.VAng, this.CameraDist );
    g.SetCameraZoom( this.FEZoom );

    // set left viewport and map camera to this viewport
    // camera clipping clips graphic that extends behind the camera sensor plane
    g.SetViewport( g.CanvasWidth/2, 0, g.CanvasWidth/2, g.CanvasHeight, false, true );
    g.SetViewportRel( 2, 0, 0, 0 );
    g.SetWindowToCameraScreen();
    g.SetCameraClipping( 0.001 );

    g.SetPlane( [0,0,0], [0,1,0], [-1,0,0] );
    EarthMap.FEMode = 2;
    if (this.ShowContinents) {
      g.SetAlpha( this.Alpha );
      EarthMap.DrawFlatEarth( g );
      g.SetAlpha( 1 );
    }

    if (this.ShowGrid) {
      g.SetLineAttr( 'black', 1 );
      g.SetAlpha( 0.15 );
      EarthMap.DrawFlatEarthGrid( g );

      g.SetAlpha( 0.5 );
      EarthMap.DrawFlatEarthEquator( g, 0 );
      EarthMap.DrawFlatEarthMeridian( g, 0 );
      EarthMap.DrawFlatEarthBorder( g );
    }

    if (this.ShowCoordSystem) {
      this.DrawFECoordSys( g, this.CoordSysSize * 0.95 * this.Rfe );
    }

    if (this.ShowVector) {
      this.DrawVectorFE( g );
    }

    if (this.ShowCartCoord) {
      this.DrawFECartCoord( g );
    }

    if (this.ShowGeoCoord) {
      this.DrawFEGeoCoord( g );
    }

    if (this.ShowData) {
      this.DrawFEData( g );
    }

  },

  DrawGlobeCoordSys: function( g, size ) {
    g.SetAlpha( 1 );

    g.SetMarkerAttr( 'Arrow1', 12, 'black', 'black', 1 );
    g.Arrow3D( [-this.R,0,0], [size,0,0], 1, 3 );
    g.Arrow3D( [0,-this.R,0], [0,size,0], 1, 3 );
    g.Arrow3D( [0,0,0], [0,0, size], 1, 3 );

    var textDist = 1.1 * size;
    g.SetTextAttrS( 'Arial', 12, 'black', 'C' );
    g.Text3D( 'X', [ textDist, 0, 0 ] );
    g.Text3D( 'Y', [ 0, textDist, 0 ] );
    g.Text3D( 'Z', [ 0, 0, textDist ] );
  },

  DrawFECoordSys: function( g, size ) {
    g.SetAlpha( 1 );

    g.SetMarkerAttr( 'Arrow1', 12, 'black', 'black', 1 );
    g.Arrow3D( [0,0,0], [size,0,0], 1, 3 );
    g.Arrow3D( [0,0,0], [0,size,0], 1, 3 );
    g.Arrow3D( [0,0,0], [0,0, size/2], 1, 3 );

    var textDist = 1.1 * size;
    g.SetTextAttrS( 'Arial', 12, 'black', 'C' );
    g.Text3D( 'X', [ textDist, 0, 0 ] );
    g.Text3D( 'Y', [ 0, textDist, 0 ] );
    g.Text3D( 'Z', [ 0, 0, textDist/2 ] );
  },

  DrawVectorGlobe: function( g ) {
    g.SetAlpha( 0.5 );
    g.SetLineAttr( 'gray', 3 );
    g.Line3D( [0,0,0], [this.Xglobe,this.Yglobe,0] );
    g.SetAlpha( 1 );
    g.SetMarkerAttr( 'Arrow1', 6, 'magenta', 'magenta', 3 );
    g.Arrow3D( [0,0,0], [this.Xglobe,this.Yglobe,this.Zglobe], 1, 3 );
  },

  DrawVectorFE: function( g ) {
    g.SetAlpha( 0.5 );
    g.SetLineAttr( 'gray', 3 );
    g.Line3D( [0,0,0], [this.Xfe,this.Yfe,0] );
    g.SetAlpha( 1 );
    g.SetMarkerAttr( 'Arrow1', 6, 'magenta', 'magenta', 3 );
    g.Arrow3D( [0,0,0], [ this.Xfe, this.Yfe, this.Zfe ], 1, 3 );
  },

  DrawGlobeCartCoord: function( g ) {
    g.SetAlpha( 1 );
    g.SetLineAttr( 'red', 2 );
    g.Line3D( [ 0, 0, 0 ], [ this.Xglobe, 0, 0 ] );
    g.SetLineAttr( 'blue', 2 );
    g.Line3D( [ this.Xglobe, 0, 0 ], [ this.Xglobe, this.Yglobe, 0 ] );
    g.SetLineAttr( 'green', 2 );
    g.Line3D( [ this.Xglobe, this.Yglobe, 0 ], [ this.Xglobe, this.Yglobe, this.Zglobe ] );
  },

  DrawFECartCoord: function( g ) {
    g.SetAlpha( 1 );
    g.SetLineAttr( 'red', 2 );
    g.Line3D( [ 0, 0, 0 ], [ this.Xfe, 0, 0 ] );
    g.SetLineAttr( 'blue', 2 );
    g.Line3D( [ this.Xfe, 0, 0 ], [ this.Xfe, this.Yfe, 0 ] );
    g.SetLineAttr( 'green', 2 );
    g.Line3D( [ this.Xfe, this.Yfe, 0 ], [ this.Xfe, this.Yfe, this.Zfe ] );
  },

  DrawGlobeGeoCoord: function( g ) {
    g.SetPlane( [0,0,0], [1,0,0], [0,0,1] );

    g.SetLineAttr( 'red', 2 );
    var r = this.R;
    if (this.Lat < 0) r *= -1;
    g.ArcOnPlane( 0, 0, r, 0, this.Lat, 1 );

    var z = this.R * sin( ToRad( this.Lat ) );
    var r = this.R * cos( ToRad( this.Lat ) );
    if (this.Long < 0) r *= -1;
    g.SetPlane( [0,0,z], [1,0,0], [0,1,0] );
    g.SetLineAttr( 'blue', 2 );
    g.ArcOnPlane( 0, 0, r, 0, this.Long, 1 );

    var vect = [ this.Xglobe, this.Yglobe, this.Zglobe ];
    var len = V3.Length( vect );
    if (len > 0) {
      var dir = V3.Norm( vect );
      var surf = V3.Scale( dir, this.R );
      g.SetLineAttr( 'green', 2 );
      g.Line3D( surf, vect );
    }
  },

  DrawFEGeoCoord: function( g ) {
    eq = this.Rfe / 2;
    var rlat = eq * (1 - (this.Lat / 90));
    g.SetLineAttr( 'red', 2 );
    g.Line3D( [eq,0,0], [rlat,0,0] );

    if (this.Long < 0) rlat *= -1;
    g.SetPlane( [0,0,0], [1,0,0], [0,1,0] );
    g.SetLineAttr( 'blue', 2 );
    g.ArcOnPlane( 0, 0, rlat, 0, this.Long, 1 );

    g.SetLineAttr( 'green', 2 );
    g.Line3D( [ this.Xfe, this.Yfe, 0 ], [ this.Xfe, this.Yfe, this.Zfe ] );
  },

  DrawGlobeData: function( g ) {
    data = {
      name: 'V',
      sub: 'ecef',
      model: 'globe',
      label: [ 'X', 'Y', 'Z' ],
      value: [ this.Xglobe/1000, this.Yglobe/1000, this.Zglobe/1000 ],
      leng: this.VlenGlobe / 1000,
      digits: [ 4, 4, 4 ],
      fix: [ 0, 0, 0 ],
    };

    this.DrawData( g, data, 0.2 * g.VpWidth, this.ShowCartCoord );

    data = {
      name: 'V',
      sub: 'geo',
      model: 'globe',
      label: [ 'Lat', 'Long', 'Alt' ],
      value: [ this.Lat, this.Long, this.Alt/1000 ],
      leng: null,
      digits: [ 3, 4, 4 ],
      fix: [ 1, 1, 0 ],
    };

    this.DrawData( g, data, 0.7 * g.VpWidth, this.ShowGeoCoord );
  },

  DrawFEData: function( g ) {
    data = {
      name: 'V',
      sub: 'geo',
      model: 'fe',
      label: [ 'Lat', 'Long', 'Alt' ],
      value: [ this.Lat, this.Long, this.Alt/1000 ],
      leng: null,
      digits: [ 3, 4, 4 ],
      fix: [ 1, 1, 0 ],
    };

    this.DrawData( g, data, 0.3 * g.VpWidth, this.ShowGeoCoord );

    data = {
      name: 'V',
      sub: 'cart',
      model: 'fe',
      label: [ 'X', 'Y', 'Z' ],
      value: [ this.Xfe/1000, this.Yfe/1000, this.Zfe/1000 ],
      leng: this.VlenFe / 1000,
      digits: [ 4, 4, 4 ],
      fix: [ 0, 0, 0 ],
    };

    this.DrawData( g, data, 0.8 * g.VpWidth, this.ShowCartCoord );
  },

  DrawData: function( g, data, x, colorize ) {

    function col( color ) {
      return colorize ? color : 'black';
    }

    var oldTrans = g.SelectTrans( 'viewport' );
    g.SetTextRendering( 'html' );

    var fontSize = 16;
    var charWidth = g.ScalePix( fontSize * 0.65 );
    var cellHeight = 1.2 * g.ScalePix( fontSize );
    var cellWidthName = (data.name.length + data.sub.length) * charWidth * 0.7;
    var cellWidthLabel = Math.max( data.label[0].length, data.label[1].length, data.label[2].length ) * charWidth;
    var cellWidthEq = 3 * charWidth;
    var cellWidthValue = Math.max( data.digits[0], data.digits[1], data.digits[2] ) * charWidth;
    var totWidth = cellWidthName + 2*cellWidthEq + cellWidthLabel + cellWidthValue;
    var y = g.VpHeight - 4 * cellHeight;

    g.SetTextAttrS( 'Arial', fontSize, 'black' );

    if (data.leng != null) {
      g.SetTextColor( 'magenta' );
      g.Text( '|'+data.name+'<sub>'+data.model+'</sub>| = '+data.leng.toFixed(3), x, y - 1.5 * cellHeight );
    }

    var txtX = x - totWidth/2 + cellWidthName/2;
    var txtY = y + 1.5 * cellHeight;
    g.SetTextColor( 'black' );
    g.Text( data.name+'<sub>'+data.sub+'</sub>', txtX, txtY );

    txtX += (cellWidthName + cellWidthEq) / 2;
    g.Text( '=', txtX, txtY );

    txtX += (cellWidthEq + cellWidthLabel) / 2;
    txtY = y + cellHeight / 2;
    g.SetTextColor( col('red') );
    g.Text( data.label[0], txtX, txtY );
    txtY += cellHeight;
    g.SetTextColor( col('blue') );
    g.Text( data.label[1], txtX, txtY );
    txtY += cellHeight;
    g.SetTextColor( col('green') );
    g.Text( data.label[2], txtX, txtY );

    txtY = y + 1.5 * cellHeight;
    txtX += (cellWidthLabel + cellWidthEq) / 2;
    g.SetTextColor( 'black' );
    g.Text( '=', txtX, txtY );

    txtX += (cellWidthEq + cellWidthValue) / 2;
    txtY = y + cellHeight / 2;
    g.SetTextColor( col('red') );
    g.Text( data.value[0].toFixed(data.fix[0]), txtX, txtY );
    txtY += cellHeight;
    g.SetTextColor( col('blue') );
    g.Text( data.value[1].toFixed(data.fix[1]), txtX, txtY );
    txtY += cellHeight;
    g.SetTextColor( col('green') );
    g.Text( data.value[2].toFixed(data.fix[2]), txtX, txtY );

    g.SetTextColor( 'black' );
    g.SelectTrans( oldTrans );
  },

};

Model.Create();

///////////////////////////////////////////////////////////////////////////

ControlPanels.NewSliderPanel( {
  Name: 'VectorSliderPanel',
  ModelRef: 'Model',
  NCols: 1,
  OnModelChange: function(field){ Model.Update(field); },
  Format: 'fix0',
  Digits: 3,
  PanelFormat: 'InputSmallerWidth',

} ).AddValueSliderField( {
  Name: 'Lat',
  Color: 'red',
  Min: -90,
  Max: 90,
  LowerLimit: -90,
  UpperLimit: 90,
  SnapTo: 0,
  Units: '°',
  Inc: 10,

} ).AddValueSliderField( {
  Name: 'Long',
  Color: 'red',
  Min: -180,
  Max: 180,
  LowerLimit: -180,
  UpperLimit: 180,
  SnapTo: [ -90, 0, 90 ],
  Units: '°',
  Inc: 10,

} ).AddValueSliderField( {
  Name: 'Alt',
  Color: 'red',
  Min: -6371000,
  Max: 1e7,
  LowerLimit: -6371000,
  UpperLimit: 1e7,
  SnapTo: 0,
  Mult: 1000,
  Units: 'km',
  Inc: 10,

} ).AddValueSliderField( {
  Name: 'Xglobe',
  Color: 'blue',
  Min: -1e7,
  Max: 1e7,
  LowerLimit: -1e7,
  UpperLimit: 1e7,
  SnapTo: [ -6371000, 0, 6371000 ],
  Mult: 1000,
  Units: 'km',
  Inc: 10,

} ).AddValueSliderField( {
  Name: 'Yglobe',
  Color: 'blue',
  Min: -1e7,
  Max: 1e7,
  LowerLimit: -1e7,
  UpperLimit: 1e7,
  SnapTo: [ -6371000, 0, 6371000 ],
  Mult: 1000,
  Units: 'km',
  Inc: 10,

} ).AddValueSliderField( {
  Name: 'Zglobe',
  Color: 'blue',
  Min: -1e7,
  Max: 1e7,
  LowerLimit: -1e7,
  UpperLimit: 1e7,
  SnapTo: [ -6371000, 0, 6371000 ],
  Mult: 1000,
  Units: 'km',
  Inc: 10,

} ).AddValueSliderField( {
  Name: 'Xfe',
  Color: 'green',
  Min: -3e7,
  Max: 3e7,
  LowerLimit: -3e7,
  UpperLimit: 3e7,
  SnapTo: [ -Model.Rfe, - Model.Rfe/2, 0, Model.Rfe/2, Model.Rfe ],
  Mult: 1000,
  Units: 'km',
  Inc: 10,

} ).AddValueSliderField( {
  Name: 'Yfe',
  Color: 'green',
  Min: -3e7,
  Max: 3e7,
  LowerLimit: -3e7,
  UpperLimit: 3e7,
  SnapTo: [ -Model.Rfe, - Model.Rfe/2, 0, Model.Rfe/2, Model.Rfe ],
  Mult: 1000,
  Units: 'km',
  Inc: 10,

} ).AddValueSliderField( {
  Name: 'Zfe',
  Color: 'green',
  Min: -1e7,
  Max: 1e7,
  LowerLimit: -1e7,
  UpperLimit: 1e7,
  SnapTo: 0,
  Mult: 1000,
  Units: 'km',
  Inc: 10,

} ).Render();

ControlPanels.NewSliderPanel( {
  Name: 'ViewSliderPanel',
  ModelRef: 'Model',
  NCols: 1,
  OnModelChange: function(field){ Model.Update(field); },
  Format: 'fix0',
  Digits: 2,
  PanelFormat: 'InputSmallerWidth',

} ).AddValueSliderField( {
  Name: 'HAng',
  Color: 'green',
  Min: -360,
  Max: 360,
  LowerLimit: -360,
  UpperLimit: 360,
  Units: '°',

} ).AddValueSliderField( {
  Name: 'VAng',
  Color: 'green',
  Min: -90,
  Max: 90,
  LowerLimit: -90,
  UpperLimit: 90,
  SnapTo: 0,
  Units: '°',

} ).AddValueSliderField( {
  Name: 'GlobeZoom',
  Color: 'blue',
  Min: 0.25,
  Max: 2,
  SnapTo: [ 0.7, 1 ],

} ).AddValueSliderField( {
  Name: 'FEZoom',
  Color: 'blue',
  Min: 0.25,
  Max: 2,
  SnapTo: 1,

} ).AddValueSliderField( {
  Name: 'Alpha',
  Color: 'gray',
  Min: 0.05,
  Max: 1,

} ).Render();

////////////////////////////////////////

ControlPanels.NewPanel( {
  Name: 'OptionsPanel',
  ModelRef: 'Model',
  NCols: 1,
  OnModelChange: function(field) { Model.Update(field); },

} ).AddCheckboxField( {
  Name: 'OptionField',
  Label: 'Show',
  NCols: 7,
  ColSpan: 3,
  Items: [
    {
      Name: 'ShowContinents',
      Text: 'Continents',
    }, {
      Name: 'ShowGrid',
      Text: 'Grid',
    }, {
      Name: 'ShowCoordSystem',
      Text: 'CoordSystem',
    }, {
      Name: 'ShowVector',
      Text: 'Vector',
    }, {
      Name: 'ShowCartCoord',
      Text: 'CartCoord',
    }, {
      Name: 'ShowGeoCoord',
      Text: 'GeoCoord',
    }, {
      Name: 'ShowData',
      Text: 'Data',
    }
  ]

} ).Render();

</jscript>

WGS84 Model

The WGS84 Module contains many more functions, but only EcefToGeo() and GeoToEcef() are used by the App.

// 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 abbreviations 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; }

// 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)) {
      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;

  },

  :

};

FlatEarthMath Module

The same functions as in the WGS84 module, but for Flat Earth geometry.

var FlatEarthMath = {

  precision: 1e-13,

  // constants

  // radius of the flat earth = distance north pole to south pole
  R: 6371000 * Pi,

  // equator radius of the flat earth = distance north pole to equator
  Re: 6371000 * Pi2,

  // radius of globe earth
  Rglobe: 6371000,

  // gravitational acceleration https://en.wikipedia.org/wiki/Earth
  g: 9.80665,

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

  GeoToEcef: function( lat, lng, h ) {
    // or GeoToEcef( { lat, lng, h } )
    //
    // Converts lat, lng, h to [ x, y, z ].
    // All angles in degrees, h is optional (default = 0)

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

    var lat = rad( lat );
    var lng = rad( lng );
    var r = (Pi2 - lat) * this.Rglobe;
    var x = r * cos( lng );
    var y = r * sin( lng );

    return [ x, y, h ];

  },

  EcefToGeo: function( pos, lngDefault ) {
    // Calculates latitude, longitude and height of pos.
    //
    // If pos is on the north pole, longitude is not defined and can be set by lngDefault.
    // lngDefault is set to 0 if not defined.
    //
    // return { lat, lng, h }
    //
    // All angles are in degrees.

    var lat = 0;
    var lng = 0;
    lngDefault = xDefNum( lngDefault, 0 );

    var posXY = [ pos[0], pos[1], 0 ];
    var r = V3.Length( posXY );
    if (r == 0) {

      // pos is up, so lng is undefined -> set lng = lngDefault
      lat = 90;
      lng = lngDefault;

    } else {

      // assert r > 0 [0..Pi*Rglobe(=20000km)..unlimited)
      lat = deg( Pi2 - r / this.Rglobe );
      lng = acos( limit1( posXY[0] / r ) );
      if (posXY[1] < 0) lng *= -1;
      lng = deg( lng );

    }

    return { lat: lat, lng: lng, h: pos[2] };

  },

  :

};

More Page Infos / Sitemap
Created Thursday, September 30, 2021
Scroll to Top of Page
Changed Sunday, November 7, 2021