diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index dd86e5ec..e449156f 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -221,6 +221,9 @@ class Vector: def __sub__(self, other): return Vector(self.x - other.x, self.y - other.y, self.z - other.z, self.t) + def __neg__(self): + return Vector(-self.x, -self.y, -self.z, self.t) + def format(self, coord_format): """Returns a custom format string representation of the vector.""" layout = '({:' + coord_format + '}, {:' + coord_format + '}, {:' + coord_format + '}, {})' @@ -8506,9 +8509,12 @@ def _CalcShadow(body_radius_km, time, target, sdir): def _EarthShadow(time): - e = _CalcEarth(time) + # This function helps find when the Earth's shadow falls upon the Moon. + # Light-travel and aberration corrected vector from the Earth to the Sun. + # The negative vector -s is thus the path of sunlight through the center of the Earth. + s = GeoVector(Body.Sun, time, True) m = GeoMoon(time) - return _CalcShadow(_EARTH_ECLIPSE_RADIUS_KM, time, m, e) + return _CalcShadow(_EARTH_ECLIPSE_RADIUS_KM, time, m, -s) def _MoonShadow(time): @@ -8628,6 +8634,12 @@ class LunarEclipseInfo: The `kind` field thus holds one of the values `EclipseKind.Penumbral`, `EclipseKind.Partial`, or `EclipseKind.Total`, depending on the kind of lunar eclipse found. + The `obscuration` field holds a value in the range [0, 1] that indicates what fraction + of the Moon's apparent disc area is covered by the Earth's umbra at the eclipse's peak. + This indicates how dark the peak eclipse appears. For penumbral eclipses, the obscuration + is 0, because the Moon does not pass through the Earth's umbra. For partial eclipses, + the obscuration is somewhere between 0 and 1. For total lunar eclipses, the obscuration is 1. + Field `peak` holds the date and time of the peak of the eclipse, when it is at its peak. Fields `sd_penum`, `sd_partial`, and `sd_total` hold the semi-duration of each phase @@ -8640,6 +8652,8 @@ class LunarEclipseInfo: ---------- kind : EclipseKind The type of lunar eclipse found. + obscuration : float + The peak fraction of the Moon's apparent disc that is covered by the Earth's umbra. peak : Time The time of the eclipse at its peak. sd_penum : float @@ -8649,16 +8663,18 @@ class LunarEclipseInfo: sd_total : float The semi-duration of the penumbral phase in minutes, or 0.0 if none. """ - def __init__(self, kind, peak, sd_penum, sd_partial, sd_total): + def __init__(self, kind, obscuration, peak, sd_penum, sd_partial, sd_total): self.kind = kind + self.obscuration = obscuration self.peak = peak self.sd_penum = sd_penum self.sd_partial = sd_partial self.sd_total = sd_total def __repr__(self): - return 'LunarEclipseInfo({}, peak={}, sd_penum={}, sd_partial={}, sd_total={})'.format( + return 'LunarEclipseInfo({}, obscuration={}, peak={}, sd_penum={}, sd_partial={}, sd_total={})'.format( self.kind, + self.obscuration, repr(self.peak), self.sd_penum, self.sd_partial, @@ -8989,6 +9005,47 @@ def _MoonEclipticLatitudeDegrees(time): moon = _CalcMoon(time) return math.degrees(moon.geo_eclip_lat) +def _Obscuration(a, b, c): + # a = radius of first disc + # b = radius of second disc + # c = distance between the centers of the discs + if a <= 0.0: + raise Error('Radius of first disc must be positive.') + + if b <= 0.0: + raise Error('Radius of second disc must be positive.') + + if c < 0.0: + raise Error('Distance between discs is not allowed to be negative.') + + if c >= a + b: + # The discs are too far apart to have any overlapping area + return 0.0 + + if c == 0.0: + # The discs have a common center. Therefore, one disc is inside the other. + return 1.0 if (a <= b) else (b*b)/(a*a) + + x = (a*a - b*b + c*c) / (2.0*c) + radicand = a*a - x*x + if radicand <= 0.0: + # The circumferences do not intersect, or are tangent. + # We already ruled out the case of non-overlapping discs. + # Therefore, one disc is inside the other. + return 1.0 if (a <= b) else (b*b)/(a*a) + + # The discs overlap fractionally in a pair of lens-shaped areas. + y = math.sqrt(radicand) + + # Return the overlapping fractional area. + # There are two lens-shaped areas, one to the left of x, the other to the right of x. + # Each part is calculated by subtracting a triangular area from a sector's area. + lens1 = a*a*math.acos(x/a) - x*y + lens2 = b*b*math.acos((c-x)/b) - (c-x)*y + + # Find the fractional area with respect to the first disc. + return (lens1 + lens2) / (math.pi*a*a) + def SearchLunarEclipse(startTime): """Searches for a lunar eclipse. @@ -9028,6 +9085,7 @@ def SearchLunarEclipse(startTime): if shadow.r < shadow.p + _MOON_MEAN_RADIUS_KM: # This is at least a penumbral eclipse. We will return a result. kind = EclipseKind.Penumbral + obscuration = 0.0 sd_total = 0.0 sd_partial = 0.0 sd_penum = _ShadowSemiDurationMinutes(shadow.time, shadow.p + _MOON_MEAN_RADIUS_KM, 200.0) @@ -9040,9 +9098,12 @@ def SearchLunarEclipse(startTime): if shadow.r + _MOON_MEAN_RADIUS_KM < shadow.k: # This is a total eclipse. kind = EclipseKind.Total + obscuration = 1.0 sd_total = _ShadowSemiDurationMinutes(shadow.time, shadow.k - _MOON_MEAN_RADIUS_KM, sd_partial) + else: + obscuration = _Obscuration(_MOON_MEAN_RADIUS_KM, shadow.k, shadow.r) - return LunarEclipseInfo(kind, shadow.time, sd_penum, sd_partial, sd_total) + return LunarEclipseInfo(kind, obscuration, shadow.time, sd_penum, sd_partial, sd_total) # We didn't find an eclipse on this full moon, so search for the next one. fmtime = fullmoon.AddDays(10) diff --git a/demo/python/correct/lunar_eclipse.txt b/demo/python/correct/lunar_eclipse.txt index aa3f5f1f..f50bf7eb 100644 --- a/demo/python/correct/lunar_eclipse.txt +++ b/demo/python/correct/lunar_eclipse.txt @@ -1,50 +1,50 @@ -1988-03-03T16:04:42.484Z Partial eclipse begins. -1988-03-03T16:13:29.227Z Peak of partial eclipse. -1988-03-03T16:22:15.970Z Partial eclipse ends. +1988-03-03T16:05:01.550Z Partial eclipse begins. +1988-03-03T16:12:44.158Z Peak of partial eclipse. +1988-03-03T16:20:26.766Z Partial eclipse ends. -1988-08-27T10:07:57.338Z Partial eclipse begins. -1988-08-27T11:05:04.753Z Peak of partial eclipse. -1988-08-27T12:02:12.168Z Partial eclipse ends. +1988-08-27T10:07:28.658Z Partial eclipse begins. +1988-08-27T11:04:30.915Z Peak of partial eclipse. +1988-08-27T12:01:33.172Z Partial eclipse ends. -1989-02-20T13:44:07.827Z Partial eclipse begins. -1989-02-20T14:56:15.430Z Total eclipse begins. -1989-02-20T15:36:04.478Z Peak of total eclipse. -1989-02-20T16:15:53.526Z Total eclipse ends. -1989-02-20T17:28:01.129Z Partial eclipse ends. +1989-02-20T13:43:24.718Z Partial eclipse begins. +1989-02-20T14:55:34.839Z Total eclipse begins. +1989-02-20T15:35:19.956Z Peak of total eclipse. +1989-02-20T16:15:05.072Z Total eclipse ends. +1989-02-20T17:27:15.194Z Partial eclipse ends. -1989-08-17T01:21:19.201Z Partial eclipse begins. -1989-08-17T02:20:31.653Z Total eclipse begins. -1989-08-17T03:08:46.679Z Peak of total eclipse. -1989-08-17T03:57:01.705Z Total eclipse ends. -1989-08-17T04:56:14.157Z Partial eclipse ends. +1989-08-17T01:20:43.868Z Partial eclipse begins. +1989-08-17T02:19:56.976Z Total eclipse begins. +1989-08-17T03:08:10.847Z Peak of total eclipse. +1989-08-17T03:56:24.719Z Total eclipse ends. +1989-08-17T04:55:37.826Z Partial eclipse ends. -1990-02-09T17:29:19.371Z Partial eclipse begins. -1990-02-09T18:50:02.143Z Total eclipse begins. -1990-02-09T19:11:46.639Z Peak of total eclipse. -1990-02-09T19:33:31.135Z Total eclipse ends. -1990-02-09T20:54:13.906Z Partial eclipse ends. +1990-02-09T17:28:36.915Z Partial eclipse begins. +1990-02-09T18:49:12.803Z Total eclipse begins. +1990-02-09T19:11:06.009Z Peak of total eclipse. +1990-02-09T19:32:59.215Z Total eclipse ends. +1990-02-09T20:53:35.103Z Partial eclipse ends. -1990-08-06T12:44:54.290Z Partial eclipse begins. -1990-08-06T14:13:00.287Z Peak of partial eclipse. -1990-08-06T15:41:06.284Z Partial eclipse ends. +1990-08-06T12:44:10.959Z Partial eclipse begins. +1990-08-06T14:12:20.199Z Peak of partial eclipse. +1990-08-06T15:40:29.439Z Partial eclipse ends. -1991-12-21T10:00:27.419Z Partial eclipse begins. -1991-12-21T10:33:39.370Z Peak of partial eclipse. -1991-12-21T11:06:51.320Z Partial eclipse ends. +1991-12-21T10:00:02.660Z Partial eclipse begins. +1991-12-21T10:33:04.084Z Peak of partial eclipse. +1991-12-21T11:06:05.507Z Partial eclipse ends. -1992-06-15T03:27:15.071Z Partial eclipse begins. -1992-06-15T04:57:38.272Z Peak of partial eclipse. -1992-06-15T06:28:01.474Z Partial eclipse ends. +1992-06-15T03:26:36.541Z Partial eclipse begins. +1992-06-15T04:56:56.392Z Peak of partial eclipse. +1992-06-15T06:27:16.242Z Partial eclipse ends. -1992-12-09T21:59:59.280Z Partial eclipse begins. -1992-12-09T23:07:16.402Z Total eclipse begins. -1992-12-09T23:44:42.637Z Peak of total eclipse. -1992-12-10T00:22:08.873Z Total eclipse ends. -1992-12-10T01:29:25.994Z Partial eclipse ends. +1992-12-09T21:59:22.080Z Partial eclipse begins. +1992-12-09T23:06:41.497Z Total eclipse begins. +1992-12-09T23:44:04.198Z Peak of total eclipse. +1992-12-10T00:21:26.899Z Total eclipse ends. +1992-12-10T01:28:46.317Z Partial eclipse ends. -1993-06-04T11:11:48.333Z Partial eclipse begins. -1993-06-04T12:12:49.449Z Total eclipse begins. -1993-06-04T13:01:01.725Z Peak of total eclipse. -1993-06-04T13:49:14.001Z Total eclipse ends. -1993-06-04T14:50:15.117Z Partial eclipse ends. +1993-06-04T11:11:10.310Z Partial eclipse begins. +1993-06-04T12:12:10.636Z Total eclipse begins. +1993-06-04T13:00:24.284Z Peak of total eclipse. +1993-06-04T13:48:37.931Z Total eclipse ends. +1993-06-04T14:49:38.257Z Partial eclipse ends. diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 7c908617..2af45adb 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -221,6 +221,9 @@ class Vector: def __sub__(self, other): return Vector(self.x - other.x, self.y - other.y, self.z - other.z, self.t) + def __neg__(self): + return Vector(-self.x, -self.y, -self.z, self.t) + def format(self, coord_format): """Returns a custom format string representation of the vector.""" layout = '({:' + coord_format + '}, {:' + coord_format + '}, {:' + coord_format + '}, {})' @@ -6013,9 +6016,12 @@ def _CalcShadow(body_radius_km, time, target, sdir): def _EarthShadow(time): - e = _CalcEarth(time) + # This function helps find when the Earth's shadow falls upon the Moon. + # Light-travel and aberration corrected vector from the Earth to the Sun. + # The negative vector -s is thus the path of sunlight through the center of the Earth. + s = GeoVector(Body.Sun, time, True) m = GeoMoon(time) - return _CalcShadow(_EARTH_ECLIPSE_RADIUS_KM, time, m, e) + return _CalcShadow(_EARTH_ECLIPSE_RADIUS_KM, time, m, -s) def _MoonShadow(time): @@ -6135,6 +6141,12 @@ class LunarEclipseInfo: The `kind` field thus holds one of the values `EclipseKind.Penumbral`, `EclipseKind.Partial`, or `EclipseKind.Total`, depending on the kind of lunar eclipse found. + The `obscuration` field holds a value in the range [0, 1] that indicates what fraction + of the Moon's apparent disc area is covered by the Earth's umbra at the eclipse's peak. + This indicates how dark the peak eclipse appears. For penumbral eclipses, the obscuration + is 0, because the Moon does not pass through the Earth's umbra. For partial eclipses, + the obscuration is somewhere between 0 and 1. For total lunar eclipses, the obscuration is 1. + Field `peak` holds the date and time of the peak of the eclipse, when it is at its peak. Fields `sd_penum`, `sd_partial`, and `sd_total` hold the semi-duration of each phase @@ -6147,6 +6159,8 @@ class LunarEclipseInfo: ---------- kind : EclipseKind The type of lunar eclipse found. + obscuration : float + The peak fraction of the Moon's apparent disc that is covered by the Earth's umbra. peak : Time The time of the eclipse at its peak. sd_penum : float @@ -6156,16 +6170,18 @@ class LunarEclipseInfo: sd_total : float The semi-duration of the penumbral phase in minutes, or 0.0 if none. """ - def __init__(self, kind, peak, sd_penum, sd_partial, sd_total): + def __init__(self, kind, obscuration, peak, sd_penum, sd_partial, sd_total): self.kind = kind + self.obscuration = obscuration self.peak = peak self.sd_penum = sd_penum self.sd_partial = sd_partial self.sd_total = sd_total def __repr__(self): - return 'LunarEclipseInfo({}, peak={}, sd_penum={}, sd_partial={}, sd_total={})'.format( + return 'LunarEclipseInfo({}, obscuration={}, peak={}, sd_penum={}, sd_partial={}, sd_total={})'.format( self.kind, + self.obscuration, repr(self.peak), self.sd_penum, self.sd_partial, @@ -6496,6 +6512,47 @@ def _MoonEclipticLatitudeDegrees(time): moon = _CalcMoon(time) return math.degrees(moon.geo_eclip_lat) +def _Obscuration(a, b, c): + # a = radius of first disc + # b = radius of second disc + # c = distance between the centers of the discs + if a <= 0.0: + raise Error('Radius of first disc must be positive.') + + if b <= 0.0: + raise Error('Radius of second disc must be positive.') + + if c < 0.0: + raise Error('Distance between discs is not allowed to be negative.') + + if c >= a + b: + # The discs are too far apart to have any overlapping area + return 0.0 + + if c == 0.0: + # The discs have a common center. Therefore, one disc is inside the other. + return 1.0 if (a <= b) else (b*b)/(a*a) + + x = (a*a - b*b + c*c) / (2.0*c) + radicand = a*a - x*x + if radicand <= 0.0: + # The circumferences do not intersect, or are tangent. + # We already ruled out the case of non-overlapping discs. + # Therefore, one disc is inside the other. + return 1.0 if (a <= b) else (b*b)/(a*a) + + # The discs overlap fractionally in a pair of lens-shaped areas. + y = math.sqrt(radicand) + + # Return the overlapping fractional area. + # There are two lens-shaped areas, one to the left of x, the other to the right of x. + # Each part is calculated by subtracting a triangular area from a sector's area. + lens1 = a*a*math.acos(x/a) - x*y + lens2 = b*b*math.acos((c-x)/b) - (c-x)*y + + # Find the fractional area with respect to the first disc. + return (lens1 + lens2) / (math.pi*a*a) + def SearchLunarEclipse(startTime): """Searches for a lunar eclipse. @@ -6535,6 +6592,7 @@ def SearchLunarEclipse(startTime): if shadow.r < shadow.p + _MOON_MEAN_RADIUS_KM: # This is at least a penumbral eclipse. We will return a result. kind = EclipseKind.Penumbral + obscuration = 0.0 sd_total = 0.0 sd_partial = 0.0 sd_penum = _ShadowSemiDurationMinutes(shadow.time, shadow.p + _MOON_MEAN_RADIUS_KM, 200.0) @@ -6547,9 +6605,12 @@ def SearchLunarEclipse(startTime): if shadow.r + _MOON_MEAN_RADIUS_KM < shadow.k: # This is a total eclipse. kind = EclipseKind.Total + obscuration = 1.0 sd_total = _ShadowSemiDurationMinutes(shadow.time, shadow.k - _MOON_MEAN_RADIUS_KM, sd_partial) + else: + obscuration = _Obscuration(_MOON_MEAN_RADIUS_KM, shadow.k, shadow.r) - return LunarEclipseInfo(kind, shadow.time, sd_penum, sd_partial, sd_total) + return LunarEclipseInfo(kind, obscuration, shadow.time, sd_penum, sd_partial, sd_total) # We didn't find an eclipse on this full moon, so search for the next one. fmtime = fullmoon.AddDays(10) diff --git a/generate/test.py b/generate/test.py index c7d7ecac..4a64678c 100755 --- a/generate/test.py +++ b/generate/test.py @@ -1403,22 +1403,30 @@ def LunarEclipse(): return 1 partial_minutes = float(token[0]) total_minutes = float(token[1]) - valid = False + sd_valid = False + frac_valid = False # Verify that the calculated eclipse semi-durations are consistent with the kind. + # Verify that obscurations also make sense for the kind. if eclipse.kind == astronomy.EclipseKind.Penumbral: - valid = (eclipse.sd_penum > 0.0) and (eclipse.sd_partial == 0.0) and (eclipse.sd_total == 0.0) + sd_valid = (eclipse.sd_penum > 0.0) and (eclipse.sd_partial == 0.0) and (eclipse.sd_total == 0.0) + frac_valid = (eclipse.obscuration == 0.0) elif eclipse.kind == astronomy.EclipseKind.Partial: - valid = (eclipse.sd_penum > 0.0) and (eclipse.sd_partial > 0.0) and (eclipse.sd_total == 0.0) + sd_valid = (eclipse.sd_penum > 0.0) and (eclipse.sd_partial > 0.0) and (eclipse.sd_total == 0.0) + frac_valid = (0.0 < eclipse.obscuration < 1.0) elif eclipse.kind == astronomy.EclipseKind.Total: - valid = (eclipse.sd_penum > 0.0) and (eclipse.sd_partial > 0.0) and (eclipse.sd_total > 0.0) + sd_valid = (eclipse.sd_penum > 0.0) and (eclipse.sd_partial > 0.0) and (eclipse.sd_total > 0.0) + frac_valid = (eclipse.obscuration == 1.0) else: print('PY LunarEclipse({} line {}): invalid eclipse kind {}.'.format(filename, lnum, eclipse.kind)) return 1 - if not valid: + if not sd_valid: print('PY LunarEclipse({} line {}): invalid semidurations.'.format(filename, lnum)) return 1 + if not frac_valid: + print('PY LunarEclipse({} line {}): invalid obscuration {:0.8f} for eclipsekind {}.'.format(filename, lnum, eclipse.obscuration, eclipse.kind)) + # Check eclipse peak time. diff_days = eclipse.peak.ut - peak_time.ut @@ -2876,6 +2884,51 @@ def DatesIssue250(): #----------------------------------------------------------------------------------------------------------- +def LunarFractionCase(year, month, day, obscuration): + time = astronomy.Time.Make(year, month, day, 0, 0, 0.0) + eclipse = astronomy.SearchLunarEclipse(time) + # This should be a partial lunar eclipse. + if eclipse.kind != astronomy.EclipseKind.Partial: + print('PY LunarFractionCase({:04d}-{:02d}-{:02d}) FAIL: expected partial eclipse, but found {}.'.format(year, month, day, eclipse.kind)) + return 1 + + # The partial eclipse should always happen within 24 hours of the given date. + dt = v(eclipse.peak.ut - time.ut) + if dt < 0.0 or dt > 1.0: + print('PY LunarFractionCase({:04d}-{:02d}-{:02d}) FAIL: eclipse occurs {:0.4f} days after predicted date.'.format(year, month, day, dt)) + return 1 + + diff = v(eclipse.obscuration - obscuration) + if abs(diff) > 0.00763: + print('PY LunarFractionCase({:04d}-{:02d}-{:02d}) FAIL: excessive obscuration diff = {:0.8f}, expected = {:0.8f}, actual = {:0.8f}'.format(year, month, day, diff, obscuration, eclipse.obscuration)) + return 1 + + Debug('PY LunarFractionCase({:04d}-{:02d}-{:02d}): obscuration diff = {:11.8f}'.format(year, month, day, diff)) + return 0 + + +def LunarFraction(): + # Verify calculation of the fraction of the Moon's disc covered by the Earth's umbra during a partial eclipse. + # Data for this is more tedious to gather, because Espenak data does not contain it. + # We already verify fraction=0.0 for penumbral eclipses and fraction=1.0 for total eclipses in LunarEclipseTest. + return ( + LunarFractionCase(2010, 6, 26, 0.506) or # https://www.timeanddate.com/eclipse/lunar/2010-june-26 + LunarFractionCase(2012, 6, 4, 0.304) or # https://www.timeanddate.com/eclipse/lunar/2012-june-4 + LunarFractionCase(2013, 4, 25, 0.003) or # https://www.timeanddate.com/eclipse/lunar/2013-april-25 + LunarFractionCase(2017, 8, 7, 0.169) or # https://www.timeanddate.com/eclipse/lunar/2017-august-7 + LunarFractionCase(2019, 7, 16, 0.654) or # https://www.timeanddate.com/eclipse/lunar/2019-july-16 + LunarFractionCase(2021, 11, 19, 0.991) or # https://www.timeanddate.com/eclipse/lunar/2021-november-19 + LunarFractionCase(2023, 10, 28, 0.060) or # https://www.timeanddate.com/eclipse/lunar/2023-october-28 + LunarFractionCase(2024, 9, 18, 0.035) or # https://www.timeanddate.com/eclipse/lunar/2024-september-18 + LunarFractionCase(2026, 8, 28, 0.962) or # https://www.timeanddate.com/eclipse/lunar/2026-august-28 + LunarFractionCase(2028, 1, 12, 0.024) or # https://www.timeanddate.com/eclipse/lunar/2028-january-12 + LunarFractionCase(2028, 7, 6, 0.325) or # https://www.timeanddate.com/eclipse/lunar/2028-july-6 + LunarFractionCase(2030, 6, 15, 0.464) or # https://www.timeanddate.com/eclipse/lunar/2030-june-15 + Pass('LunarFraction') + ) + +#----------------------------------------------------------------------------------------------------------- + UnitTests = { 'aberration': Aberration, 'axis': Axis, @@ -2895,6 +2948,7 @@ UnitTests = { 'lunar_apsis': LunarApsis, 'lunar_eclipse': LunarEclipse, 'lunar_eclipse_78': LunarEclipseIssue78, + 'lunar_fraction': LunarFraction, 'magnitude': Magnitude, 'moon': GeoMoon, 'moon_nodes': MoonNodes, diff --git a/source/python/README.md b/source/python/README.md index eecb0ac8..1ca733cd 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -702,6 +702,11 @@ Partial eclipses occur when part, but not all, of the Moon touches the Earth's u Total eclipses occur when the entire Moon passes into the Earth's umbra. The `kind` field thus holds one of the values `EclipseKind.Penumbral`, `EclipseKind.Partial`, or `EclipseKind.Total`, depending on the kind of lunar eclipse found. +The `obscuration` field holds a value in the range [0, 1] that indicates what fraction +of the Moon's apparent disc area is covered by the Earth's umbra at the eclipse's peak. +This indicates how dark the peak eclipse appears. For penumbral eclipses, the obscuration +is 0, because the Moon does not pass through the Earth's umbra. For partial eclipses, +the obscuration is somewhere between 0 and 1. For total lunar eclipses, the obscuration is 1. Field `peak` holds the date and time of the peak of the eclipse, when it is at its peak. Fields `sd_penum`, `sd_partial`, and `sd_total` hold the semi-duration of each phase of the eclipse, which is half of the amount of time the eclipse spends in each @@ -712,6 +717,7 @@ may determine the date and time of the beginning/end of each eclipse phase. | Type | Attribute | Description | | --- | --- | --- | | [`EclipseKind`](#EclipseKind) | `kind` | The type of lunar eclipse found. | +| `float` | `obscuration` | The peak fraction of the Moon's apparent disc that is covered by the Earth's umbra. | | [`Time`](#Time) | `peak` | The time of the eclipse at its peak. | | `float` | `sd_penum` | The semi-duration of the penumbral phase in minutes. | | `float` | `sd_partial` | The semi-duration of the penumbral phase in minutes, or 0.0 if none. | diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index dd86e5ec..e449156f 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -221,6 +221,9 @@ class Vector: def __sub__(self, other): return Vector(self.x - other.x, self.y - other.y, self.z - other.z, self.t) + def __neg__(self): + return Vector(-self.x, -self.y, -self.z, self.t) + def format(self, coord_format): """Returns a custom format string representation of the vector.""" layout = '({:' + coord_format + '}, {:' + coord_format + '}, {:' + coord_format + '}, {})' @@ -8506,9 +8509,12 @@ def _CalcShadow(body_radius_km, time, target, sdir): def _EarthShadow(time): - e = _CalcEarth(time) + # This function helps find when the Earth's shadow falls upon the Moon. + # Light-travel and aberration corrected vector from the Earth to the Sun. + # The negative vector -s is thus the path of sunlight through the center of the Earth. + s = GeoVector(Body.Sun, time, True) m = GeoMoon(time) - return _CalcShadow(_EARTH_ECLIPSE_RADIUS_KM, time, m, e) + return _CalcShadow(_EARTH_ECLIPSE_RADIUS_KM, time, m, -s) def _MoonShadow(time): @@ -8628,6 +8634,12 @@ class LunarEclipseInfo: The `kind` field thus holds one of the values `EclipseKind.Penumbral`, `EclipseKind.Partial`, or `EclipseKind.Total`, depending on the kind of lunar eclipse found. + The `obscuration` field holds a value in the range [0, 1] that indicates what fraction + of the Moon's apparent disc area is covered by the Earth's umbra at the eclipse's peak. + This indicates how dark the peak eclipse appears. For penumbral eclipses, the obscuration + is 0, because the Moon does not pass through the Earth's umbra. For partial eclipses, + the obscuration is somewhere between 0 and 1. For total lunar eclipses, the obscuration is 1. + Field `peak` holds the date and time of the peak of the eclipse, when it is at its peak. Fields `sd_penum`, `sd_partial`, and `sd_total` hold the semi-duration of each phase @@ -8640,6 +8652,8 @@ class LunarEclipseInfo: ---------- kind : EclipseKind The type of lunar eclipse found. + obscuration : float + The peak fraction of the Moon's apparent disc that is covered by the Earth's umbra. peak : Time The time of the eclipse at its peak. sd_penum : float @@ -8649,16 +8663,18 @@ class LunarEclipseInfo: sd_total : float The semi-duration of the penumbral phase in minutes, or 0.0 if none. """ - def __init__(self, kind, peak, sd_penum, sd_partial, sd_total): + def __init__(self, kind, obscuration, peak, sd_penum, sd_partial, sd_total): self.kind = kind + self.obscuration = obscuration self.peak = peak self.sd_penum = sd_penum self.sd_partial = sd_partial self.sd_total = sd_total def __repr__(self): - return 'LunarEclipseInfo({}, peak={}, sd_penum={}, sd_partial={}, sd_total={})'.format( + return 'LunarEclipseInfo({}, obscuration={}, peak={}, sd_penum={}, sd_partial={}, sd_total={})'.format( self.kind, + self.obscuration, repr(self.peak), self.sd_penum, self.sd_partial, @@ -8989,6 +9005,47 @@ def _MoonEclipticLatitudeDegrees(time): moon = _CalcMoon(time) return math.degrees(moon.geo_eclip_lat) +def _Obscuration(a, b, c): + # a = radius of first disc + # b = radius of second disc + # c = distance between the centers of the discs + if a <= 0.0: + raise Error('Radius of first disc must be positive.') + + if b <= 0.0: + raise Error('Radius of second disc must be positive.') + + if c < 0.0: + raise Error('Distance between discs is not allowed to be negative.') + + if c >= a + b: + # The discs are too far apart to have any overlapping area + return 0.0 + + if c == 0.0: + # The discs have a common center. Therefore, one disc is inside the other. + return 1.0 if (a <= b) else (b*b)/(a*a) + + x = (a*a - b*b + c*c) / (2.0*c) + radicand = a*a - x*x + if radicand <= 0.0: + # The circumferences do not intersect, or are tangent. + # We already ruled out the case of non-overlapping discs. + # Therefore, one disc is inside the other. + return 1.0 if (a <= b) else (b*b)/(a*a) + + # The discs overlap fractionally in a pair of lens-shaped areas. + y = math.sqrt(radicand) + + # Return the overlapping fractional area. + # There are two lens-shaped areas, one to the left of x, the other to the right of x. + # Each part is calculated by subtracting a triangular area from a sector's area. + lens1 = a*a*math.acos(x/a) - x*y + lens2 = b*b*math.acos((c-x)/b) - (c-x)*y + + # Find the fractional area with respect to the first disc. + return (lens1 + lens2) / (math.pi*a*a) + def SearchLunarEclipse(startTime): """Searches for a lunar eclipse. @@ -9028,6 +9085,7 @@ def SearchLunarEclipse(startTime): if shadow.r < shadow.p + _MOON_MEAN_RADIUS_KM: # This is at least a penumbral eclipse. We will return a result. kind = EclipseKind.Penumbral + obscuration = 0.0 sd_total = 0.0 sd_partial = 0.0 sd_penum = _ShadowSemiDurationMinutes(shadow.time, shadow.p + _MOON_MEAN_RADIUS_KM, 200.0) @@ -9040,9 +9098,12 @@ def SearchLunarEclipse(startTime): if shadow.r + _MOON_MEAN_RADIUS_KM < shadow.k: # This is a total eclipse. kind = EclipseKind.Total + obscuration = 1.0 sd_total = _ShadowSemiDurationMinutes(shadow.time, shadow.k - _MOON_MEAN_RADIUS_KM, sd_partial) + else: + obscuration = _Obscuration(_MOON_MEAN_RADIUS_KM, shadow.k, shadow.r) - return LunarEclipseInfo(kind, shadow.time, sd_penum, sd_partial, sd_total) + return LunarEclipseInfo(kind, obscuration, shadow.time, sd_penum, sd_partial, sd_total) # We didn't find an eclipse on this full moon, so search for the next one. fmtime = fullmoon.AddDays(10)