Skip to content

Commit 5d0dafa

Browse files
parzivailparzivail
committed
improve sun and moon
1 parent d4bbd8c commit 5d0dafa

6 files changed

Lines changed: 538 additions & 129 deletions

File tree

‎SGP.NET.Tests/MoonPositionTests.cs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,7 @@ public void MeeusExample_1992April12_MatchesMeanValues()
6363
var distance = moon.Position.Length;
6464

6565
// Assert mean position matches reference implementation
66-
// Longitude matches PyMeeus geocentric_ecliptical_pos: 133.162655°
66+
// Longitude matches geocentric ecliptical: 133.162655°
6767
Assert.AreEqual(133.163, lambdaDeg, 0.01, "Ecliptic longitude mismatch");
6868
// Latitude is ~-3.23° for mean position (apparent is +13.77° with nutation)
6969
Assert.AreEqual(-3.23, betaDeg, 0.01, "Ecliptic latitude mismatch");

‎SGP.NET/Propagation/Bodies/Moon.cs‎

Lines changed: 208 additions & 123 deletions
Original file line numberDiff line numberDiff line change
@@ -165,129 +165,214 @@ private readonly struct LatitudeTerm(int d, int m, int mp, int f, int sigmaB)
165165
new(4, -1, 0, -1, 115),
166166
new(2, -2, 0, 1, 107),
167167
];
168-
169-
/// <summary>
170-
/// Calculates the Moon's position in Earth-Centered Inertial (ECI) coordinates
171-
/// using the truncated ELP 2000/82 lunar theory.
172-
/// </summary>
173-
/// <param name="time">The time of observation (UTC).</param>
174-
/// <returns>
175-
/// An EciCoordinate representing the Moon's position. The Position vector is in kilometers
176-
/// relative to Earth's center. Velocity is zero (not computed by this algorithm).
177-
/// </returns>
178-
/// <remarks>
179-
/// <para>
180-
/// This implementation follows Jean Meeus, "Astronomical Algorithms", 2nd Edition,
181-
/// Chapter 47 "Position of the Moon" (pp. 337–342). The algorithm is a truncated
182-
/// version of the ELP 2000/82 lunar theory by Chapront-Touzé and Chapront.
183-
/// </para>
184-
/// <para>
185-
/// Accuracy is approximately 10 arcseconds in latitude and 4 arcseconds in longitude
186-
/// for dates between 1900 and 2100. This does not account for nutation or the
187-
/// conversion from Terrestrial Time (TT) to Universal Time (UT1), which introduces
188-
/// an additional error of ~30 arcseconds when using UTC directly.
189-
/// </para>
190-
/// <para>
191-
/// The returned ECI coordinate can be converted to geodetic coordinates via
192-
/// <c>eci.ToGeodetic()</c> to obtain the sublunar point, or used with
193-
/// <c>GroundStation.Observe()</c> for moonrise/moonset and lunar track calculations.
194-
/// </para>
195-
/// </remarks>
196-
public static EciCoordinate Predict(DateTime time)
197-
{
198-
var jde = time.ToJulian();
199-
var t = (jde - 2451545.0) / 36525.0;
200-
201-
var degToRad = Math.PI / 180.0;
202-
203-
var lPrimeDeg = Horner(t, 218.3164477, 481267.88123421, -0.0015786, 1.0 / 538841.0, -1.0 / 65194000.0);
204-
var dDeg = Horner(t, 297.8501921, 445267.1114034, -0.0018819, 1.0 / 545868.0, -1.0 / 113065000.0);
205-
var mDeg = Horner(t, 357.5291092, 35999.0502909, -0.0001535, 1.0 / 24490000.0);
206-
var mPrimeDeg = Horner(t, 134.9633964, 477198.8675055, 0.0087414, 1.0 / 69699.0, -1.0 / 14712000.0);
207-
var fDeg = Horner(t, 93.2720950, 483202.0175233, -0.0036539, -1.0 / 3526000.0, 1.0 / 863310000.0);
208-
209-
var lPrime = MathUtil.Wrap360(lPrimeDeg) * degToRad;
210-
var d = MathUtil.Wrap360(dDeg) * degToRad;
211-
var m = MathUtil.Wrap360(mDeg) * degToRad;
212-
var mPrime = MathUtil.Wrap360(mPrimeDeg) * degToRad;
213-
var f = MathUtil.Wrap360(fDeg) * degToRad;
214-
215-
var a1 = (119.75 + 131.849 * t) * degToRad;
216-
var a2 = (53.09 + 479264.29 * t) * degToRad;
217-
var a3 = (313.45 + 481266.484 * t) * degToRad;
218-
219-
var e = 1.0 - 0.002516 * t - 0.0000074 * t * t;
220-
var e2 = e * e;
221-
222-
var sigmaL = 3958.0 * Math.Sin(a1)
223-
+ 1962.0 * Math.Sin(lPrime - f)
224-
+ 318.0 * Math.Sin(a2);
225-
var sigmaR = 0.0;
226-
var sigmaB = -2235.0 * Math.Sin(lPrime)
227-
+ 382.0 * Math.Sin(a3)
228-
+ 175.0 * Math.Sin(a1 - f)
229-
+ 175.0 * Math.Sin(a1 + f)
230-
+ 127.0 * Math.Sin(lPrime - mPrime)
231-
- 115.0 * Math.Sin(lPrime + mPrime);
232-
233-
for (var i = 0; i < LongitudeTerms.Length; i++)
234-
{
235-
ref var term = ref LongitudeTerms[i];
236-
var arg = d * term.D + m * term.M + mPrime * term.MP + f * term.F;
237-
var sa = Math.Sin(arg);
238-
var ca = Math.Cos(arg);
239-
240-
var factor = term.M switch
241-
{
242-
0 => 1.0,
243-
1 or -1 => e,
244-
2 or -2 => e2,
245-
_ => 1.0
246-
};
247-
248-
sigmaL += term.SigmaL * sa * factor;
249-
sigmaR += term.SigmaR * ca * factor;
250-
}
251-
252-
for (var i = 0; i < LatitudeTerms.Length; i++)
253-
{
254-
ref var term = ref LatitudeTerms[i];
255-
var arg = d * term.D + m * term.M + mPrime * term.MP + f * term.F;
256-
var sb = Math.Sin(arg);
257-
258-
var factor = term.M switch
259-
{
260-
0 => 1.0,
261-
1 or -1 => e,
262-
2 or -2 => e2,
263-
_ => 1.0
264-
};
265-
266-
sigmaB += term.SigmaB * sb * factor;
267-
}
268-
269-
var lambdaRad = WrapPi(lPrime + sigmaL * 1e-6 * degToRad);
270-
var betaRad = sigmaB * 1e-6 * degToRad;
271-
var distanceKm = 385000.56 + sigmaR * 1e-3;
272-
273-
var cosLambda = Math.Cos(lambdaRad);
274-
var sinLambda = Math.Sin(lambdaRad);
275-
var cosBeta = Math.Cos(betaRad);
276-
var sinBeta = Math.Sin(betaRad);
277-
278-
var epsilonDeg = 23.439 - 0.0000004 * (jde - 2451545.0);
279-
var epsilon = epsilonDeg * degToRad;
280-
var cosEpsilon = Math.Cos(epsilon);
281-
var sinEpsilon = Math.Sin(epsilon);
282-
283-
var xKm = distanceKm * cosBeta * cosLambda;
284-
var yKm = distanceKm * (cosEpsilon * cosBeta * sinLambda - sinEpsilon * sinBeta);
285-
var zKm = distanceKm * (sinEpsilon * cosBeta * sinLambda + cosEpsilon * sinBeta);
286-
287-
return new EciCoordinate(time, new Vector3(xKm, yKm, zKm));
288-
}
289-
290-
private static double Horner(double t, params double[] coeffs)
168+
169+
/// <summary>
170+
/// Calculates the Moon's geocentric position using the truncated ELP 2000/82
171+
/// lunar theory as presented by Meeus, "Astronomical Algorithms" 2nd Ed., Chapter 47.
172+
/// </summary>
173+
/// <param name="time">The time of observation (UTC).</param>
174+
/// <returns>
175+
/// An EciCoordinate representing the Moon's position as a geocentric equatorial
176+
/// Cartesian vector in the mean equator and mean equinox of date, in kilometers.
177+
/// Velocity is zero (not computed by this algorithm).
178+
/// </returns>
179+
/// <remarks>
180+
/// <para>
181+
/// The periodic lunar series returns geocentric ecliptic longitude, latitude,
182+
/// and distance referred to the mean equinox/ecliptic of date. This method
183+
/// converts those to a geocentric equatorial Cartesian vector of date.
184+
/// </para>
185+
/// <para>
186+
/// This is not a J2000/GCRF/TEME state vector. Feeding this into observation
187+
/// code that assumes a different ECI frame may produce plausible-looking but
188+
/// incorrect results.
189+
/// </para>
190+
/// <para>
191+
/// Input UTC is converted to TT (Terrestrial Time) using the IERS leap second
192+
/// table before evaluating the lunar series. Using UTC directly would introduce
193+
/// a date-dependent error of tens of arcseconds for modern dates (~69s offset
194+
/// as of 2026, or ~40 arcseconds in lunar position).
195+
/// </para>
196+
/// <para>
197+
/// Nutation, topocentric parallax, light-time, and frame transformations to
198+
/// Earth-fixed coordinates are not included. Meeus quotes approximate accuracy
199+
/// of about 10 arcseconds in longitude and 4 arcseconds in latitude for the
200+
/// truncated lunar series.
201+
/// </para>
202+
/// </remarks>
203+
public static EciCoordinate Predict(DateTime time)
204+
{
205+
// Convert UTC to TT for the ephemeris argument.
206+
// TT = UTC + (TAI - UTC) + 32.184s
207+
var jdUtc = time.ToJulian();
208+
var jdTt = LeapSeconds.UtcToTt(jdUtc);
209+
var t = (jdTt - 2451545.0) / 36525.0;
210+
211+
var degToRad = Math.PI / 180.0;
212+
213+
var lPrimeDeg = Horner(t, 218.3164477, 481267.88123421, -0.0015786, 1.0 / 538841.0, -1.0 / 65194000.0);
214+
var dDeg = Horner(t, 297.8501921, 445267.1114034, -0.0018819, 1.0 / 545868.0, -1.0 / 113065000.0);
215+
var mDeg = Horner(t, 357.5291092, 35999.0502909, -0.0001536, 1.0 / 24490000.0);
216+
var mPrimeDeg = Horner(t, 134.9633964, 477198.8675055, 0.0087414, 1.0 / 69699.0, -1.0 / 14712000.0);
217+
var fDeg = Horner(t, 93.2720950, 483202.0175233, -0.0036539, -1.0 / 3526000.0, 1.0 / 863310000.0);
218+
219+
var lPrime = MathUtil.Wrap360(lPrimeDeg) * degToRad;
220+
var d = MathUtil.Wrap360(dDeg) * degToRad;
221+
var m = MathUtil.Wrap360(mDeg) * degToRad;
222+
var mPrime = MathUtil.Wrap360(mPrimeDeg) * degToRad;
223+
var f = MathUtil.Wrap360(fDeg) * degToRad;
224+
225+
var a1 = (119.75 + 131.849 * t) * degToRad;
226+
var a2 = (53.09 + 479264.29 * t) * degToRad;
227+
var a3 = (313.45 + 481266.484 * t) * degToRad;
228+
229+
var e = 1.0 - 0.002516 * t - 0.0000074 * t * t;
230+
var e2 = e * e;
231+
232+
var sigmaL = 3958.0 * Math.Sin(a1)
233+
+ 1962.0 * Math.Sin(lPrime - f)
234+
+ 318.0 * Math.Sin(a2);
235+
var sigmaR = 0.0;
236+
var sigmaB = -2235.0 * Math.Sin(lPrime)
237+
+ 382.0 * Math.Sin(a3)
238+
+ 175.0 * Math.Sin(a1 - f)
239+
+ 175.0 * Math.Sin(a1 + f)
240+
+ 127.0 * Math.Sin(lPrime - mPrime)
241+
- 115.0 * Math.Sin(lPrime + mPrime);
242+
243+
for (var i = 0; i < LongitudeTerms.Length; i++)
244+
{
245+
ref var term = ref LongitudeTerms[i];
246+
var arg = d * term.D + m * term.M + mPrime * term.MP + f * term.F;
247+
var sa = Math.Sin(arg);
248+
var ca = Math.Cos(arg);
249+
250+
var factor = term.M switch
251+
{
252+
0 => 1.0,
253+
1 or -1 => e,
254+
2 or -2 => e2,
255+
_ => 1.0
256+
};
257+
258+
sigmaL += term.SigmaL * sa * factor;
259+
sigmaR += term.SigmaR * ca * factor;
260+
}
261+
262+
for (var i = 0; i < LatitudeTerms.Length; i++)
263+
{
264+
ref var term = ref LatitudeTerms[i];
265+
var arg = d * term.D + m * term.M + mPrime * term.MP + f * term.F;
266+
var sb = Math.Sin(arg);
267+
268+
var factor = term.M switch
269+
{
270+
0 => 1.0,
271+
1 or -1 => e,
272+
2 or -2 => e2,
273+
_ => 1.0
274+
};
275+
276+
sigmaB += term.SigmaB * sb * factor;
277+
}
278+
279+
// Geocentric ecliptic longitude, latitude, and distance (mean of date).
280+
var lambdaRad = WrapPi(lPrime + sigmaL * 1e-6 * degToRad);
281+
var betaRad = sigmaB * 1e-6 * degToRad;
282+
var distanceKm = 385000.56 + sigmaR * 1e-3;
283+
284+
// Use Meeus mean obliquity polynomial (not the rough linear approximation).
285+
var epsilonRad = MathUtil.DegreesToRadians(Obliquity.MeanObliquityDeg(jdTt));
286+
var cosEpsilon = Math.Cos(epsilonRad);
287+
var sinEpsilon = Math.Sin(epsilonRad);
288+
289+
var cosLambda = Math.Cos(lambdaRad);
290+
var sinLambda = Math.Sin(lambdaRad);
291+
var cosBeta = Math.Cos(betaRad);
292+
var sinBeta = Math.Sin(betaRad);
293+
294+
// Rotate from ecliptic to equatorial (mean of date).
295+
var xKm = distanceKm * cosBeta * cosLambda;
296+
var yKm = distanceKm * (cosEpsilon * cosBeta * sinLambda - sinEpsilon * sinBeta);
297+
var zKm = distanceKm * (sinEpsilon * cosBeta * sinLambda + cosEpsilon * sinBeta);
298+
299+
return new EciCoordinate(time, new Vector3(xKm, yKm, zKm));
300+
}
301+
302+
/// <summary>
303+
/// Calculates the Moon's apparent geocentric position, including nutation corrections.
304+
/// This produces coordinates closer to what would be observed visually or compared
305+
/// against JPL Horizons (excluding light-time and topocentric parallax).
306+
/// </summary>
307+
/// <param name="time">The time of observation (UTC).</param>
308+
/// <returns>
309+
/// An EciCoordinate representing the Moon's apparent position as a geocentric
310+
/// equatorial Cartesian vector in the true equator and true equinox of date,
311+
/// in kilometers. Velocity is zero.
312+
/// </returns>
313+
/// <remarks>
314+
/// <para>
315+
/// This method first computes the mean position via <see cref="Predict" />,
316+
/// then applies nutation in longitude and true obliquity (mean + nutation in
317+
/// obliquity) to produce apparent equatorial coordinates.
318+
/// </para>
319+
/// <para>
320+
/// This is still not a J2000/GCRF/TEME state vector. Light-time correction
321+
/// and topocentric parallax are not included.
322+
/// </para>
323+
/// </remarks>
324+
public static EciCoordinate PredictApparent(DateTime time)
325+
{
326+
var jdUtc = time.ToJulian();
327+
var jdTt = LeapSeconds.UtcToTt(jdUtc);
328+
329+
// Get the mean position first.
330+
var mean = Predict(time);
331+
var pos = mean.Position;
332+
333+
// Convert the mean equatorial vector back to ecliptic using mean obliquity.
334+
var epsilonMeanRad = MathUtil.DegreesToRadians(Obliquity.MeanObliquityDeg(jdTt));
335+
var cosEpsM = Math.Cos(epsilonMeanRad);
336+
var sinEpsM = Math.Sin(epsilonMeanRad);
337+
338+
var x = pos.X;
339+
var y = pos.Y;
340+
var z = pos.Z;
341+
342+
// Equatorial to ecliptic (mean of date).
343+
var xEcl = x;
344+
var yEcl = y * cosEpsM + z * sinEpsM;
345+
var zEcl = -y * sinEpsM + z * cosEpsM;
346+
347+
var rho = Math.Sqrt(xEcl * xEcl + yEcl * yEcl);
348+
var lambdaRad = Math.Atan2(yEcl, xEcl);
349+
var betaRad = Math.Atan2(zEcl, rho);
350+
var distance = pos.Length;
351+
352+
// Apply nutation in longitude.
353+
var deltaPsiRad = MathUtil.DegreesToRadians(Nutation.LongitudeDeg(jdTt));
354+
var lambdaApparent = lambdaRad + deltaPsiRad;
355+
356+
// Use true obliquity (mean + nutation in obliquity) for apparent coordinates.
357+
var deltaEpsilonRad = MathUtil.DegreesToRadians(Nutation.ObliquityDeg(jdTt));
358+
var epsilonTrueRad = epsilonMeanRad + deltaEpsilonRad;
359+
var cosEpsT = Math.Cos(epsilonTrueRad);
360+
var sinEpsT = Math.Sin(epsilonTrueRad);
361+
362+
var cosLambdaA = Math.Cos(lambdaApparent);
363+
var sinLambdaA = Math.Sin(lambdaApparent);
364+
var cosBeta = Math.Cos(betaRad);
365+
var sinBeta = Math.Sin(betaRad);
366+
367+
// Rotate apparent ecliptic to apparent equatorial.
368+
var xApp = distance * cosBeta * cosLambdaA;
369+
var yApp = distance * (cosEpsT * cosBeta * sinLambdaA - sinEpsT * sinBeta);
370+
var zApp = distance * (sinEpsT * cosBeta * sinLambdaA + cosEpsT * sinBeta);
371+
372+
return new EciCoordinate(time, new Vector3(xApp, yApp, zApp));
373+
}
374+
375+
private static double Horner(double t, params double[] coeffs)
291376
{
292377
var result = 0.0;
293378
for (var i = coeffs.Length - 1; i >= 0; i--)

‎SGP.NET/Propagation/Bodies/Sun.cs‎

Lines changed: 13 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,9 @@
44

55
namespace SGPdotNET.Propagation.Bodies;
66

7+
/// <summary>
8+
/// Provides methods to calculate the position of the Sun.
9+
/// </summary>
710
public static class Sun
811
{
912
/// <summary>
@@ -25,6 +28,11 @@ public static class Sun
2528
/// Accuracy is approximately 0.01° (36 arcseconds) for dates between 1950 and 2050.
2629
/// </para>
2730
/// <para>
31+
/// The returned vector is in the mean equator and mean equinox of date frame,
32+
/// not J2000/GCRF/TEME. The mean obliquity of the ecliptic is computed using
33+
/// the Meeus polynomial (Ch. 22) rather than a linear approximation.
34+
/// </para>
35+
/// <para>
2836
/// The returned ECI coordinate can be converted to geodetic coordinates via
2937
/// <c>eci.ToGeodetic()</c> to obtain the subsolar point (latitude/longitude where
3038
/// the Sun is directly overhead), or used with <c>GroundStation.Observe()</c> for
@@ -42,15 +50,15 @@ public static EciCoordinate Predict(DateTime time)
4250
var lambdaDeg = l + 1.915 * Math.Sin(g) + 0.020 * Math.Sin(2.0 * g);
4351
var lambda = MathUtil.DegreesToRadians(lambdaDeg);
4452

45-
var epsilonDeg = 23.439 - 0.0000004 * n;
46-
var epsilon = MathUtil.DegreesToRadians(epsilonDeg);
53+
// Use Meeus mean obliquity polynomial (Ch. 22) instead of the rough linear approximation.
54+
var epsilonRad = MathUtil.DegreesToRadians(Obliquity.MeanObliquityDeg(time.ToJulian()));
55+
var cosEpsilon = Math.Cos(epsilonRad);
56+
var sinEpsilon = Math.Sin(epsilonRad);
4757

4858
var r = 1.00014 - 0.01671 * Math.Cos(g) - 0.00014 * Math.Cos(2.0 * g);
4959

5060
var cosLambda = Math.Cos(lambda);
5161
var sinLambda = Math.Sin(lambda);
52-
var cosEpsilon = Math.Cos(epsilon);
53-
var sinEpsilon = Math.Sin(epsilon);
5462

5563
var xAu = r * cosLambda;
5664
var yAu = r * cosEpsilon * sinLambda;
@@ -62,4 +70,4 @@ public static EciCoordinate Predict(DateTime time)
6270

6371
return new EciCoordinate(time, new Vector3(xKm, yKm, zKm));
6472
}
65-
}
73+
}

0 commit comments

Comments
 (0)