Source Code for JavaScript Mapping Tools: geo.js


// Copyright 2006-2008 (c) Paul Demers  <paul@acscdg.com>

// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.

// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA., or visit one 
// of the links here:
//  http://www.gnu.org/licenses/old-licenses/gpl-2.0.txt
//  http://www.acscdg.com/LICENSE.txt

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

//
// JavaScript geography functions.
// Computes great circle courses, rhumb line courses, and
//   rhumb line approximations to great circles.

// Web site with this code running: http://www.acscdg.com/
//
// Dependency on other modules:
//  Google maps (for geo point object).
//  index.js for _optsNotClickable shared object.

// Angle constants.
var piOver4 = Math.PI / 4.0;
var twoPi = Math.PI * 2.0;
var piOver2 = Math.PI / 2.0;
var radiansToDegrees = 180.0 / Math.PI;

// Distance constants.  From Bowditch. 
var kilometersPerNM = 1.852;
var milesPerNM = 1.150779448;

// Drawing settings.
var vectorColor = "#ff0000";
var vectorWeight = 3;     // Was 4
var vectorOpacity = 0.8;  // Default is supposed to be 1.0


//
//// Enumerations to pass around the user's distance measure preference.
var NAUTICAL_MILES_UNITS = 0;
var MILES_UNITS = 1;
var KILOMETERS_UNITS = 2;

//
//// Converts distance measure enumeration to actual units.
function getDistanceMultiplier(distanceUnits)
{
  var distanceMultiplier = 1.0;
  if (distanceUnits == MILES_UNITS)
    distanceMultiplier = milesPerNM;
  else if (distanceUnits == KILOMETERS_UNITS)
    distanceMultiplier = kilometersPerNM;

  return distanceMultiplier;
}



//
//// Reduce an angle to (-PI/2, PI/2), for latitudes.
function reduceLat(lat)  // old name: fng(x)
{
  var reducedLat = lat - Math.PI * Math.floor((lat + piOver2) / Math.PI);
  return reducedLat;
}

//
//// Reduce and angle to (-PI, PI), for longitudes.
function reduceLng(lng)  // old name: fnl(x)
{
  var reducedLng = lng - twoPi * Math.floor((lng + Math.PI ) / twoPi);
  return reducedLng;
}

//
//// Reduce an angle to (0, 2*PI), for direction and azimuth.
function reduceAzimuth(azimuth)
{
  var reducedAzimuth = azimuth - twoPi * Math.floor(azimuth / twoPi);
  return reducedAzimuth;
}


////////////////// Earth Vector Object /////////////////////

// Debugging method.
//function earthVectorToString()
//{
//  return new String("Distance (NM): " + this.distanceNM + ", Azimuth: " + this.azimuthDegrees);
//}


//
//// A vector object contains distance and azimuth.
//// Distance is in nautical miles, azimuth in radians (north is 0, east is PI/2)
function EarthVector(distanceNM, azimuth)
{
  this.distanceNM = distanceNM;
  this.azimuth = reduceAzimuth(azimuth);
  this.azimuthDegrees = azimuth * radiansToDegrees; 
//  this.toString = earthVectorToString;
}


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

// 
//// Utility getter method.
function getVector()
{
  return this.vector;
}


//
//// Computes rhumbline course and distance.  Round Earth formula, for speed.
//// Coordinate system: west is negative, south is negative.
//// Also called the rhumb line direct solution.
function computeRhumbLineVector()
{
  this.deltaLng = this.endPoint.lngRadians() - this.startPoint.lngRadians();

  var gb = Math.tan(piOver4 + (this.endPoint.latRadians() / 2.0));
  var ga = Math.tan(piOver4 + (this.startPoint.latRadians() /2.0));
  var X = Math.log(gb) - Math.log(ga);
  var azimuth = Math.atan2(this.deltaLng, X);

  var distanceNM = 0.0;

  var cosineAzimuth = Math.cos(azimuth);
  if (Math.abs(cosineAzimuth) > 0.0001)
  {
    var deltaLatDegrees = (this.endPoint.lat() - this.startPoint.lat());
    distanceNM = 60.0 * deltaLatDegrees / cosineAzimuth;   // NMi
  }
  else // Too close to straight east or straight west to use cosine.
  {
    var deltaLngDegrees = Math.abs(this.endPoint.lng() - this.startPoint.lng());
    distanceNM = 60.0 * deltaLngDegrees * Math.cos(this.startPoint.lat());
  }

  this.vector = new EarthVector(distanceNM, azimuth);
}


