Compute sun and moon positions

Source/Core/Cartesian3.js

@@ -642,6 +642,43 @@ define([
         return result;
     };
 
+    /**
+     * Produces a Cartesian3 representing this instance which results from rotating
+     * the original axes used to represent this instance by the provided Quaternion rotation.
+     * @memberof Cartesian3
+     *
+     * @param {Cartesian3} cartesian
+     * @param {Quaternion} quaternion
+     * @return {Cartesian3} The result of the rotataion
+     *
+     * @exception {DeveloperError} left is required.
+     * @exception {DeveloperError} right is required.
+     */
+    Cartesian3.rotate = function(cartesian, quaternion, result) {
+        if (typeof cartesian === 'undefined') {
+            throw new DeveloperError('cartesian is required');
+        }
+        if (typeof quaternion === 'undefined') {
+            throw new DeveloperError('quaternion is required');
+        }
+        var w = quaternion.w;
+        var difference = w * w - quaternion.x * quaternion.x - quaternion.y * quaternion.y - quaternion.z * quaternion.z;
+        var dot = cartesian.x * quaternion.x + cartesian.y * quaternion.y + cartesian.z * quaternion.z;
+
+        var x = difference * cartesian.x + 2.0 * (w * (cartesian.y * quaternion.z - cartesian.z * quaternion.y) + dot * quaternion.x);
+        var y = difference * cartesian.y + 2.0 * (w * (cartesian.z * quaternion.x - cartesian.x * quaternion.z) + dot * quaternion.y);
+        var z = difference * cartesian.z + 2.0 * (w * (cartesian.x * quaternion.y - cartesian.y * quaternion.x) + dot * quaternion.z);
+
+        if (typeof result === 'undefined') {
+            return new Cartesian3(x, y, z);
+        }
+
+        result.x = x;
+        result.y = y;
+        result.z = z;
+        return result;
+    };
+
     /**
      * An immutable Cartesian3 instance initialized to (0.0, 0.0, 0.0).
      * @memberof Cartesian3
@@ -926,5 +963,18 @@ define([
         return Cartesian3.cross(this, right, result);
     };
 
+    /**
+     * Produces a Cartesian3 representing this instance which results from rotating
+     * the original axes used to represent this instance by the provided Quaternion rotation.
+     * @memberof Cartesian3
+     *
+     * @param {Quaternion}
+     * @return {Cartesian3} The result of the rotation
+     *
+     */
+    Cartesian3.prototype.rotate = function(quaternion, result) {
+        return Cartesian3.rotate(this, quaternion, result);
+    };
+
     return Cartesian3;
 });

Source/Core/KeplerianElements.js

