/*global define*/
define([
'./Cartesian3',
'./Cartographic',
'./defaultValue',
'./defined',
'./defineProperties',
'./DeveloperError',
'./Ellipsoid',
'./Math'
], function(
Cartesian3,
Cartographic,
defaultValue,
defined,
defineProperties,
DeveloperError,
Ellipsoid,
CesiumMath) {
'use strict';
function setConstants(ellipsoidGeodesic) {
var uSquared = ellipsoidGeodesic._uSquared;
var a = ellipsoidGeodesic._ellipsoid.maximumRadius;
var b = ellipsoidGeodesic._ellipsoid.minimumRadius;
var f = (a - b) / a;
var cosineHeading = Math.cos(ellipsoidGeodesic._startHeading);
var sineHeading = Math.sin(ellipsoidGeodesic._startHeading);
var tanU = (1 - f) * Math.tan(ellipsoidGeodesic._start.latitude);
var cosineU = 1.0 / Math.sqrt(1.0 + tanU * tanU);
var sineU = cosineU * tanU;
var sigma = Math.atan2(tanU, cosineHeading);
var sineAlpha = cosineU * sineHeading;
var sineSquaredAlpha = sineAlpha * sineAlpha;
var cosineSquaredAlpha = 1.0 - sineSquaredAlpha;
var cosineAlpha = Math.sqrt(cosineSquaredAlpha);
var u2Over4 = uSquared / 4.0;
var u4Over16 = u2Over4 * u2Over4;
var u6Over64 = u4Over16 * u2Over4;
var u8Over256 = u4Over16 * u4Over16;
var a0 = (1.0 + u2Over4 - 3.0 * u4Over16 / 4.0 + 5.0 * u6Over64 / 4.0 - 175.0 * u8Over256 / 64.0);
var a1 = (1.0 - u2Over4 + 15.0 * u4Over16 / 8.0 - 35.0 * u6Over64 / 8.0);
var a2 = (1.0 - 3.0 * u2Over4 + 35.0 * u4Over16 / 4.0);
var a3 = (1.0 - 5.0 * u2Over4);
var distanceRatio = a0 * sigma - a1 * Math.sin(2.0 * sigma) * u2Over4 / 2.0 - a2 * Math.sin(4.0 * sigma) * u4Over16 / 16.0 -
a3 * Math.sin(6.0 * sigma) * u6Over64 / 48.0 - Math.sin(8.0 * sigma) * 5.0 * u8Over256 / 512;
var constants = ellipsoidGeodesic._constants;
constants.a = a;
constants.b = b;
constants.f = f;
constants.cosineHeading = cosineHeading;
constants.sineHeading = sineHeading;
constants.tanU = tanU;
constants.cosineU = cosineU;
constants.sineU = sineU;
constants.sigma = sigma;
constants.sineAlpha = sineAlpha;
constants.sineSquaredAlpha = sineSquaredAlpha;
constants.cosineSquaredAlpha = cosineSquaredAlpha;
constants.cosineAlpha = cosineAlpha;
constants.u2Over4 = u2Over4;
constants.u4Over16 = u4Over16;
constants.u6Over64 = u6Over64;
constants.u8Over256 = u8Over256;
constants.a0 = a0;
constants.a1 = a1;
constants.a2 = a2;
constants.a3 = a3;
constants.distanceRatio = distanceRatio;
}
function computeC(f, cosineSquaredAlpha) {
return f * cosineSquaredAlpha * (4.0 + f * (4.0 - 3.0 * cosineSquaredAlpha)) / 16.0;
}
function computeDeltaLambda(f, sineAlpha, cosineSquaredAlpha, sigma, sineSigma, cosineSigma, cosineTwiceSigmaMidpoint) {
var C = computeC(f, cosineSquaredAlpha);
return (1.0 - C) * f * sineAlpha * (sigma + C * sineSigma * (cosineTwiceSigmaMidpoint +
C * cosineSigma * (2.0 * cosineTwiceSigmaMidpoint * cosineTwiceSigmaMidpoint - 1.0)));
}
function vincentyInverseFormula(ellipsoidGeodesic, major, minor, firstLongitude, firstLatitude, secondLongitude, secondLatitude) {
var eff = (major - minor) / major;
var l = secondLongitude - firstLongitude;
var u1 = Math.atan((1 - eff) * Math.tan(firstLatitude));
var u2 = Math.atan((1 - eff) * Math.tan(secondLatitude));
var cosineU1 = Math.cos(u1);
var sineU1 = Math.sin(u1);
var cosineU2 = Math.cos(u2);
var sineU2 = Math.sin(u2);
var cc = cosineU1 * cosineU2;
var cs = cosineU1 * sineU2;
var ss = sineU1 * sineU2;
var sc = sineU1 * cosineU2;
var lambda = l;
var lambdaDot = CesiumMath.TWO_PI;
var cosineLambda = Math.cos(lambda);
var sineLambda = Math.sin(lambda);
var sigma;
var cosineSigma;
var sineSigma;
var cosineSquaredAlpha;
var cosineTwiceSigmaMidpoint;
do {
cosineLambda = Math.cos(lambda);
sineLambda = Math.sin(lambda);
var temp = cs - sc * cosineLambda;
sineSigma = Math.sqrt(cosineU2 * cosineU2 * sineLambda * sineLambda + temp * temp);
cosineSigma = ss + cc * cosineLambda;
sigma = Math.atan2(sineSigma, cosineSigma);
var sineAlpha;
if (sineSigma === 0.0) {
sineAlpha = 0.0;
cosineSquaredAlpha = 1.0;
} else {
sineAlpha = cc * sineLambda / sineSigma;
cosineSquaredAlpha = 1.0 - sineAlpha * sineAlpha;
}
lambdaDot = lambda;
cosineTwiceSigmaMidpoint = cosineSigma - 2.0 * ss / cosineSquaredAlpha;
if (isNaN(cosineTwiceSigmaMidpoint)) {
cosineTwiceSigmaMidpoint = 0.0;
}
lambda = l + computeDeltaLambda(eff, sineAlpha, cosineSquaredAlpha,
sigma, sineSigma, cosineSigma, cosineTwiceSigmaMidpoint);
} while (Math.abs(lambda - lambdaDot) > CesiumMath.EPSILON12);
var uSquared = cosineSquaredAlpha * (major * major - minor * minor) / (minor * minor);
var A = 1.0 + uSquared * (4096.0 + uSquared * (uSquared * (320.0 - 175.0 * uSquared) - 768.0)) / 16384.0;
var B = uSquared * (256.0 + uSquared * (uSquared * (74.0 - 47.0 * uSquared) - 128.0)) / 1024.0;
var cosineSquaredTwiceSigmaMidpoint = cosineTwiceSigmaMidpoint * cosineTwiceSigmaMidpoint;
var deltaSigma = B * sineSigma * (cosineTwiceSigmaMidpoint + B * (cosineSigma *
(2.0 * cosineSquaredTwiceSigmaMidpoint - 1.0) - B * cosineTwiceSigmaMidpoint *
(4.0 * sineSigma * sineSigma - 3.0) * (4.0 * cosineSquaredTwiceSigmaMidpoint - 3.0) / 6.0) / 4.0);
var distance = minor * A * (sigma - deltaSigma);
var startHeading = Math.atan2(cosineU2 * sineLambda, cs - sc * cosineLambda);
var endHeading = Math.atan2(cosineU1 * sineLambda, cs * cosineLambda - sc);
ellipsoidGeodesic._distance = distance;
ellipsoidGeodesic._startHeading = startHeading;
ellipsoidGeodesic._endHeading = endHeading;
ellipsoidGeodesic._uSquared = uSquared;
}
function computeProperties(ellipsoidGeodesic, start, end, ellipsoid) {
var firstCartesian = Cartesian3.normalize(ellipsoid.cartographicToCartesian(start, scratchCart2), scratchCart1);
var lastCartesian = Cartesian3.normalize(ellipsoid.cartographicToCartesian(end, scratchCart2), scratchCart2);
//>>includeStart('debug', pragmas.debug);
if (Math.abs(Math.abs(Cartesian3.angleBetween(firstCartesian, lastCartesian)) - Math.PI) < 0.0125) {
throw new DeveloperError('geodesic position is not unique');
}
//>>includeEnd('debug');
vincentyInverseFormula(ellipsoidGeodesic, ellipsoid.maximumRadius, ellipsoid.minimumRadius,
start.longitude, start.latitude, end.longitude, end.latitude);
ellipsoidGeodesic._start = Cartographic.clone(start, ellipsoidGeodesic._start);
ellipsoidGeodesic._end = Cartographic.clone(end, ellipsoidGeodesic._end);
ellipsoidGeodesic._start.height = 0;
ellipsoidGeodesic._end.height = 0;
setConstants(ellipsoidGeodesic);
}
var scratchCart1 = new Cartesian3();
var scratchCart2 = new Cartesian3();
/**
* Initializes a geodesic on the ellipsoid connecting the two provided planetodetic points.
*
* @alias EllipsoidGeodesic
* @constructor
*
* @param {Cartographic} [start] The initial planetodetic point on the path.
* @param {Cartographic} [end] The final planetodetic point on the path.
* @param {Ellipsoid} [ellipsoid=Ellipsoid.WGS84] The ellipsoid on which the geodesic lies.
*/
function EllipsoidGeodesic(start, end, ellipsoid) {
var e = defaultValue(ellipsoid, Ellipsoid.WGS84);
this._ellipsoid = e;
this._start = new Cartographic();
this._end = new Cartographic();
this._constants = {};
this._startHeading = undefined;
this._endHeading = undefined;
this._distance = undefined;
this._uSquared = undefined;
if (defined(start) && defined(end)) {
computeProperties(this, start, end, e);
}
}
defineProperties(EllipsoidGeodesic.prototype, {
/**
* Gets the ellipsoid.
* @memberof EllipsoidGeodesic.prototype
* @type {Ellipsoid}
* @readonly
*/
ellipsoid : {
get : function() {
return this._ellipsoid;
}
},
/**
* Gets the surface distance between the start and end point
* @memberof EllipsoidGeodesic.prototype
* @type {Number}
* @readonly
*/
surfaceDistance : {
get : function() {
//>>includeStart('debug', pragmas.debug);
if (!defined(this._distance)) {
throw new DeveloperError('set end positions before getting surfaceDistance');
}
//>>includeEnd('debug');
return this._distance;
}
},
/**
* Gets the initial planetodetic point on the path.
* @memberof EllipsoidGeodesic.prototype
* @type {Cartographic}
* @readonly
*/
start : {
get : function() {
return this._start;
}
},
/**
* Gets the final planetodetic point on the path.
* @memberof EllipsoidGeodesic.prototype
* @type {Cartographic}
* @readonly
*/
end : {
get : function() {
return this._end;
}
},
/**
* Gets the heading at the initial point.
* @memberof EllipsoidGeodesic.prototype
* @type {Number}
* @readonly
*/
startHeading : {
get : function() {
//>>includeStart('debug', pragmas.debug);
if (!defined(this._distance)) {
throw new DeveloperError('set end positions before getting startHeading');
}
//>>includeEnd('debug');
return this._startHeading;
}
},
/**
* Gets the heading at the final point.
* @memberof EllipsoidGeodesic.prototype
* @type {Number}
* @readonly
*/
endHeading : {
get : function() {
//>>includeStart('debug', pragmas.debug);
if (!defined(this._distance)) {
throw new DeveloperError('set end positions before getting endHeading');
}
//>>includeEnd('debug');
return this._endHeading;
}
}
});
/**
* Sets the start and end points of the geodesic
*
* @param {Cartographic} start The initial planetodetic point on the path.
* @param {Cartographic} end The final planetodetic point on the path.
*/
EllipsoidGeodesic.prototype.setEndPoints = function(start, end) {
//>>includeStart('debug', pragmas.debug);
if (!defined(start)) {
throw new DeveloperError('start cartographic position is required');
}
if (!defined(end)) {
throw new DeveloperError('end cartgraphic position is required');
}
//>>includeEnd('debug');
computeProperties(this, start, end, this._ellipsoid);
};
/**
* Provides the location of a point at the indicated portion along the geodesic.
*
* @param {Number} fraction The portion of the distance between the initial and final points.
* @returns {Cartographic} The location of the point along the geodesic.
*/
EllipsoidGeodesic.prototype.interpolateUsingFraction = function(fraction, result) {
return this.interpolateUsingSurfaceDistance(this._distance * fraction, result);
};
/**
* Provides the location of a point at the indicated distance along the geodesic.
*
* @param {Number} distance The distance from the inital point to the point of interest along the geodesic
* @returns {Cartographic} The location of the point along the geodesic.
*
* @exception {DeveloperError} start and end must be set before calling function interpolateUsingSurfaceDistance
*/
EllipsoidGeodesic.prototype.interpolateUsingSurfaceDistance = function(distance, result) {
//>>includeStart('debug', pragmas.debug);
if (!defined(this._distance)) {
throw new DeveloperError('start and end must be set before calling function interpolateUsingSurfaceDistance');
}
//>>includeEnd('debug');
var constants = this._constants;
var s = constants.distanceRatio + distance / constants.b;
var cosine2S = Math.cos(2.0 * s);
var cosine4S = Math.cos(4.0 * s);
var cosine6S = Math.cos(6.0 * s);
var sine2S = Math.sin(2.0 * s);
var sine4S = Math.sin(4.0 * s);
var sine6S = Math.sin(6.0 * s);
var sine8S = Math.sin(8.0 * s);
var s2 = s * s;
var s3 = s * s2;
var u8Over256 = constants.u8Over256;
var u2Over4 = constants.u2Over4;
var u6Over64 = constants.u6Over64;
var u4Over16 = constants.u4Over16;
var sigma = 2.0 * s3 * u8Over256 * cosine2S / 3.0 +
s * (1.0 - u2Over4 + 7.0 * u4Over16 / 4.0 - 15.0 * u6Over64 / 4.0 + 579.0 * u8Over256 / 64.0 -
(u4Over16 - 15.0 * u6Over64 / 4.0 + 187.0 * u8Over256 / 16.0) * cosine2S -
(5.0 * u6Over64 / 4.0 - 115.0 * u8Over256 / 16.0) * cosine4S -
29.0 * u8Over256 * cosine6S / 16.0) +
(u2Over4 / 2.0 - u4Over16 + 71.0 * u6Over64 / 32.0 - 85.0 * u8Over256 / 16.0) * sine2S +
(5.0 * u4Over16 / 16.0 - 5.0 * u6Over64 / 4.0 + 383.0 * u8Over256 / 96.0) * sine4S -
s2 * ((u6Over64 - 11.0 * u8Over256 / 2.0) * sine2S + 5.0 * u8Over256 * sine4S / 2.0) +
(29.0 * u6Over64 / 96.0 - 29.0 * u8Over256 / 16.0) * sine6S +
539.0 * u8Over256 * sine8S / 1536.0;
var theta = Math.asin(Math.sin(sigma) * constants.cosineAlpha);
var latitude = Math.atan(constants.a / constants.b * Math.tan(theta));
// Redefine in terms of relative argument of latitude.
sigma = sigma - constants.sigma;
var cosineTwiceSigmaMidpoint = Math.cos(2.0 * constants.sigma + sigma);
var sineSigma = Math.sin(sigma);
var cosineSigma = Math.cos(sigma);
var cc = constants.cosineU * cosineSigma;
var ss = constants.sineU * sineSigma;
var lambda = Math.atan2(sineSigma * constants.sineHeading, cc - ss * constants.cosineHeading);
var l = lambda - computeDeltaLambda(constants.f, constants.sineAlpha, constants.cosineSquaredAlpha,
sigma, sineSigma, cosineSigma, cosineTwiceSigmaMidpoint);
if (defined(result)) {
result.longitude = this._start.longitude + l;
result.latitude = latitude;
result.height = 0.0;
return result;
}
return new Cartographic(this._start.longitude + l, latitude, 0.0);
};
return EllipsoidGeodesic;
});