//
//// Creates a polyline overlay object of a rhumbline course.
////  Assumes the line will overlay a mercator map.
function rhumbLineLine()
{
  var rhLine = new GPolyline( [this.startPoint, this.endPoint],
                   vectorColor, vectorWeight, vectorOpacity, _optsNotClickable);
  return rhLine;
}


//
//// Constructor for a rhumbline course object.
function RhumbLineCourse(startPoint, endPoint)
{
   this.startPoint = startPoint;
   this.endPoint = endPoint;

   this.computeRhumbLineVector = computeRhumbLineVector;
   this.getLine = rhumbLineLine;
   this.getVector = getVector;

   this.computeRhumbLineVector();
}


//
//// Computes great circle course and distance.
////  Round Earth formula, for speed.
////  Coordinate system: west is negative, south is negative.
////  Also called the great circle indirect solution.
function computeGreatCircleVector()
{
  this.deltaLng = reduceLng(this.endPoint.lngRadians() - this.startPoint.lngRadians());

  var sinDeltaLng = Math.sin(this.deltaLng);
  var cosDeltaLng = Math.cos(this.deltaLng);

  this.sinStartLat = Math.sin(this.startPoint.latRadians());
  this.cosStartLat = Math.cos(this.startPoint.latRadians());
  this.sinEndLat = Math.sin(this.endPoint.latRadians());
  this.cosEndLat = Math.cos(this.endPoint.latRadians());
  // This is the fastest formula for great circles, but is not the best.
  // The best is the Haversine formula.  But is much more computation.
  this.distanceRadians = Math.acos(this.sinStartLat*this.sinEndLat +
                                   this.cosStartLat*this.cosEndLat*cosDeltaLng);
  this.distanceDegrees = radiansToDegrees * this.distanceRadians;

  var distanceNM = 60.0 * this.distanceDegrees;

  var Y = sinDeltaLng * this.cosEndLat;
  var X = this.cosStartLat * this.sinEndLat - this.sinStartLat * this.cosEndLat * cosDeltaLng;
  azimuth = Math.atan2(Y, X);

  this.vector = new EarthVector(distanceNM, azimuth);
}


//
//// Get latitude from longitude along a great circle course
//// Needed to compute way points along a rhumbline approximation to a great circle.
//// Returns latitude (in radians)
function getLatFromLng(givenLng)
{
  var deltaLng = givenLng - this.startPoint.lngRadians() ;

  var sinDeltaLng = Math.sin(deltaLng);
  var cosDeltaLng = Math.cos(deltaLng);

  var Y = this.sinStartLat * cosDeltaLng * this.sinAzimuth + sinDeltaLng * this.cosAzimuth;
  var X = this.cosStartLat * this.sinAzimuth;
  var foundLat = reduceLat(Math.atan2(Y, X));

  return(foundLat);
}


//
//// Computes the vertex of a great circle course. 
//// The vertex is the point furthest north (or south).
//// Returns a Google LatLng object.
//// TODO: This doesn't always work.  The vertex is on the wrong side
////   of the international date line when crossing from east to west.
function computeVertex()
{
  if (this.sinAzimuth == 0.0)
    return(new GLatLng(twoPi, this.startPoint.lngRadians()));

  var vertexLng = reduceLng(
     Math.atan(1.0/this.sinStartLat/Math.tan(this.vector.azimuth)) + this.startPoint.lngRadians());
  var vertexLat = this.getLatFromLng(vertexLng);

  this.vertex = new GLatLng(vertexLat*radiansToDegrees, vertexLng*radiansToDegrees);

//GLog.write("start: " + this.startPoint + " vertext: " + this.vertex);
}