@@ -0,0 +1,175 @@
+/*global define*/
+define([
+        './DeveloperError',
+        './Math',
+        './Cartesian3',
+        './Matrix3'],
+    function(
+        DeveloperError,
+        CesiumMath,
+        Cartesian3,
+        Matrix3) {
+    "use strict";
+
+    /**
+     * Functions related to classical Keplerian elements representing an elliptical orbit.
+     *
+     * @see ModifiedKeplerianElements
+     */
+
+    var KeplerianElements = {};
+
+    /**
+     * Calculates the true anomaly given the mean anomaly and the eccentricity.
+     * @param {Number} meanAnomaly The mean anomaly (radians)
+     * @param {Number} eccentricity The eccentricity
+     *
+     * @exception {DeveloperError} This method assumes a circular or elliptical orbit, so the eccentricity must be between zero and unity.
+
+     *
+     * @returns {Number} The true anomaly (radians)
+     */
+    KeplerianElements.MeanAnomalyToTrueAnomaly = function(meanAnomaly, eccentricity) {
+        if (eccentricity < 0.0 || eccentricity >= 1.0)
+        {
+            throw new DeveloperError('eccentricity out of range.');
+        }
+        var eccentricAnomaly = this.MeanAnomalyToEccentricAnomaly(meanAnomaly, eccentricity);
+        return this.EccentricAnomalyToTrueAnomaly(eccentricAnomaly, eccentricity);
+    };
+
+    /**
+     * Calculates the eccentric anomaly given the mean anomaly and the eccentricity.
+     * @param {Number} meanAnomaly The mean anomaly (radians)
+     * @param {Number} eccentricity The eccentricity
+     *
+     * @exception {DeveloperError} This method assumes a circular or elliptical orbit, so the eccentricity must be between zero and unity.
+     *
+     * @returns {Number} The eccentric anomaly (radians)
+     */
+    KeplerianElements.MeanAnomalyToEccentricAnomaly = function(meanAnomaly, eccentricity)
+    {
+        if (eccentricity < 0.0 || eccentricity >= 1.0)
+        {
+            throw new DeveloperError('eccentricity out of range.');
+        }
+
+        // Maximum number of times to iterate in on Kepler's equation
+        var maxIterationCount = 50;
+
+        // Maximum difference to be considered convergence of Kepler's equation
+        var keplerEqConvergence = CesiumMath.EPSILON8;
+
+        var revs = Math.floor(meanAnomaly / CesiumMath.TWO_PI);
+
+        // Find angle in current revolution
+        meanAnomaly -= revs * CesiumMath.TWO_PI;
+
+        // calculate starting value for iteration sequence
+        var iterationValue = meanAnomaly + (eccentricity * Math.sin(meanAnomaly)) /
+            (1.0 - Math.sin(meanAnomaly + eccentricity) + Math.sin(meanAnomaly));
+
+        // Perform Newton-Raphson iteration on Kepler's equation
+        var eccentricAnomaly = Number.MAX_VALUE;
+
+        var count;
+        for (count = 0;
+            count < maxIterationCount && Math.abs(eccentricAnomaly - iterationValue) > keplerEqConvergence;
+            ++count)
+        {
+            eccentricAnomaly = iterationValue;
+            var NRfunction = eccentricAnomaly - eccentricity * Math.sin(eccentricAnomaly) - meanAnomaly;
+            var dNRfunction = 1 - eccentricity * Math.cos(eccentricAnomaly);
+            iterationValue = eccentricAnomaly - NRfunction / dNRfunction;
+        }
+
+        if (count >= maxIterationCount) {
+            throw new DeveloperError('Kepler equation did not converge');
+            //TODO: Port 'DoubleFunctionExplorer' from components
+        }
+
+        eccentricAnomaly = iterationValue + revs * CesiumMath.TWO_PI;
+        return eccentricAnomaly;
+    };
+
+    /**
+     * Calculates the true anomaly given the eccentric anomaly and the eccentricity.
+     * @param {Number} eccentricAnomaly The eccentric anomaly (radians)
+     * @param {Number} eccentricity The eccentricity
+     *
+     * @exception {DeveloperError} This method assumes a circular or elliptical orbit, so the eccentricity must be between zero and unity.
+     *
+     * @returns {Number} The true anomaly (radians)
+     */
+    KeplerianElements.EccentricAnomalyToTrueAnomaly = function(eccentricAnomaly, eccentricity) {
+        if (eccentricity < 0.0 || eccentricity >= 1.0)
+        {
+            throw new DeveloperError('eccentricity out of range.');
+        }
+
+        // Calculate the number of previous revolutions
+        var revs = Math.floor(eccentricAnomaly / CesiumMath.TWO_PI);
+
+        // Find angle in current revolution
+        eccentricAnomaly -= revs * CesiumMath.TWO_PI;
+
+        // Calculate true anomaly from eccentric anomaly
+        var trueAnomalyX = Math.cos(eccentricAnomaly) - eccentricity;
+        var trueAnomalyY = Math.sin(eccentricAnomaly) * Math.sqrt(1 - eccentricity * eccentricity);
+
+        var trueAnomaly = Math.atan2(trueAnomalyY, trueAnomalyX);
+
+        // Ensure the correct quadrant
+        trueAnomaly = CesiumMath.zeroToTwoPi(trueAnomaly);
+        if (eccentricAnomaly < 0)
+        {
+            trueAnomaly -= CesiumMath.TWO_PI;
+        }
+
+        // Add on previous revolutions
+        trueAnomaly += revs * CesiumMath.TWO_PI;
+
+        return trueAnomaly;
+    };
+
+    /**
+     * Calculates the transformation matrix to convert from the perifocal (PQW) coordinate
+     * system to inertial cartesian coordinates.
+     *
+     * @param {Number} argumentOfPeriapsis The argument of periapsis (radians)
+     * @param {Number} inclination The inclination (radians)
+     * @param {Number} rightAscension The right ascension of the ascending node (radians)
+     *
+     * @returns {Matrix3} The transformation matrix to convert from the perifocal coordinate system to inertial cartesian coordinates.
+     */
+    KeplerianElements.PerifocalToCartesianMatrix = function(argumentOfPeriapsis, inclination, rightAscension) {
+        if (inclination < 0 || inclination > CesiumMath.PI)
+        {
+            throw new DeveloperError('inclination out of range');
+        }
+        var cosap = Math.cos(argumentOfPeriapsis);
+        var sinap = Math.sin(argumentOfPeriapsis);
+
+        var cosi = Math.cos(inclination);
+        var sini = Math.sin(inclination);
+
+        var cosraan = Math.cos(rightAscension);
+        var sinraan = Math.sin(rightAscension);
+
+        return new Matrix3(
+            cosraan * cosap - sinraan * sinap * cosi,
+            -cosraan * sinap - sinraan * cosap * cosi,
+            sinraan * sini,
+
+            sinraan * cosap + cosraan * sinap * cosi,
+            -sinraan * sinap + cosraan * cosap * cosi,
+            -cosraan * sini,
+
+            sinap * sini,
+            cosap * sini,
+            cosi);
+    };
+
+    return KeplerianElements;
+
+});

Source/Core/Math.js

