diff --git a/CHANGES.md b/CHANGES.md
index bc887b57aa9b..084711fe8c12 100644
--- a/CHANGES.md
+++ b/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.
diff --git a/Source/Core/Math.js b/Source/Core/Math.js
index da2c8c7a1a88..005d55cf30c6 100644
--- a/Source/Core/Math.js
+++ b/Source/Core/Math.js
@@ -397,7 +397,7 @@ define([
};
/**
- * Alters the value of input x such that -CesiumMath.PI <= x <= CesiumMath.PI
+ * 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 ()-CesiumMath.PI, CesiumMath.PI).
*/
@@ -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 , CesiumMath.TWO_PI).
+ */
+ 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
*/
diff --git a/Source/Core/Simon1994PlanetaryPositions.js b/Source/Core/Simon1994PlanetaryPositions.js
new file mode 100644
index 000000000000..6ff71046f6f0
--- /dev/null
+++ b/Source/Core/Simon1994PlanetaryPositions.js
@@ -0,0 +1,521 @@
+/*global define*/
+define(['./Cartesian3',
+ './DeveloperError',
+ './JulianDate',
+ './Math',
+ './Matrix3',
+ './TimeConstants',
+ './TimeStandard',
+ './Transforms'
+ ],
+ function(
+ Cartesian3,
+ DeveloperError,
+ JulianDate,
+ CesiumMath,
+ Matrix3,
+ TimeConstants,
+ TimeStandard,
+ Transforms) {
+ "use strict";
+
+ /**
+ * Contains functions for finding the Cartesian coordinates of the sun and the moon in the Earth-centered inertial frame.
+ * @exports Simon1994PlanetaryPositions
+ */
+ var Simon1994PlanetaryPositions = {};
+
+ function computeTdbMinusTtSpice(daysSinceJ2000InTerrestrialTime) {
+ /* STK Comments ------------------------------------------------------
+ * This function uses constants designed to be consistent with
+ * the SPICE Toolkit from JPL version N0051 (unitim.c)
+ * M0 = 6.239996
+ * M0Dot = 1.99096871e-7 rad/s = 0.01720197 rad/d
+ * EARTH_ECC = 1.671e-2
+ * TDB_AMPL = 1.657e-3 secs
+ *--------------------------------------------------------------------*/
+
+ //* Values taken as specified in STK Comments except: 0.01720197 rad/day = 1.99096871e-7 rad/sec
+ //* Here we use the more precise value taken from the SPICE value 1.99096871e-7 rad/sec converted to rad/day
+ //* All other constants are consistent with the SPICE implementation of the TDB conversion
+ //* except where we treat the independent time parameter to be in TT instead of TDB.
+ //* This is an approximation made to facilitate performance due to the higher prevalance of
+ //* the TT2TDB conversion over TDB2TT in order to avoid having to iterate when converting to TDB for the JPL ephemeris.
+ //* Days are used instead of seconds to provide a slight improvement in numerical precision.
+
+ //* For more information see:
+ //* http://www.cv.nrao.edu/~rfisher/Ephemerides/times.html#TDB
+ //* ftp://ssd.jpl.nasa.gov/pub/eph/planets/ioms/ExplSupplChap8.pdf
+
+ var g = 6.239996 + (0.0172019696544) * daysSinceJ2000InTerrestrialTime;
+ return 1.657e-3 * Math.sin(g + 1.671e-2 * Math.sin(g));
+ }
+
+ var TdtMinusTai = 32.184;
+ var J2000d = 2451545;
+ function taiToTdb(date, result) {
+ //Converts TAI to TT
+ result = date.addSeconds(TdtMinusTai, result);
+
+ //Converts TT to TDB
+ var days = result.getTotalDays() - J2000d;
+ result = result.addSeconds(computeTdbMinusTtSpice(days), result);
+
+ return result;
+ }
+
+ var epoch = JulianDate.fromTotalDays(2451545.0, TimeStandard.TAI); //Actually TDB (not TAI)
+ var GravitationalParameterOfEarth = 3.98600435e14;
+ var GravitationalParameterOfSun = GravitationalParameterOfEarth * (1.0 + 0.012300034) * 328900.56;
+ var MetersPerKilometer = 1000.0;
+ var RadiansPerDegree = CesiumMath.RADIANS_PER_DEGREE;
+ var RadiansPerArcSecond = CesiumMath.RADIANS_PER_ARCSECOND;
+ var MetersPerAstronomicalUnit = 1.49597870e+11; // IAU 1976 value
+
+ var perifocalToEquatorial = new Matrix3();
+ function elementsToCartesian(semimajorAxis, eccentricity, inclination, longitudeOfPerigee, longitudeOfNode, meanLongitude, gravitationalParameter, result) {
+ if (inclination < 0.0) {
+ inclination = -inclination;
+ longitudeOfNode += CesiumMath.PI;
+ }
+ if (inclination < 0 || inclination > CesiumMath.PI) {
+ throw new DeveloperError('The inclination is out of range. Inclination must be greater than or equal to zero and less than or equal to Pi radians.');
+ }
+
+ var radiusOfPeriapsis = semimajorAxis * (1.0 - eccentricity);
+ var argumentOfPeriapsis = longitudeOfPerigee - longitudeOfNode;
+ var rightAscensionOfAscendingNode = longitudeOfNode;
+ var trueAnomaly = meanAnomalyToTrueAnomaly(meanLongitude - longitudeOfPerigee, eccentricity);
+ var type = chooseOrbit(eccentricity, 0.0);
+ if (type === 'Hyperbolic' && Math.abs(CesiumMath.NegativePiToPi(trueAnomaly)) >= Math.acos(- 1.0 / eccentricity)) {
+ throw new DeveloperError('The true anomaly of the hyperbolic orbit lies outside of the bounds of the hyperbola.');
+ }
+ perifocalToCartesianMatrix(argumentOfPeriapsis, inclination, rightAscensionOfAscendingNode, perifocalToEquatorial);
+ var semilatus = radiusOfPeriapsis * (1.0 + eccentricity);
+ var costheta = Math.cos(trueAnomaly);
+ var sintheta = Math.sin(trueAnomaly);
+
+ var denom = (1.0 + eccentricity * costheta);
+ if (denom <= CesiumMath.Epsilon10) {
+ throw new DeveloperError('elements cannot be converted to cartesian');
+ }
+
+ var radius = semilatus / denom;
+ if (typeof result === 'undefined') {
+ result = new Cartesian3(radius * costheta, radius * sintheta, 0.0);
+ } else {
+ result.x = radius * costheta;
+ result.y = radius * sintheta;
+ result.z = 0.0;
+ }
+
+ return perifocalToEquatorial.multiplyByVector(result, result);
+ }
+
+ 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';
+ }
+ }
+
+ // Calculates the true anomaly given the mean anomaly and the eccentricity.
+ function meanAnomalyToTrueAnomaly(meanAnomaly, eccentricity) {
+ if (eccentricity < 0.0 || eccentricity >= 1.0) {
+ throw new DeveloperError('eccentricity out of range.');
+ }
+ var eccentricAnomaly = meanAnomalyToEccentricAnomaly(meanAnomaly, eccentricity);
+ return eccentricAnomalyToTrueAnomaly(eccentricAnomaly, eccentricity);
+ }
+
+ var maxIterationCount = 50;
+ var keplerEqConvergence = CesiumMath.EPSILON8;
+ // Calculates the eccentric anomaly given the mean anomaly and the eccentricity.
+ function meanAnomalyToEccentricAnomaly(meanAnomaly, eccentricity) {
+ if (eccentricity < 0.0 || eccentricity >= 1.0) {
+ throw new DeveloperError('eccentricity out of range.');
+ }
+ 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');
+ // STK Components uses a numerical method to find the eccentric anomaly in the case that Kepler's
+ // equation does not converge. We don't expect that to ever be necessary for the reasonable orbits used here.
+ }
+
+ eccentricAnomaly = iterationValue + revs * CesiumMath.TWO_PI;
+ return eccentricAnomaly;
+ }
+
+ // Calculates the true anomaly given the eccentric anomaly and the eccentricity.
+ function eccentricAnomalyToTrueAnomaly(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.
+ function perifocalToCartesianMatrix(argumentOfPeriapsis, inclination, rightAscension, result) {
+ 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);
+ if (typeof result === 'undefined') {
+ result = 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);
+ } else {
+ result[0] = cosraan * cosap - sinraan * sinap * cosi;
+ result[1] = sinraan * cosap + cosraan * sinap * cosi;
+ result[2] = sinap * sini;
+ result[3] = -cosraan * sinap - sinraan * cosap * cosi;
+ result[4] = -sinraan * sinap + cosraan * cosap * cosi;
+ result[5] = cosap * sini;
+ result[6] = sinraan * sini;
+ result[7] = -cosraan * sini;
+ result[8] = cosi;
+ }
+ return result;
+ }
+
+ // From section 5.8
+ var semiMajorAxis0 = 1.0000010178 * MetersPerAstronomicalUnit;
+ var meanLongitude0 = 100.46645683 * RadiansPerDegree;
+ var meanLongitude1 = 1295977422.83429 * RadiansPerArcSecond;
+
+ // From table 6
+ var p1u = 16002;
+ var p2u = 21863;
+ var p3u = 32004;
+ var p4u = 10931;
+ var p5u = 14529;
+ var p6u = 16368;
+ var p7u = 15318;
+ var p8u = 32794;
+
+ var Ca1 = 64 * 1e-7 * MetersPerAstronomicalUnit;
+ var Ca2 = -152 * 1e-7 * MetersPerAstronomicalUnit;
+ var Ca3 = 62 * 1e-7 * MetersPerAstronomicalUnit;
+ var Ca4 = -8 * 1e-7 * MetersPerAstronomicalUnit;
+ var Ca5 = 32 * 1e-7 * MetersPerAstronomicalUnit;
+ var Ca6 = -41 * 1e-7 * MetersPerAstronomicalUnit;
+ var Ca7 = 19 * 1e-7 * MetersPerAstronomicalUnit;
+ var Ca8 = -11 * 1e-7 * MetersPerAstronomicalUnit;
+
+ var Sa1 = -150 * 1e-7 * MetersPerAstronomicalUnit;
+ var Sa2 = -46 * 1e-7 * MetersPerAstronomicalUnit;
+ var Sa3 = 68 * 1e-7 * MetersPerAstronomicalUnit;
+ var Sa4 = 54 * 1e-7 * MetersPerAstronomicalUnit;
+ var Sa5 = 14 * 1e-7 * MetersPerAstronomicalUnit;
+ var Sa6 = 24 * 1e-7 * MetersPerAstronomicalUnit;
+ var Sa7 = -28 * 1e-7 * MetersPerAstronomicalUnit;
+ var Sa8 = 22 * 1e-7 * MetersPerAstronomicalUnit;
+
+ var q1u = 10;
+ var q2u = 16002;
+ var q3u = 21863;
+ var q4u = 10931;
+ var q5u = 1473;
+ var q6u = 32004;
+ var q7u = 4387;
+ var q8u = 73;
+
+ var Cl1 = -325 * 1e-7;
+ var Cl2 = -322 * 1e-7;
+ var Cl3 = -79 * 1e-7;
+ var Cl4 = 232 * 1e-7;
+ var Cl5 = -52 * 1e-7;
+ var Cl6 = 97 * 1e-7;
+ var Cl7 = 55 * 1e-7;
+ var Cl8 = -41 * 1e-7;
+
+ var Sl1 = -105 * 1e-7;
+ var Sl2 = -137 * 1e-7;
+ var Sl3 = 258 * 1e-7;
+ var Sl4 = 35 * 1e-7;
+ var Sl5 = -116 * 1e-7;
+ var Sl6 = -88 * 1e-7;
+ var Sl7 = -112 * 1e-7;
+ var Sl8 = -80 * 1e-7;
+
+ var scratchDate = new JulianDate();
+ /**
+ * Gets a point describing the motion of the Earth-Moon barycenter according to the equations
+ * described in section 6.
+ */
+
+ function computeSimonEarthMoonBarycenter(date, result) {
+
+ // t is thousands of years from J2000 TDB
+ taiToTdb(date, scratchDate);
+ var x = (scratchDate.getJulianDayNumber() - epoch.getJulianDayNumber()) + ((scratchDate.getSecondsOfDay() - epoch.getSecondsOfDay())/TimeConstants.SECONDS_PER_DAY);
+ var t = x / (TimeConstants.DAYS_PER_JULIAN_CENTURY * 10.0);
+
+ var u = 0.35953620 * t;
+ var semimajorAxis = semiMajorAxis0 +
+ Ca1 * Math.cos(p1u * u) + Sa1 * Math.sin(p1u * u) +
+ Ca2 * Math.cos(p2u * u) + Sa2 * Math.sin(p2u * u) +
+ Ca3 * Math.cos(p3u * u) + Sa3 * Math.sin(p3u * u) +
+ Ca4 * Math.cos(p4u * u) + Sa4 * Math.sin(p4u * u) +
+ Ca5 * Math.cos(p5u * u) + Sa5 * Math.sin(p5u * u) +
+ Ca6 * Math.cos(p6u * u) + Sa6 * Math.sin(p6u * u) +
+ Ca7 * Math.cos(p7u * u) + Sa7 * Math.sin(p7u * u) +
+ Ca8 * Math.cos(p8u * u) + Sa8 * Math.sin(p8u * u);
+ var meanLongitude = meanLongitude0 + meanLongitude1 * t +
+ Cl1 * Math.cos(q1u * u) + Sl1 * Math.sin(q1u * u) +
+ Cl2 * Math.cos(q2u * u) + Sl2 * Math.sin(q2u * u) +
+ Cl3 * Math.cos(q3u * u) + Sl3 * Math.sin(q3u * u) +
+ Cl4 * Math.cos(q4u * u) + Sl4 * Math.sin(q4u * u) +
+ Cl5 * Math.cos(q5u * u) + Sl5 * Math.sin(q5u * u) +
+ Cl6 * Math.cos(q6u * u) + Sl6 * Math.sin(q6u * u) +
+ Cl7 * Math.cos(q7u * u) + Sl7 * Math.sin(q7u * u) +
+ Cl8 * Math.cos(q8u * u) + Sl8 * Math.sin(q8u * u);
+
+ // All constants in this part are from section 5.8
+ var eccentricity = 0.0167086342 - 0.0004203654 * t;
+ var longitudeOfPerigee = 102.93734808 * RadiansPerDegree + 11612.35290 * RadiansPerArcSecond * t;
+ var inclination = 469.97289 * RadiansPerArcSecond * t;
+ var longitudeOfNode = 174.87317577 * RadiansPerDegree - 8679.27034 * RadiansPerArcSecond * t;
+
+ return elementsToCartesian(semimajorAxis, eccentricity, inclination, longitudeOfPerigee,
+ longitudeOfNode, meanLongitude, GravitationalParameterOfSun, result);
+ }
+
+ /**
+ * Gets a point describing the position of the moon according to the equations described in section 4.
+ */
+ function computeSimonMoon(date, result) {
+ taiToTdb(date, scratchDate);
+ var x = (scratchDate.getJulianDayNumber() - epoch.getJulianDayNumber()) + ((scratchDate.getSecondsOfDay() - epoch.getSecondsOfDay())/TimeConstants.SECONDS_PER_DAY);
+ var t = x / (TimeConstants.DAYS_PER_JULIAN_CENTURY);
+ var t2 = t * t;
+ var t3 = t2 * t;
+ var t4 = t3 * t;
+
+ // Terms from section 3.4 (b.1)
+ var semimajorAxis = 383397.7725 + 0.0040 * t;
+ var eccentricity = 0.055545526 - 0.000000016 * t;
+ var inclinationConstant = 5.15668983 * RadiansPerDegree;
+ var inclinationSecPart = -0.00008 * t + 0.02966 * t2 -
+ 0.000042 * t3 - 0.00000013 * t4;
+ var longitudeOfPerigeeConstant = 83.35324312 * RadiansPerDegree;
+ var longitudeOfPerigeeSecPart = 14643420.2669 * t - 38.2702 * t2 -
+ 0.045047 * t3 + 0.00021301 * t4;
+ var longitudeOfNodeConstant = 125.04455501 * RadiansPerDegree;
+ var longitudeOfNodeSecPart = -6967919.3631 * t + 6.3602 * t2 +
+ 0.007625 * t3 - 0.00003586 * t4;
+ var meanLongitudeConstant = 218.31664563 * RadiansPerDegree;
+ var meanLongitudeSecPart = 1732559343.48470 * t - 6.3910 * t2 +
+ 0.006588 * t3 - 0.00003169 * t4;
+
+ // Delaunay arguments from section 3.5 b
+ var D = 297.85019547 * RadiansPerDegree + RadiansPerArcSecond *
+ (1602961601.2090 * t - 6.3706 * t2 + 0.006593 * t3 - 0.00003169 * t4);
+ var F = 93.27209062 * RadiansPerDegree + RadiansPerArcSecond *
+ (1739527262.8478 * t - 12.7512 * t2 - 0.001037 * t3 + 0.00000417 * t4);
+ var l = 134.96340251 * RadiansPerDegree + RadiansPerArcSecond *
+ (1717915923.2178 * t + 31.8792 * t2 + 0.051635 * t3 - 0.00024470 * t4);
+ var lprime = 357.52910918 * RadiansPerDegree + RadiansPerArcSecond *
+ (129596581.0481 * t - 0.5532 * t2 + 0.000136 * t3 - 0.00001149 * t4);
+ var psi = 310.17137918 * RadiansPerDegree - RadiansPerArcSecond *
+ (6967051.4360 * t + 6.2068 * t2 + 0.007618 * t3 - 0.00003219 * t4);
+
+ // Add terms from Table 4
+ var twoD = 2.0 * D;
+ var fourD = 4.0 * D;
+ var sixD = 6.0 * D;
+ var twol = 2.0 * l;
+ var threel = 3.0 * l;
+ var fourl = 4.0 * l;
+ var twoF = 2.0 * F;
+ semimajorAxis += 3400.4 * Math.cos(twoD) - 635.6 * Math.cos(twoD - l) -
+ 235.6 * Math.cos(l) + 218.1 * Math.cos(twoD - lprime) +
+ 181.0 * Math.cos(twoD + l);
+ eccentricity += 0.014216 * Math.cos(twoD - l) + 0.008551 * Math.cos(twoD - twol) -
+ 0.001383 * Math.cos(l) + 0.001356 * Math.cos(twoD + l) -
+ 0.001147 * Math.cos(fourD - threel) - 0.000914 * Math.cos(fourD - twol) +
+ 0.000869 * Math.cos(twoD - lprime - l) - 0.000627 * Math.cos(twoD) -
+ 0.000394 * Math.cos(fourD - fourl) + 0.000282 * Math.cos(twoD - lprime - twol) -
+ 0.000279 * Math.cos(D - l) - 0.000236 * Math.cos(twol) +
+ 0.000231 * Math.cos(fourD) + 0.000229 * Math.cos(sixD - fourl) -
+ 0.000201 * Math.cos(twol - twoF);
+ inclinationSecPart += 486.26 * Math.cos(twoD - twoF) - 40.13 * Math.cos(twoD) +
+ 37.51 * Math.cos(twoF) + 25.73 * Math.cos(twol - twoF) +
+ 19.97 * Math.cos(twoD - lprime - twoF);
+ longitudeOfPerigeeSecPart += -55609 * Math.sin(twoD - l) - 34711 * Math.sin(twoD - twol) -
+ 9792 * Math.sin(l) + 9385 * Math.sin(fourD - threel) +
+ 7505 * Math.sin(fourD - twol) + 5318 * Math.sin(twoD + l) +
+ 3484 * Math.sin(fourD - fourl) - 3417 * Math.sin(twoD - lprime - l) -
+ 2530 * Math.sin(sixD - fourl) - 2376 * Math.sin(twoD) -
+ 2075 * Math.sin(twoD - threel) - 1883 * Math.sin(twol) -
+ 1736 * Math.sin(sixD - 5.0 * l) + 1626 * Math.sin(lprime) -
+ 1370 * Math.sin(sixD - threel);
+ longitudeOfNodeSecPart += -5392 * Math.sin(twoD - twoF) - 540 * Math.sin(lprime) -
+ 441 * Math.sin(twoD) + 423 * Math.sin(twoF) -
+ 288 * Math.sin(twol - twoF);
+ meanLongitudeSecPart += -3332.9 * Math.sin(twoD) + 1197.4 * Math.sin(twoD - l) -
+ 662.5 * Math.sin(lprime) + 396.3 * Math.sin(l) -
+ 218.0 * Math.sin(twoD - lprime);
+
+ // Add terms from Table 5
+ var twoPsi = 2.0 * psi;
+ var threePsi = 3.0 * psi;
+ inclinationSecPart += 46.997 * Math.cos(psi) * t - 0.614 * Math.cos(twoD - twoF + psi) * t +
+ 0.614 * Math.cos(twoD - twoF - psi) * t - 0.0297 * Math.cos(twoPsi) * t2 -
+ 0.0335 * Math.cos(psi) * t2 + 0.0012 * Math.cos(twoD - twoF + twoPsi) * t2 -
+ 0.00016 * Math.cos(psi) * t3 + 0.00004 * Math.cos(threePsi) * t3 +
+ 0.00004 * Math.cos(twoPsi) * t3;
+ var perigeeAndMean = 2.116 * Math.sin(psi) * t - 0.111 * Math.sin(twoD - twoF - psi) * t -
+ 0.0015 * Math.sin(psi) * t2;
+ longitudeOfPerigeeSecPart += perigeeAndMean;
+ meanLongitudeSecPart += perigeeAndMean;
+ longitudeOfNodeSecPart += -520.77 * Math.sin(psi) * t + 13.66 * Math.sin(twoD - twoF + psi) * t +
+ 1.12 * Math.sin(twoD - psi) * t - 1.06 * Math.sin(twoF - psi) * t +
+ 0.660 * Math.sin(twoPsi) * t2 + 0.371 * Math.sin(psi) * t2 -
+ 0.035 * Math.sin(twoD - twoF + twoPsi) * t2 - 0.015 * Math.sin(twoD - twoF + psi) * t2 +
+ 0.0014 * Math.sin(psi) * t3 - 0.0011 * Math.sin(threePsi) * t3 -
+ 0.0009 * Math.sin(twoPsi) * t3;
+
+ // Add constants and convert units
+ semimajorAxis *= MetersPerKilometer;
+ var inclination = inclinationConstant + inclinationSecPart * RadiansPerArcSecond;
+ var longitudeOfPerigee = longitudeOfPerigeeConstant + longitudeOfPerigeeSecPart * RadiansPerArcSecond;
+ var meanLongitude = meanLongitudeConstant + meanLongitudeSecPart * RadiansPerArcSecond;
+ var longitudeOfNode = longitudeOfNodeConstant + longitudeOfNodeSecPart * RadiansPerArcSecond;
+
+ return elementsToCartesian(semimajorAxis, eccentricity, inclination, longitudeOfPerigee,
+ longitudeOfNode, meanLongitude, GravitationalParameterOfEarth, result);
+ }
+
+ /**
+ * Gets a point describing the motion of the Earth. This point uses the Moon point and
+ * the 1992 mu value (ratio between Moon and Earth masses) in Table 2 of the paper in order
+ * to determine the position of the Earth relative to the Earth-Moon barycenter.
+ */
+ var moonEarthMassRatio = 0.012300034; // From 1992 mu value in Table 2
+ var factor = moonEarthMassRatio / (moonEarthMassRatio + 1.0) * -1;
+ function computeSimonEarth(date, result) {
+ var moon = computeSimonMoon(date);
+ result = moon.multiplyByScalar(factor, result);
+ return result;
+ }
+
+ // Values for the axesTransformation needed for the rotation were found using the STK Components
+ // GreographicTransformer on the position of the sun center of mass point and the earth J2000 frame.
+
+ var axesTransformation = new Matrix3(1.0000000000000002, 5.619723173785822e-16, 4.690511510146299e-19,
+ -5.154129427414611e-16, 0.9174820620691819, -0.39777715593191376,
+ -2.23970096136568e-16, 0.39777715593191376, 0.9174820620691819);
+ var translation = new Cartesian3();
+ /**
+ * Computes the position of the Sun in the Earth-centered inertial frame
+ *
+ * @param {JulianDate} [julianDate] The time at which to compute the Sun's position, if not provided the current system time is used.
+ * @param {Cartesian3} [result] The object onto which to store the result.
+ * @returns {Cartesian3} Calculated sun position
+ */
+ Simon1994PlanetaryPositions.ComputeSunPositionInEarthInertialFrame= function(date, result){
+ if (typeof date === 'undefined') {
+ date = new JulianDate();
+ }
+ //first forward transformation
+ translation = computeSimonEarthMoonBarycenter(date, translation);
+ result = translation.negate(result);
+
+ //second forward transformation
+ computeSimonEarth(date, translation);
+
+ result.subtract(translation, result);
+ axesTransformation.multiplyByVector(result, result);
+
+ return result;
+ };
+
+ /**
+ * Computes the position of the Moon in the Earth-centered inertial frame
+ *
+ * @param {JulianDate} [julianDate] The time at which to compute the Sun's position, if not provided the current system time is used.
+ * @param {Cartesian3} [result] The object onto which to store the result.
+ * @returns {Cartesian3} Calculated moon position
+ */
+ Simon1994PlanetaryPositions.ComputeMoonPositionInEarthInertialFrame = function(date, result){
+ if (typeof date === 'undefined') {
+ date = new JulianDate();
+ }
+ result = computeSimonMoon(date, result);
+ axesTransformation.multiplyByVector(result, result);
+
+ return result;
+ };
+
+ return Simon1994PlanetaryPositions;
+});
\ No newline at end of file
diff --git a/Source/Core/Transforms.js b/Source/Core/Transforms.js
index f7f684feecd6..4a4a9879df5c 100644
--- a/Source/Core/Transforms.js
+++ b/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);
diff --git a/Source/Core/computeSunPosition.js b/Source/Core/computeSunPosition.js
deleted file mode 100644
index 44be1e0fc3b3..000000000000
--- a/Source/Core/computeSunPosition.js
+++ /dev/null
@@ -1,116 +0,0 @@
-/*global define*/
-define([
- './isLeapYear',
- './DeveloperError',
- './Math',
- './Cartesian3',
- './Cartographic',
- './Matrix3',
- './JulianDate'
- ],
- function(
- isLeapYear,
- DeveloperError,
- CesiumMath,
- Cartesian3,
- Cartographic,
- Matrix3,
- JulianDate) {
- "use strict";
-
- var offSets = [0, // January
- 31, // February
- 59, // March
- 90, // April
- 120,// May
- 151,// June
- 181,// July
- 212,// August
- 243,// September
- 273,// October
- 304,// November
- 334 // December
- ];
-
- var scratch = new Cartesian3();
- var transform = new Matrix3(0.0, 0.0, 0.0,
- 0.0, 1.0, 0.0,
- 0.0, 0.0, 0.0);
-
- var AU_TO_METERS = 149597870700.0;
-
- var direction = new Cartesian3();
-
- /**
- * Computes the position of the Sun in Earth's fixed frame.
- * @exports computeSunPosition
- *
- * @param {JulianDate} [julianDate] The time at which to compute the Sun's position, if not provided the current system time is used.
- * @param {Cartesian3} [result] The object onto which to store the result.
- * @return {Cartesian3} The modified result parameter or a new Cartesian3 instance if none was provided.
- *
- * @exception {DeveloperError} julianDate is required.
- *
- * @example
- * var sunPosition = computeSunPosition(new JulianDate());
- */
- var computeSunPosition = function(julianDate, result) {
- if (typeof julianDate === 'undefined') {
- julianDate = new JulianDate();
- }
-
- var T = (julianDate.getTotalDays() - 2451545.0) / 36525;
- var meanAnomaly = CesiumMath.convertLongitudeRange(CesiumMath.toRadians(357.5277233 + 35999.05034 * T));
- var distanceToSunInAU = 1.000140612 - 0.016708617 * Math.cos(meanAnomaly) - 0.000139589 * Math.cos(2 * meanAnomaly);
-
- var date = julianDate.toDate();
- var month = date.getUTCMonth();
- var dayOfYear = date.getUTCDate() + offSets[month];
-
- if (isLeapYear(date.getUTCFullYear()) && month > 1) {
- dayOfYear++;
- }
-
- var temp = CesiumMath.toRadians((360.0 / 365.0) * (dayOfYear - 81.0));
- var equationOfTime = 9.87 * Math.sin(2 * temp) - 7.53 * Math.cos(temp) - 1.5 * Math.sin(temp);
- var timeFraction = julianDate.getJulianTimeFraction();
- var localTime;
- if (timeFraction >= 0.5) {
- localTime = timeFraction * 24.0 - 12.0;
- } else {
- localTime = 12.0 + timeFraction * 24.0;
- }
- var localSolarTime = localTime + (equationOfTime / 60.0);
- var hourAngle = CesiumMath.toRadians(15.0 * (12.0 - localSolarTime));
- var declinationAngle = Math.asin(0.39795 * Math.cos(CesiumMath.toRadians(0.98563 * (dayOfYear - 173.0))));
-
- var latitudeAngle = 0.0;
-
- var cosLatitudeAngle = Math.cos(latitudeAngle);
- var sinLatitudeAngle = Math.sin(latitudeAngle);
-
- //Since some of these are constant,
- //there's no need to set them every time.
- transform[0] = cosLatitudeAngle;
- //transform[1] = 0.0;
- transform[2] = -sinLatitudeAngle;
- //transform[3] = 0.0;
- //transform[4] = 1.0;
- //transform[5] = 0.0;
- transform[6] = sinLatitudeAngle;
- //transform[7] = 0.0;
- transform[8] = cosLatitudeAngle;
-
- var cosDeclinationAngle = Math.cos(declinationAngle);
- scratch.x = cosDeclinationAngle * Math.cos(hourAngle);
- scratch.y = cosDeclinationAngle * Math.sin(hourAngle);
- scratch.z = Math.sin(declinationAngle);
-
- Matrix3.multiplyByVector(transform, scratch, direction);
- var distance = distanceToSunInAU * AU_TO_METERS;
-
- return direction.multiplyByScalar(distance, result);
- };
-
- return computeSunPosition;
-});
diff --git a/Source/Renderer/UniformState.js b/Source/Renderer/UniformState.js
index 7f87635968bd..0dd811a761c8 100644
--- a/Source/Renderer/UniformState.js
+++ b/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);
}
/**
diff --git a/Specs/Core/Simon1994PlanetaryPositionsSpec.js b/Specs/Core/Simon1994PlanetaryPositionsSpec.js
new file mode 100644
index 000000000000..4d2b4740960f
--- /dev/null
+++ b/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);
+ });
+
+});
diff --git a/Specs/Core/computeSunPositionSpec.js b/Specs/Core/computeSunPositionSpec.js
deleted file mode 100644
index d23853e13835..000000000000
--- a/Specs/Core/computeSunPositionSpec.js
+++ /dev/null
@@ -1,43 +0,0 @@
-/*global defineSuite*/
-defineSuite([
- 'Core/computeSunPosition',
- 'Core/Cartesian3',
- 'Core/JulianDate',
- 'Core/Math'
- ], function(
- computeSunPosition,
- Cartesian3,
- JulianDate,
- CesiumMath) {
- "use strict";
- /*global jasmine,describe,xdescribe,it,xit,expect,beforeEach,afterEach,beforeAll,afterAll,spyOn,runs,waits,waitsFor*/
-
- 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 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++) {
- var position = computeSunPosition(timesOfDay[i]);
- 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 with a result parameter', function() {
- var result = new Cartesian3();
- var returnedResult = computeSunPosition(new JulianDate(), result);
- expect(result).toBe(returnedResult);
- });
-
- it('works without a time', function() {
- computeSunPosition(undefined);
- });
-});
\ No newline at end of file