//
//// Determines the points of the rhumbline approximation to a great circle course.
//// It works by moving a fixed number of longitude degrees, then computing the latitude 
////  along the great circle at that longitude. 
function createWayPoints()
{
  this.sinAzimuth = Math.sin(this.vector.azimuth);
  this.cosAzimuth = Math.cos(this.vector.azimuth);

  this.longitudeIncrementRadians = this.longitudeIncrement/radiansToDegrees;
  var numberOfPointsD = Math.abs(this.deltaLng)/this.longitudeIncrementRadians;

  var numberOfPointsI = Math.floor(numberOfPointsD);

  var wayPointsList = new Array();
  wayPointsList.push(this.startPoint);

  if (numberOfPointsI > 1)
  {
    if (this.deltaLng < 0.0)
      this.longitudeIncrementRadians *= -1.0;  // Always move west to east.

    var nextPointNumber;
    var lastLng = this.startPoint.lngRadians();
    for (nextPointNumber = 0; nextPointNumber < numberOfPointsI; nextPointNumber++)
    {
      var nextLng = reduceLng(lastLng + this.longitudeIncrementRadians);
      var nextLat = this.getLatFromLng(nextLng);

      var nextPoint = new GLatLng(nextLat*radiansToDegrees, nextLng*radiansToDegrees);
  
      wayPointsList.push(nextPoint);

      lastLng = nextLng;

    } // For loop.
  } // If number of points > 1.

  wayPointsList.push(this.endPoint);

  var wayPointsLine = new GPolyline(wayPointsList, 
                       vectorColor, vectorWeight, vectorOpacity, _optsNotClickable) ;
  return(wayPointsLine);
}



//
//// TODO: Combine with rhumb line object.  Have an abstract base, and two sub classes.
function GreatCircleCourse(startPoint, endPoint)
{
   this.longitudeIncrement = 3; // TODO: Perhaps base increment on map scale?

   this.startPoint = startPoint;
   this.endPoint = endPoint;

   this.computeGreatCircleVector = computeGreatCircleVector;
   this.getLatFromLng = getLatFromLng;
   this.computeVertex = computeVertex;
   this.getLine = createWayPoints;
   this.getVector = getVector;

   this.computeGreatCircleVector();
}

////////////////////////// Best Course Object //////////////////////////
// The idea of best course is to use different algorithms based on the 
//  distance between the points.  For example, the haversine formula is 
//  more accurate for small distances.

//
////
function bestCourseGetLine()
{
  if (this.courseLine == null)
    this.courseLine = this.course.getLine();

  return(this.courseLine);
}

//
////
function bestCourseGetVector()
{
  return(this.course.vector);
}

//
//// Chooses the best algorithm, based on distance between points.
function BestCourse(startPoint, endPoint)
{

  this.getLine = bestCourseGetLine;
  this.getVector = bestCourseGetVector;

//  var minDeltaLng = 3;

//  var deltaLng = reduceLng(startPoint.lngRadians() - endPoint.lngRadians());
//  var deltaLngDegrees = Math.abs(deltaLng * radiansToDegrees);

//  if (deltaLngDegrees < minDeltaLng)
//  {
    // The great circle formula used is very fast, but inaccurate 
    //   for short distances. 

//  }
//  else
//  {
      this.course = new GreatCircleCourse(startPoint, endPoint);
//  }
}


//
//// Compute the great circle direct solution.  
//// That is, given a start point, a distance, and a course (azimith),
//// compute the end point. 
//// Returns a Google GLatLng object with the end point.
function directSolution(startPoint, azimuthR, distanceNM)
{
  var distanceR = (distanceNM / 60) / radiansToDegrees;

  var sinStartLat = Math.sin(startPoint.latRadians());
  var cosStartLat = Math.cos(startPoint.latRadians());

  var sinAzimuth = Math.sin(azimuthR);
  var cosAzimuth = Math.cos(azimuthR);

  var sinDistance = Math.sin(distanceR);
  var cosDistance = Math.cos(distanceR);

  var X = (cosStartLat * cosDistance) - (sinStartLat * cosAzimuth * sinDistance);
  var Y = sinDistance * sinAzimuth;

  var endLat = Math.asin((sinStartLat * cosDistance) + (cosStartLat * sinDistance * cosAzimuth));
  var endLng = startPoint.lngRadians() - Math.atan2(Y, X);

// The Google contructor doesn't need reduced angles.
//  endLat = reduceLat(endLat);
//  endLng = reduceLng(endLng);

  return(new GLatLng(endLat*radiansToDegrees, endLng*radiansToDegrees));
}