@@ -417,6 +417,16 @@ define([
         return x > pi ? pi : x;
     };
 
+    /**
+     * Alters the value of input x such that 0 <= x <= <code>CesiumMath.TWO_PI</code>
+     * @param {Number} angle in radians
+     * @return {Number} The angle in the range (0 , <code>CesiumMath.TWO_PI</code>).
+     */
+    CesiumMath.zeroToTwoPi = function(x) {
+        var value = x % CesiumMath.TWO_PI;
+        return (value < 0.0) ? (value + CesiumMath.TWO_PI) % CesiumMath.TWO_PI : value;
+    };
+
     /**
      * DOC_TBA
      */

Source/Core/ModifiedKeplerianElements.js

@@ -0,0 +1,116 @@
+/*global define*/
+define([
+        './DeveloperError',
+        './Math',
+        './Cartesian3',
+        './Matrix3',
+        './KeplerianElements'
+    ],
+    function(
+        DeveloperError,
+        CesiumMath,
+        Cartesian3,
+        Matrix3,
+        KeplerianElements) {
+    "use strict";
+
+    /**
+     * Modified Keplerian orbital elements.  These are the same as the Classical/Keplerian orbital elements
+     * except that Radius of Periapsis and the inverse of Semimajor Axis are used instead of
+     * Semimajor Axis and Eccentricity.  This is useful because the Radius of Periapsis is well defined
+     * for all but rectilinear orbits.
+     * @see KeplerianElements
+     */
+
+
+    /**
+     *  Initialize a set of modified Keplerian elements.
+     * @alias ModifiedKeperianElements
+     * @constructor
+     *
+     * @param {Number} radiusOfPeriapsis Radius of periapsis (distance)
+     * @param {Number} inverseSemimajorAxis The inverse of semimajor axis (distance)
+     * @param {Number} inclination Inclination (radians)
+     * @param {Number} argumentOfPeriapsis Argument of periapsis (radians)
+     * @param {Number} rightAscensionOfAscendingNode Right ascension of the ascending node (radians)
+     * @param {Number} trueAnomaly True anomaly (radians)
+     * @param {Number} gravitationalParameter The gravitational parameter associated with these elements (distance cubed per time squared)
+     *
+     */
+    var ModifiedKeplerianElements = function(radiusOfPeriapsis, inverseSemimajorAxis, inclination, argumentOfPeriapsis,
+            rightAscensionOfAscendingNode, trueAnomaly, gravitationalParameter) {
+        if (inclination < 0 || inclination > CesiumMath.PI) {
+            throw new DeveloperError('inclination out of range.');
+        }
+        this.radiusOfPeriapsis = radiusOfPeriapsis;
+        this.inverseSemimajorAxis = inverseSemimajorAxis;
+        this.inclination = inclination;
+        this.eccentricity = 1.0 - radiusOfPeriapsis * inverseSemimajorAxis;
+        this.argumentOfPeriapsis = argumentOfPeriapsis;
+        this.rightAscensionOfAscendingNode = rightAscensionOfAscendingNode;
+        this.trueAnomaly = trueAnomaly;
+        this.gravitationalParameter = gravitationalParameter;
+        this.type = chooseOrbit(this.eccentricity, 0.0);
+        if (this.type === 'Hyperbolic' && Math.abs(CesiumMath.NegativePiToPi(this.trueAnomaly)) >= Math.acos(- 1.0 / this.eccentricity)) {
+            throw new DeveloperError('invalid trueAnomaly.');
+        }
+    };
+
+    function chooseOrbit(eccentricity, tolerance)
+    {
+        if (eccentricity < 0)
+        {
+            throw new DeveloperError('eccentricity cannot be negative.');
+        }
+        else if (eccentricity <= tolerance)
+        {
+            return 'Circular';
+        }
+        else if (eccentricity < 1.0 - tolerance)
+        {
+            return 'Elliptical';
+        }
+        else if (eccentricity <= 1.0 + tolerance)
+        {
+            return 'Parabolic';
+        }
+        else
+        {
+            return 'Hyperbolic';
+        }
+    }
+
+    /**
+     * Returns a Cartesian representation of these orbital elements.
+     * @param {ModifiedKeplerianElements} element
+     * @returns {Cartesian3} Equivalent Cartesian
+     */
+    ModifiedKeplerianElements.ToCartesian = function(element) {
+        var perifocalToEquatorial = KeplerianElements.PerifocalToCartesianMatrix(element.argumentOfPeriapsis, element.inclination, element.rightAscensionOfAscendingNode);
+        var semilatus = element.radiusOfPeriapsis * (1.0 + element.eccentricity);
+        var costheta = Math.cos(element.trueAnomaly);
+        var sintheta = Math.sin(element.trueAnomaly);
+
+        var denom = (1.0 + element.eccentricity * costheta);
+        if (denom <= CesiumMath.Epsilon10) {
+            throw new DeveloperError('elements cannot be converted to cartesian');
+        }
+
+        var radius = semilatus / denom;
+        var position = new Cartesian3(radius * costheta, radius * sintheta, 0.0);
+
+        return perifocalToEquatorial.multiplyByVector(position);
+    };
+
+
+    /**
+     * Returns a Cartesian representation of these orbital elements.
+     *
+     * @returns {Cartesian3} Equivalent Cartesian
+     */
+    ModifiedKeplerianElements.prototype.ToCartesian = function(){
+        return ModifiedKeplerianElements.ToCartesian(this);
+    };
+
+    return ModifiedKeplerianElements;
+});