Compute sun and moon positions

CHANGES.md

@@ -29,6 +29,7 @@ Beta Releases
          });
    * `CzmlCartographic` has been removed and all cartographic values are converted to Cartesian internally during CZML processing.  This improves performance and fixes interpolation of cartographic source data.  The Cartographic representation can still be retrieved if needed.
    * Removed `ComplexConicSensorVolume`, which was not documented and did not work on most platforms.  It will be brought back in a future release.  This does not affect CZML, which uses a custom sensor to approximate a complex conic.
+   * Replaced `computeSunPosition` with `Simon1994PlanetaryPosition`, which has functions to calculate the position of the sun and the moon more accurately.
 * Added wide polylines that work with and without ANGLE.
 * Polylines now use materials to describe their surface appearance. See the [Fabric](https://github.com/AnalyticalGraphicsInc/cesium/wiki/Fabric) wiki page for more details on how to create materials.
 * Added new `PolylineOutline`, `PolylineArrow`, and `Fade` materials.

Source/Core/Math.js

@@ -397,7 +397,7 @@ define([
     };
 
     /**
-     * Alters the value of input x such that <code>-CesiumMath.PI</code> <= x <= <code>CesiumMath.PI</code>
+     * Produces an angle in the range 0 <= angle <= 2Pi which is equivalent to the provided angle.
      * @param {Number} angle in radians
      * @return {Number} The angle in the range ()<code>-CesiumMath.PI</code>, <code>CesiumMath.PI</code>).
      */
@@ -417,6 +417,18 @@ define([
         return x > pi ? pi : x;
     };
 
+    /**
+     * Produces an angle in the range -Pi <= angle <= Pi which is equivalent to the provided angle.
+     * @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;
+        // We do a second modules here if we add 2Pi to ensure that we don't have any numerical issues with very
+        // small negative values.
+        return (value < 0.0) ? (value + CesiumMath.TWO_PI) % CesiumMath.TWO_PI : value;
+    };
+
     /**
      * DOC_TBA
      */

Source/Core/Transforms.js

@@ -275,7 +275,7 @@ define([
      * function updateAndRender() {
      *     var now = new JulianDate();
      *     scene.initializeFrame();
-     *     scene.setSunPosition(computeSunPosition(now));
+     *     scene.setSunPosition(Simon1994PlanetaryPositions.ComputeSunPositionInEarthInertialFrame(now));
      *     scene.getCamera().transform = Matrix4.fromRotationTranslation(Transforms.computeTemeToPseudoFixedMatrix(now), Cartesian3.ZERO);
      *     scene.render();
      *     requestAnimationFrame(updateAndRender);
@@ -413,7 +413,7 @@ define([
      * function updateAndRender() {
      *     var now = new JulianDate();
      *     scene.initializeFrame();
-     *     scene.setSunPosition(computeSunPosition(now));
+     *     scene.setSunPosition(Simon1994PlanetaryPositions.ComputeSunPositionInEarthInertialFrame(now));
      *     var icrfToFixed = Transforms.computeIcrfToFixedMatrix(now);
      *     if (typeof icrfToFixed !== 'undefined') {
      *         scene.getCamera().transform = Matrix4.fromRotationTranslation(icrfToFixed, Cartesian3.ZERO);

Source/Renderer/UniformState.js

@@ -13,7 +13,7 @@ define([
         '../Core/EncodedCartesian3',
         '../Core/BoundingRectangle',
         '../Core/Transforms',
-        '../Core/computeSunPosition',
+        '../Core/Simon1994PlanetaryPositions',
         '../Scene/SceneMode'
     ], function(
         DeveloperError,
@@ -29,7 +29,7 @@ define([
         EncodedCartesian3,
         BoundingRectangle,
         Transforms,
-        computeSunPosition,
+        Simon1994PlanetaryPositions,
         SceneMode) {
     "use strict";
 
@@ -182,18 +182,23 @@ define([
         uniformState._encodedCameraPositionMCDirty = true;
     }
 
-    var sunPositionWC = new Cartesian3();
-    var sunPositionScratch = new Cartesian3();
-
+    var position = new Cartesian3();
+    var transformMatrix = new Matrix3();
     function setSunAndMoonDirections(uniformState, frameState) {
-        computeSunPosition(frameState.time, sunPositionWC);
+        if (typeof Transforms.computeIcrfToFixedMatrix(frameState.time, transformMatrix) === 'undefined') {
+            transformMatrix = Transforms.computeTemeToPseudoFixedMatrix(frameState.time, transformMatrix);
+        }
 
-        Cartesian3.normalize(sunPositionWC, uniformState._sunDirectionWC);
-        Matrix3.multiplyByVector(uniformState.getViewRotation3D(), sunPositionWC, sunPositionScratch);
-        Cartesian3.normalize(sunPositionScratch, uniformState._sunDirectionEC);
+        position = Simon1994PlanetaryPositions.ComputeSunPositionInEarthInertialFrame(frameState.time, position);
+        transformMatrix.multiplyByVector(position, position);
+        Cartesian3.normalize(position, uniformState._sunDirectionWC);
+        Matrix3.multiplyByVector(uniformState.getViewRotation3D(), position, position);
+        Cartesian3.normalize(position, uniformState._sunDirectionEC);
 
-        // Pseudo direction for now just for lighting
-        Cartesian3.negate(uniformState._sunDirectionEC, uniformState._moonDirectionEC);
+        position = Simon1994PlanetaryPositions.ComputeMoonPositionInEarthInertialFrame(frameState.time, position);
+        transformMatrix.multiplyByVector(position, position);
+        Matrix3.multiplyByVector(uniformState.getViewRotation3D(), position, position);
+        Cartesian3.normalize(position, uniformState._moonDirectionEC);
     }
 
     /**

Specs/Core/Simon1994PlanetaryPositionsSpec.js

@@ -0,0 +1,113 @@
+/*global defineSuite*/
+defineSuite(['Core/Simon1994PlanetaryPositions',
+             'Core/JulianDate',
+             'Core/TimeStandard',
+             'Core/TimeConstants',
+             'Core/Math',
+             'Core/Matrix3',
+             'Core/Cartesian3',
+             'Core/Transforms'],
+function(PlanetaryPositions,
+         JulianDate,
+         TimeStandard,
+         TimeConstants,
+         CesiumMath,
+         Matrix3,
+         Cartesian3,
+         Transforms) {
+    "use strict";
+    /*global jasmine,describe,xdescribe,it,xit,expect,beforeEach,afterEach,beforeAll,afterAll,spyOn,runs,waits,waitsFor*/
+
+
+    // Values for the X Y and Z were found using the STK Components GeometryTransformer on the position of the
+    // sun center of mass point and the earth J2000 reference frame.
+    it('computes correct sun position', function() {
+        var date = JulianDate.fromTotalDays(2451545.0, TimeStandard.TAI);
+        var sun = PlanetaryPositions.ComputeSunPositionInEarthInertialFrame(date);
+        var X = 26500268539.790234;
+        var Y = -132756447253.27325;
+        var Z = -57556483362.533806;
+        expect(X).toEqualEpsilon(sun.x, CesiumMath.EPSILON4); //TODO
+        expect(Y).toEqualEpsilon(sun.y, CesiumMath.EPSILON4);
+        expect(Z).toEqualEpsilon(sun.z, CesiumMath.EPSILON4);
+
+        date = JulianDate.fromTotalDays(2456401.5, TimeStandard.TAI);
+        sun = PlanetaryPositions.ComputeSunPositionInEarthInertialFrame(date);
+        X = 131512388940.33589;
+        Y = 66661342667.949928;
+        Z = 28897975607.905258;
+        expect(X).toEqualEpsilon(sun.x, CesiumMath.EPSILON4);
+        expect(Y).toEqualEpsilon(sun.y, CesiumMath.EPSILON4);
+        expect(Z).toEqualEpsilon(sun.z, CesiumMath.EPSILON4);
+
+        date = JulianDate.fromTotalDays(2455998.591667, TimeStandard.TAI);
+        sun = PlanetaryPositions.ComputeSunPositionInEarthInertialFrame(date);
+        X = 147109989956.19534;
+        Y = -19599996881.217579;
+        Z = -8497578102.7696457;
+        expect(X).toEqualEpsilon(sun.x, CesiumMath.EPSILON4);
+        expect(Y).toEqualEpsilon(sun.y, CesiumMath.EPSILON4);
+        expect(Z).toEqualEpsilon(sun.z, CesiumMath.EPSILON4);
+    });
+
+    // Values for X Y and Z were found using the STK Components GeometryTransformer on the Simon 1994 moon point and the earth
+    // J2000 reference frame.
+    it('computes correct moon position', function() {
+        var date = JulianDate.fromTotalDays(2451545.0, TimeStandard.TAI);
+        var moon = PlanetaryPositions.ComputeMoonPositionInEarthInertialFrame(date);
+        var X = -291632410.61232185;
+        var Y = -266522146.36821631;
+        var Z = -75994518.081043154;
+        expect(X).toEqualEpsilon(moon.x, CesiumMath.EPSILON4);
+        expect(Y).toEqualEpsilon(moon.y, CesiumMath.EPSILON4);
+        expect(Z).toEqualEpsilon(moon.z, CesiumMath.EPSILON4);
+
+        date = JulianDate.fromTotalDays(2456401.5, TimeStandard.TAI);
+        moon = PlanetaryPositions.ComputeMoonPositionInEarthInertialFrame(date);
+        X = -223792974.4736526;
+        Y = 315772435.34490639;
+        Z = 97913011.236112773;
+        expect(X).toEqualEpsilon(moon.x, CesiumMath.EPSILON4);
+        expect(Y).toEqualEpsilon(moon.y, CesiumMath.EPSILON4);
+        expect(Z).toEqualEpsilon(moon.z, CesiumMath.EPSILON4);
+
+        date = JulianDate.fromTotalDays(2455998.591667, TimeStandard.TAI);
+        moon = PlanetaryPositions.ComputeMoonPositionInEarthInertialFrame(date);
+        X = -268426117.00202647;
+        Y = -220468861.73998192;
+        Z = -110670164.58446842;
+        expect(X).toEqualEpsilon(moon.x, CesiumMath.EPSILON4);
+        expect(Y).toEqualEpsilon(moon.y, CesiumMath.EPSILON4);
+        expect(Z).toEqualEpsilon(moon.z, CesiumMath.EPSILON4);
+    });
+
+    it('has the sun rising in the east and setting in the west', function() {
+        //Julian dates for 24 hours, starting from July 6th 2011 @ 01:00 UTC
+        var transformMatrix = new Matrix3();
+        var timesOfDay = [];
+        for ( var i = 1; i < 25; i++) {
+            var date = new Date('July 6, 2011');
+            date.setUTCHours(i, 0, 0, 0);
+            timesOfDay.push(JulianDate.fromDate(date));
+        }
+        var angles = [];
+        for (i = 0; i < 24; i++) {
+            transformMatrix = Transforms.computeIcrfToFixedMatrix(timesOfDay[i], transformMatrix);
+            if (typeof transformMatrix === 'undefined') {
+                transformMatrix = Transforms.computeTemeToPseudoFixedMatrix(timesOfDay[i], transformMatrix);
+            }
+            var position = PlanetaryPositions.ComputeSunPositionInEarthInertialFrame(timesOfDay[i]);
+            transformMatrix.multiplyByVector(position, position);
+            angles.push(CesiumMath.convertLongitudeRange(Math.atan2(position.y, position.x)));
+        }
+        //Expect a clockwise motion.
+        for (i = 1; i < 24; i++) {
+            expect(angles[i]).toBeLessThan(angles[i - 1]);
+        }
+    });
+
+    it('works without a time', function() {
+        PlanetaryPositions.ComputeSunPositionInEarthInertialFrame(undefined);
+    });
+
+});