From 9187e3e96633b4e2c385fe19ae05735e6c2e7fff Mon Sep 17 00:00:00 2001 From: Don Cross Date: Thu, 4 Jun 2020 19:01:39 -0400 Subject: [PATCH] Python: Implemented GlobalSolarEclipse. --- generate/template/astronomy.py | 270 +++++++++++++++++++++++++++++++-- generate/test.py | 123 +++++++++++++-- source/python/README.md | 65 ++++++++ source/python/astronomy.py | 270 +++++++++++++++++++++++++++++++-- 4 files changed, 685 insertions(+), 43 deletions(-) diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 72361cd3..969b26cb 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -4128,34 +4128,55 @@ class EclipseKind(enum.Enum): Total = 4 -class _EarthShadowInfo: - def __init__(self, time, u, r, k, p): +class _ShadowInfo: + def __init__(self, time, u, r, k, p, target, dir): self.time = time self.u = u # dot product of (heliocentric earth) and (geocentric moon): defines the shadow plane where the Moon is self.r = r # km distance between center of Moon and the line passing through the centers of the Sun and Earth. self.k = k # umbra radius in km, at the shadow plane self.p = p # penumbra radius in km, at the shadow plane + self.target = target # vector from center of shadow-casting body to target location that might receive the shadow + self.dir = dir # vector from center of Sun to center of shadow-casting body + + +def _CalcShadow(body_radius_km, time, target, dir): + u = (dir.x*target.x + dir.y*target.y + dir.z*target.z) / (dir.x*dir.x + dir.y*dir.y + dir.z*dir.z) + dx = (u * dir.x) - target.x + dy = (u * dir.y) - target.y + dz = (u * dir.z) - target.z + r = _KM_PER_AU * math.sqrt(dx*dx + dy*dy + dz*dz) + k = +_SUN_RADIUS_KM - (1.0 + u)*(_SUN_RADIUS_KM - body_radius_km) + p = -_SUN_RADIUS_KM + (1.0 + u)*(_SUN_RADIUS_KM + body_radius_km) + return _ShadowInfo(time, u, r, k, p, target, dir) def _EarthShadow(time): e = _CalcEarth(time) m = GeoMoon(time) - u = (e.x*m.x + e.y*m.y + e.z*m.z) / (e.x*e.x + e.y*e.y + e.z*e.z) - dx = (u * e.x) - m.x - dy = (u * e.y) - m.y - dz = (u * e.z) - m.z - r = _KM_PER_AU * math.sqrt(dx*dx + dy*dy + dz*dz) - k = +_SUN_RADIUS_KM - (1.0 + u)*(_SUN_RADIUS_KM - _EARTH_ECLIPSE_RADIUS_KM) - p = -_SUN_RADIUS_KM + (1.0 + u)*(_SUN_RADIUS_KM + _EARTH_ECLIPSE_RADIUS_KM) - return _EarthShadowInfo(time, u, r, k, p) + return _CalcShadow(_EARTH_ECLIPSE_RADIUS_KM, time, m, e) -def _EarthShadowSlope(context, time): +def _MoonShadow(time): + # This is a variation on the logic in _EarthShadow(). + # Instead of a heliocentric Earth and a geocentric Moon, + # we want a heliocentric Moon and a lunacentric Earth. + h = _CalcEarth(time) # heliocentric Earth + m = GeoMoon(time) # geocentric Moon + # Calculate lunacentric Earth. + e = Vector(-m.x, -m.y, -m.z, m.t) + # Convert geocentric moon to heliocentric Moon. + m.x += h.x + m.y += h.y + m.z += h.z + return _CalcShadow(_MOON_MEAN_RADIUS_KM, time, e, m) + + +def _ShadowDistanceSlope(shadowfunc, time): dt = 1.0 / 86400.0 t1 = time.AddDays(-dt) t2 = time.AddDays(+dt) - shadow1 = _EarthShadow(t1) - shadow2 = _EarthShadow(t2) + shadow1 = shadowfunc(t1) + shadow2 = shadowfunc(t2) return (shadow2.r - shadow1.r) / dt @@ -4163,9 +4184,17 @@ def _PeakEarthShadow(search_center_time): window = 0.03 # initial search window, in days, before/after given time t1 = search_center_time.AddDays(-window) t2 = search_center_time.AddDays(+window) - tx = Search(_EarthShadowSlope, None, t1, t2, 1.0) + tx = Search(_ShadowDistanceSlope, _EarthShadow, t1, t2, 1.0) return _EarthShadow(tx) + +def _PeakMoonShadow(search_center_time): + window = 0.03 # initial search window, in days, before/after given time + t1 = search_center_time.AddDays(-window) + t2 = search_center_time.AddDays(+window) + tx = Search(_ShadowDistanceSlope, _MoonShadow, t1, t2, 1.0) + return _MoonShadow(tx) + class _ShadowDiffContext: def __init__(self, radius_limit, direction): self.radius_limit = radius_limit @@ -4218,6 +4247,148 @@ class LunarEclipseInfo: self.sd_total = sd_total +class GlobalSolarEclipseInfo: + """Reports the time and geographic location of the peak of a solar eclipse. + + Returned by #SearchGlobalSolarEclipse or #NextGlobalSolarEclipse + to report information about a solar eclipse event. + + Field `peak` holds the date and time of the peak of the eclipse, defined as + the instant when the axis of the Moon's shadow cone passes closest to the Earth's center. + + The eclipse is classified as partial, annular, or total, depending on the + maximum amount of the Sun's disc obscured, as seen at the peak location + on the surface of the Earth. + + The `kind` field thus holds `EclipseKind.Partial`, `EclipseKind.Annular`, or `EclipseKind.Total`. + A total eclipse is when the peak observer sees the Sun completely blocked by the Moon. + An annular eclipse is like a total eclipse, but the Moon is too far from the Earth's surface + to completely block the Sun; instead, the Sun takes on a ring-shaped appearance. + A partial eclipse is when the Moon blocks part of the Sun's disc, but nobody on the Earth + observes either a total or annular eclipse. + + If `kind` is `EclipseKind.Total` or `EclipseKind.Annular`, the `latitude` and `longitude` + fields give the geographic coordinates of the center of the Moon's shadow projected + onto the daytime side of the Earth at the instant of the eclipse's peak. + If `kind` has any other value, `latitude` and `longitude` are undefined and should + not be used. + """ + def __init__(self, kind, peak, distance, latitude, longitude): + self.kind = kind + self.peak = peak + self.distance = distance + self.latitude = latitude + self.longitude = longitude + + +def _EclipseKindFromUmbra(k): + # The umbra radius tells us what kind of eclipse the observer sees. + # If the umbra radius is positive, this is a total eclipse. Otherwise, it's annular. + # HACK: I added a tiny bias (14 meters) to match Espenak test data. + if k > 0.014: + return EclipseKind.Total + return EclipseKind.Annular + + +def _GeoidIntersect(shadow): + #astro_global_solar_eclipse_t eclipse; + #astro_rotation_t rot, inv; + #astro_vector_t v, e, o; + #shadow_t surface; + #double A, B, C, radic, u, R; + #double px, py, pz, proj; + #double gast; + + kind = EclipseKind.Partial + peak = shadow.time + distance = shadow.r + latitude = longitude = math.nan + + # We want to calculate the intersection of the shadow axis with the Earth's geoid. + # First we must convert EQJ (equator of J2000) coordinates to EQD (equator of date) + # coordinates that are perfectly aligned with the Earth's equator at this + # moment in time. + rot = Rotation_EQJ_EQD(shadow.time) + v = RotateVector(rot, shadow.dir) # shadow-axis vector in equator-of-date coordinates + e = RotateVector(rot, shadow.target) # lunacentric Earth in equator-of-date coordinates + + # Convert all distances from AU to km. + # But dilate the z-coordinates so that the Earth becomes a perfect sphere. + # Then find the intersection of the vector with the sphere. + # See p 184 in Montenbruck & Pfleger's "Astronomy on the Personal Computer", second edition. + v.x *= _KM_PER_AU + v.y *= _KM_PER_AU + v.z *= _KM_PER_AU / _EARTH_FLATTENING + e.x *= _KM_PER_AU + e.y *= _KM_PER_AU + e.z *= _KM_PER_AU / _EARTH_FLATTENING + + # Solve the quadratic equation that finds whether and where + # the shadow axis intersects with the Earth in the dilated coordinate system. + R = _EARTH_EQUATORIAL_RADIUS_KM + A = v.x*v.x + v.y*v.y + v.z*v.z + B = -2.0 * (v.x*e.x + v.y*e.y + v.z*e.z) + C = (e.x*e.x + e.y*e.y + e.z*e.z) - R*R + radic = B*B - 4*A*C + + if radic > 0.0: + # Calculate the closer of the two intersection points. + # This will be on the day side of the Earth. + u = (-B - math.sqrt(radic)) / (2 * A) + + # Convert lunacentric dilated coordinates to geocentric coordinates. + px = u*v.x - e.x + py = u*v.y - e.y + pz = (u*v.z - e.z) * _EARTH_FLATTENING + + # Convert cartesian coordinates into geodetic latitude/longitude. + proj = math.sqrt(px*px + py*py) * (_EARTH_FLATTENING * _EARTH_FLATTENING) + if proj == 0.0: + if pz > 0.0: + latitude = +90.0 + else: + latitude = -90.0 + else: + latitude = math.degrees(math.atan(pz / proj)) + + # Adjust longitude for Earth's rotation at the given UT. + gast = _sidereal_time(peak) + longitude = math.fmod(math.degrees(math.atan2(py, px)) - (15*gast), 360.0) + if longitude <= -180.0: + longitude += 360.0 + elif longitude > +180.0: + longitude -= 360.0 + + # We want to determine whether the observer sees a total eclipse or an annular eclipse. + # We need to perform a series of vector calculations... + # Calculate the inverse rotation matrix, so we can convert EQD to EQJ. + inv = InverseRotation(rot) + + # Put the EQD geocentric coordinates of the observer into the vector 'o'. + # Also convert back from kilometers to astronomical units. + o = Vector(px / _KM_PER_AU, py / _KM_PER_AU, pz / _KM_PER_AU, shadow.time) + + # Rotate the observer's geocentric EQD back to the EQJ system. + o = RotateVector(inv, o) + + # Convert geocentric vector to lunacentric vector. + o.x += shadow.target.x + o.y += shadow.target.y + o.z += shadow.target.z + + # Recalculate the shadow using a vector from the Moon's center toward the observer. + surface = _CalcShadow(_MOON_POLAR_RADIUS_KM, shadow.time, o, shadow.dir) + + # If we did everything right, the shadow distance should be very close to zero. + # That's because we already determined the observer 'o' is on the shadow axis! + if surface.r > 1.0e-9 or surface.r < 0.0: + raise Error('Unexpected shadow distance from geoid intersection = {}'.format(surface.r)) + + kind = _EclipseKindFromUmbra(surface.k) + + return GlobalSolarEclipseInfo(kind, peak, distance, latitude, longitude) + + def _ShadowSemiDurationMinutes(center_time, radius_limit, window_minutes): # Search backwards and forwards from the center time until shadow axis distance crosses radius limit. window = window_minutes / (24.0 * 60.0) @@ -4310,4 +4481,73 @@ def NextLunarEclipse(prevEclipseTime): LunarEclipseInfo """ startTime = prevEclipseTime.AddDays(10.0) - return SearchLunarEclipse(startTime) \ No newline at end of file + return SearchLunarEclipse(startTime) + + +def SearchGlobalSolarEclipse(startTime): + """Searches for a solar eclipse visible anywhere on the Earth's surface. + + This function finds the first solar eclipse that occurs after `startTime`. + A solar eclipse found may be partial, annular, or total. + See #GlobalSolarEclipseInfo for more information. + To find a series of solar eclipses, call this function once, + then keep calling #NextGlobalSolarEclipse as many times as desired, + passing in the `peak` value returned from the previous call. + + Parameters + ---------- + startTime : Time + The date and time for starting the search for a solar eclipse. + + Returns + ------- + GlobalSolarEclipseInfo + """ + PruneLatitude = 1.8 # Moon's ecliptic latitude beyond which eclipse is impossible + # Iterate through consecutive new moons until we find a solar eclipse visible somewhere on Earth. + nmtime = startTime + for nmcount in range(12): + # Search for the next new moon. Any eclipse will be near it. + newmoon = SearchMoonPhase(0.0, nmtime, 40.0) + if newmoon is None: + raise Error('Cannot find new moon') + + # Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible. + mp = _CalcMoon(newmoon) + if math.degrees(abs(mp.geo_eclip_lat)) < PruneLatitude: + # Search near the new moon for the time when the center of the Earth + # is closest to the line passing through the centers of the Sun and Moon. + shadow = _PeakMoonShadow(newmoon) + if shadow.r < shadow.p + _EARTH_MEAN_RADIUS_KM: + # This is at least a partial solar eclipse visible somewhere on Earth. + # Try to find an intersection between the shadow axis and the Earth's oblate geoid. + return _GeoidIntersect(shadow) + + # We didn't find an eclipse on this new moon, so search for the next one. + nmtime = newmoon.AddDays(10.0) + + # Safety valve to prevent infinite loop. + # This should never happen, because at least 2 solar eclipses happen per year. + raise Error('Failed to find solar eclipse within 12 full moons.') + + +def NextGlobalSolarEclipse(prevEclipseTime): + """Searches for the next global solar eclipse in a series. + + After using #SearchGlobalSolarEclipse to find the first solar eclipse + in a series, you can call this function to find the next consecutive solar eclipse. + Pass in the `peak` value from the #GlobalSolarEclipseInfo returned by the + previous call to `SearchGlobalSolarEclipse` or `NextGlobalSolarEclipse` + to find the next solar eclipse. + + Parameters + ---------- + prevEclipseTime : Time + A date and time near a new moon. Solar eclipse search will start at the next new moon. + + Returns + ------- + GlobalSolarEclipseInfo + """ + startTime = prevEclipseTime.AddDays(10.0) + return SearchGlobalSolarEclipse(startTime) diff --git a/generate/test.py b/generate/test.py index 41b8b0eb..7864bc89 100755 --- a/generate/test.py +++ b/generate/test.py @@ -1221,20 +1221,117 @@ def LunarEclipse(): #----------------------------------------------------------------------------------------------------------- +def VectorFromAngles(lat, lon): + coslat = math.cos(math.radians(lat)) + return [ + math.cos(math.radians(lon)) * coslat, + math.sin(math.radians(lon)) * coslat, + math.sin(math.radians(lat)) + ] + + +def AngleDiff(alat, alon, blat, blon): + a = VectorFromAngles(alat, alon) + b = VectorFromAngles(blat, blon) + dot = a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + if dot <= -1.0: + return 180.0 + if dot >= +1.0: + return 0.0 + return math.degrees(math.acos(dot)) + + +def GlobalSolarEclipse(): + expected_count = 1180 + max_minutes = 0.0 + max_angle = 0.0 + skip_count = 0 + eclipse = astronomy.SearchGlobalSolarEclipse(astronomy.Time.Make(1701, 1, 1, 0, 0, 0)) + filename = 'eclipse/solar_eclipse.txt' + with open(filename, 'rt') as infile: + lnum = 0 + for line in infile: + lnum += 1 + # 1889-12-22T12:54:15Z -6 T -12.7 -12.8 + token = line.split() + if len(token) != 5: + print('GlobalSolarEclipse({} line {}): invalid token count = {}'.format(filename, lnum, len(token))) + return 1 + peak = astronomy.Time.Parse(token[0]) + typeChar = token[2] + lat = float(token[3]) + lon = float(token[4]) + expected_kind = { + 'P': astronomy.EclipseKind.Partial, + 'A': astronomy.EclipseKind.Annular, + 'T': astronomy.EclipseKind.Total, + 'H': astronomy.EclipseKind.Total, + }[typeChar] + + diff_days = eclipse.peak.tt - peak.tt + # Sometimes we find marginal eclipses that aren't listed in the test data. + # Ignore them if the distance between the Sun/Moon shadow axis and the Earth's center is large. + while diff_days < -25.0 and eclipse.distance > 9000.0: + skip_count += 1 + eclipse = astronomy.NextGlobalSolarEclipse(eclipse.peak) + diff_days = eclipse.peak.ut - peak.ut + + # Validate the eclipse prediction. + diff_minutes = (24 * 60) * abs(diff_days) + if diff_minutes > 6.93: + print('PY GlobalSolarEclipseTest({} line {}): EXCESSIVE TIME ERROR = {} minutes'.format(filename, lnum, diff_minutes)) + return 1 + + if diff_minutes > max_minutes: + max_minutes = diff_minutes + + # Validate the eclipse kind, but only when it is not a "glancing" eclipse. + if (eclipse.distance < 6360) and (eclipse.kind != expected_kind): + print('PY GlobalSolarEclipseTest({} line {}): WRONG ECLIPSE KIND: expected {}, found {}'.format(filename, lnum, expected_kind, eclipse.kind)) + return 1 + + if eclipse.kind == astronomy.EclipseKind.Total or eclipse.kind == astronomy.EclipseKind.Annular: + # When the distance between the Moon's shadow ray and the Earth's center is beyond 6100 km, + # it creates a glancing blow whose geographic coordinates are excessively sensitive to + # slight changes in the ray. Therefore, it is unreasonable to count large errors there. + if eclipse.distance < 6100.0: + diff_angle = AngleDiff(lat, lon, eclipse.latitude, eclipse.longitude) + if diff_angle > 0.247: + print('PY GlobalSolarEclipseTest({} line {}): EXCESSIVE GEOGRAPHIC LOCATION ERROR = {} degrees'.format(filename, lnum, diff_angle)) + return 1 + if diff_angle > max_angle: + max_angle = diff_angle + + eclipse = astronomy.NextGlobalSolarEclipse(eclipse.peak) + + if lnum != expected_count: + print('PY GlobalSolarEclipseTest: WRONG LINE COUNT = {}, expected {}'.format(lnum, expected_count)) + return 1 + + if skip_count > 2: + print('PY GlobalSolarEclipseTest: EXCESSIVE SKIP COUNT = {}'.format(skip_count)) + return 1 + + print('PY GlobalSolarEclipseTest: PASS ({} verified, {} skipped, max minutes = {}, max angle = {})'.format(lnum, skip_count, max_minutes, max_angle)) + return 0 + +#----------------------------------------------------------------------------------------------------------- + UnitTests = { - 'constellation': Constellation, - 'elongation': Elongation, - 'lunar_apsis': LunarApsis, - 'lunar_eclipse': LunarEclipse, - 'magnitude': Magnitude, - 'moon': GeoMoon, - 'moonphase': MoonPhase, - 'planet_apsis': PlanetApsis, - 'refraction': Refraction, - 'riseset': RiseSet, - 'rotation': Rotation, - 'seasons': Seasons, - 'time': AstroTime, + 'constellation': Constellation, + 'elongation': Elongation, + 'global_solar_eclipse': GlobalSolarEclipse, + 'lunar_apsis': LunarApsis, + 'lunar_eclipse': LunarEclipse, + 'magnitude': Magnitude, + 'moon': GeoMoon, + 'moonphase': MoonPhase, + 'planet_apsis': PlanetApsis, + 'refraction': Refraction, + 'riseset': RiseSet, + 'rotation': Rotation, + 'seasons': Seasons, + 'time': AstroTime, } #----------------------------------------------------------------------------------------------------------- diff --git a/source/python/README.md b/source/python/README.md index f194938c..972b0317 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -240,6 +240,32 @@ equator projected onto the sky. --- + +### class GlobalSolarEclipseInfo + +**Reports the time and geographic location of the peak of a solar eclipse.** + +Returned by [`SearchGlobalSolarEclipse`](#SearchGlobalSolarEclipse) or [`NextGlobalSolarEclipse`](#NextGlobalSolarEclipse) +to report information about a solar eclipse event. +Field `peak` holds the date and time of the peak of the eclipse, defined as +the instant when the axis of the Moon's shadow cone passes closest to the Earth's center. +The eclipse is classified as partial, annular, or total, depending on the +maximum amount of the Sun's disc obscured, as seen at the peak location +on the surface of the Earth. +The `kind` field thus holds `EclipseKind.Partial`, `EclipseKind.Annular`, or `EclipseKind.Total`. +A total eclipse is when the peak observer sees the Sun completely blocked by the Moon. +An annular eclipse is like a total eclipse, but the Moon is too far from the Earth's surface +to completely block the Sun; instead, the Sun takes on a ring-shaped appearance. +A partial eclipse is when the Moon blocks part of the Sun's disc, but nobody on the Earth +observes either a total or annular eclipse. +If `kind` is `EclipseKind.Total` or `EclipseKind.Annular`, the `latitude` and `longitude` +fields give the geographic coordinates of the center of the Moon's shadow projected +onto the daytime side of the Earth at the instant of the eclipse's peak. +If `kind` has any other value, `latitude` and `longitude` are undefined and should +not be used. + +--- + ### class HorizontalCoordinates @@ -1212,6 +1238,25 @@ Certain values of the angle have conventional definitions: --- + +### NextGlobalSolarEclipse(prevEclipseTime) + +**Searches for the next global solar eclipse in a series.** + +After using [`SearchGlobalSolarEclipse`](#SearchGlobalSolarEclipse) to find the first solar eclipse +in a series, you can call this function to find the next consecutive solar eclipse. +Pass in the `peak` value from the [`GlobalSolarEclipseInfo`](#GlobalSolarEclipseInfo) returned by the +previous call to `SearchGlobalSolarEclipse` or `NextGlobalSolarEclipse` +to find the next solar eclipse. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Time`](#Time) | `prevEclipseTime` | A date and time near a new moon. Solar eclipse search will start at the next new moon. | + +### Returns: [`GlobalSolarEclipseInfo`](#GlobalSolarEclipseInfo) + +--- + ### NextLunarApsis(apsis) @@ -1636,6 +1681,26 @@ the function returns `None`. --- + +### SearchGlobalSolarEclipse(startTime) + +**Searches for a solar eclipse visible anywhere on the Earth's surface.** + +This function finds the first solar eclipse that occurs after `startTime`. +A solar eclipse found may be partial, annular, or total. +See [`GlobalSolarEclipseInfo`](#GlobalSolarEclipseInfo) for more information. +To find a series of solar eclipses, call this function once, +then keep calling [`NextGlobalSolarEclipse`](#NextGlobalSolarEclipse) as many times as desired, +passing in the `peak` value returned from the previous call. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Time`](#Time) | `startTime` | The date and time for starting the search for a solar eclipse. | + +### Returns: [`GlobalSolarEclipseInfo`](#GlobalSolarEclipseInfo) + +--- + ### SearchHourAngle(body, observer, hourAngle, startTime) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 75e18714..472ee09d 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -6598,34 +6598,55 @@ class EclipseKind(enum.Enum): Total = 4 -class _EarthShadowInfo: - def __init__(self, time, u, r, k, p): +class _ShadowInfo: + def __init__(self, time, u, r, k, p, target, dir): self.time = time self.u = u # dot product of (heliocentric earth) and (geocentric moon): defines the shadow plane where the Moon is self.r = r # km distance between center of Moon and the line passing through the centers of the Sun and Earth. self.k = k # umbra radius in km, at the shadow plane self.p = p # penumbra radius in km, at the shadow plane + self.target = target # vector from center of shadow-casting body to target location that might receive the shadow + self.dir = dir # vector from center of Sun to center of shadow-casting body + + +def _CalcShadow(body_radius_km, time, target, dir): + u = (dir.x*target.x + dir.y*target.y + dir.z*target.z) / (dir.x*dir.x + dir.y*dir.y + dir.z*dir.z) + dx = (u * dir.x) - target.x + dy = (u * dir.y) - target.y + dz = (u * dir.z) - target.z + r = _KM_PER_AU * math.sqrt(dx*dx + dy*dy + dz*dz) + k = +_SUN_RADIUS_KM - (1.0 + u)*(_SUN_RADIUS_KM - body_radius_km) + p = -_SUN_RADIUS_KM + (1.0 + u)*(_SUN_RADIUS_KM + body_radius_km) + return _ShadowInfo(time, u, r, k, p, target, dir) def _EarthShadow(time): e = _CalcEarth(time) m = GeoMoon(time) - u = (e.x*m.x + e.y*m.y + e.z*m.z) / (e.x*e.x + e.y*e.y + e.z*e.z) - dx = (u * e.x) - m.x - dy = (u * e.y) - m.y - dz = (u * e.z) - m.z - r = _KM_PER_AU * math.sqrt(dx*dx + dy*dy + dz*dz) - k = +_SUN_RADIUS_KM - (1.0 + u)*(_SUN_RADIUS_KM - _EARTH_ECLIPSE_RADIUS_KM) - p = -_SUN_RADIUS_KM + (1.0 + u)*(_SUN_RADIUS_KM + _EARTH_ECLIPSE_RADIUS_KM) - return _EarthShadowInfo(time, u, r, k, p) + return _CalcShadow(_EARTH_ECLIPSE_RADIUS_KM, time, m, e) -def _EarthShadowSlope(context, time): +def _MoonShadow(time): + # This is a variation on the logic in _EarthShadow(). + # Instead of a heliocentric Earth and a geocentric Moon, + # we want a heliocentric Moon and a lunacentric Earth. + h = _CalcEarth(time) # heliocentric Earth + m = GeoMoon(time) # geocentric Moon + # Calculate lunacentric Earth. + e = Vector(-m.x, -m.y, -m.z, m.t) + # Convert geocentric moon to heliocentric Moon. + m.x += h.x + m.y += h.y + m.z += h.z + return _CalcShadow(_MOON_MEAN_RADIUS_KM, time, e, m) + + +def _ShadowDistanceSlope(shadowfunc, time): dt = 1.0 / 86400.0 t1 = time.AddDays(-dt) t2 = time.AddDays(+dt) - shadow1 = _EarthShadow(t1) - shadow2 = _EarthShadow(t2) + shadow1 = shadowfunc(t1) + shadow2 = shadowfunc(t2) return (shadow2.r - shadow1.r) / dt @@ -6633,9 +6654,17 @@ def _PeakEarthShadow(search_center_time): window = 0.03 # initial search window, in days, before/after given time t1 = search_center_time.AddDays(-window) t2 = search_center_time.AddDays(+window) - tx = Search(_EarthShadowSlope, None, t1, t2, 1.0) + tx = Search(_ShadowDistanceSlope, _EarthShadow, t1, t2, 1.0) return _EarthShadow(tx) + +def _PeakMoonShadow(search_center_time): + window = 0.03 # initial search window, in days, before/after given time + t1 = search_center_time.AddDays(-window) + t2 = search_center_time.AddDays(+window) + tx = Search(_ShadowDistanceSlope, _MoonShadow, t1, t2, 1.0) + return _MoonShadow(tx) + class _ShadowDiffContext: def __init__(self, radius_limit, direction): self.radius_limit = radius_limit @@ -6688,6 +6717,148 @@ class LunarEclipseInfo: self.sd_total = sd_total +class GlobalSolarEclipseInfo: + """Reports the time and geographic location of the peak of a solar eclipse. + + Returned by #SearchGlobalSolarEclipse or #NextGlobalSolarEclipse + to report information about a solar eclipse event. + + Field `peak` holds the date and time of the peak of the eclipse, defined as + the instant when the axis of the Moon's shadow cone passes closest to the Earth's center. + + The eclipse is classified as partial, annular, or total, depending on the + maximum amount of the Sun's disc obscured, as seen at the peak location + on the surface of the Earth. + + The `kind` field thus holds `EclipseKind.Partial`, `EclipseKind.Annular`, or `EclipseKind.Total`. + A total eclipse is when the peak observer sees the Sun completely blocked by the Moon. + An annular eclipse is like a total eclipse, but the Moon is too far from the Earth's surface + to completely block the Sun; instead, the Sun takes on a ring-shaped appearance. + A partial eclipse is when the Moon blocks part of the Sun's disc, but nobody on the Earth + observes either a total or annular eclipse. + + If `kind` is `EclipseKind.Total` or `EclipseKind.Annular`, the `latitude` and `longitude` + fields give the geographic coordinates of the center of the Moon's shadow projected + onto the daytime side of the Earth at the instant of the eclipse's peak. + If `kind` has any other value, `latitude` and `longitude` are undefined and should + not be used. + """ + def __init__(self, kind, peak, distance, latitude, longitude): + self.kind = kind + self.peak = peak + self.distance = distance + self.latitude = latitude + self.longitude = longitude + + +def _EclipseKindFromUmbra(k): + # The umbra radius tells us what kind of eclipse the observer sees. + # If the umbra radius is positive, this is a total eclipse. Otherwise, it's annular. + # HACK: I added a tiny bias (14 meters) to match Espenak test data. + if k > 0.014: + return EclipseKind.Total + return EclipseKind.Annular + + +def _GeoidIntersect(shadow): + #astro_global_solar_eclipse_t eclipse; + #astro_rotation_t rot, inv; + #astro_vector_t v, e, o; + #shadow_t surface; + #double A, B, C, radic, u, R; + #double px, py, pz, proj; + #double gast; + + kind = EclipseKind.Partial + peak = shadow.time + distance = shadow.r + latitude = longitude = math.nan + + # We want to calculate the intersection of the shadow axis with the Earth's geoid. + # First we must convert EQJ (equator of J2000) coordinates to EQD (equator of date) + # coordinates that are perfectly aligned with the Earth's equator at this + # moment in time. + rot = Rotation_EQJ_EQD(shadow.time) + v = RotateVector(rot, shadow.dir) # shadow-axis vector in equator-of-date coordinates + e = RotateVector(rot, shadow.target) # lunacentric Earth in equator-of-date coordinates + + # Convert all distances from AU to km. + # But dilate the z-coordinates so that the Earth becomes a perfect sphere. + # Then find the intersection of the vector with the sphere. + # See p 184 in Montenbruck & Pfleger's "Astronomy on the Personal Computer", second edition. + v.x *= _KM_PER_AU + v.y *= _KM_PER_AU + v.z *= _KM_PER_AU / _EARTH_FLATTENING + e.x *= _KM_PER_AU + e.y *= _KM_PER_AU + e.z *= _KM_PER_AU / _EARTH_FLATTENING + + # Solve the quadratic equation that finds whether and where + # the shadow axis intersects with the Earth in the dilated coordinate system. + R = _EARTH_EQUATORIAL_RADIUS_KM + A = v.x*v.x + v.y*v.y + v.z*v.z + B = -2.0 * (v.x*e.x + v.y*e.y + v.z*e.z) + C = (e.x*e.x + e.y*e.y + e.z*e.z) - R*R + radic = B*B - 4*A*C + + if radic > 0.0: + # Calculate the closer of the two intersection points. + # This will be on the day side of the Earth. + u = (-B - math.sqrt(radic)) / (2 * A) + + # Convert lunacentric dilated coordinates to geocentric coordinates. + px = u*v.x - e.x + py = u*v.y - e.y + pz = (u*v.z - e.z) * _EARTH_FLATTENING + + # Convert cartesian coordinates into geodetic latitude/longitude. + proj = math.sqrt(px*px + py*py) * (_EARTH_FLATTENING * _EARTH_FLATTENING) + if proj == 0.0: + if pz > 0.0: + latitude = +90.0 + else: + latitude = -90.0 + else: + latitude = math.degrees(math.atan(pz / proj)) + + # Adjust longitude for Earth's rotation at the given UT. + gast = _sidereal_time(peak) + longitude = math.fmod(math.degrees(math.atan2(py, px)) - (15*gast), 360.0) + if longitude <= -180.0: + longitude += 360.0 + elif longitude > +180.0: + longitude -= 360.0 + + # We want to determine whether the observer sees a total eclipse or an annular eclipse. + # We need to perform a series of vector calculations... + # Calculate the inverse rotation matrix, so we can convert EQD to EQJ. + inv = InverseRotation(rot) + + # Put the EQD geocentric coordinates of the observer into the vector 'o'. + # Also convert back from kilometers to astronomical units. + o = Vector(px / _KM_PER_AU, py / _KM_PER_AU, pz / _KM_PER_AU, shadow.time) + + # Rotate the observer's geocentric EQD back to the EQJ system. + o = RotateVector(inv, o) + + # Convert geocentric vector to lunacentric vector. + o.x += shadow.target.x + o.y += shadow.target.y + o.z += shadow.target.z + + # Recalculate the shadow using a vector from the Moon's center toward the observer. + surface = _CalcShadow(_MOON_POLAR_RADIUS_KM, shadow.time, o, shadow.dir) + + # If we did everything right, the shadow distance should be very close to zero. + # That's because we already determined the observer 'o' is on the shadow axis! + if surface.r > 1.0e-9 or surface.r < 0.0: + raise Error('Unexpected shadow distance from geoid intersection = {}'.format(surface.r)) + + kind = _EclipseKindFromUmbra(surface.k) + + return GlobalSolarEclipseInfo(kind, peak, distance, latitude, longitude) + + def _ShadowSemiDurationMinutes(center_time, radius_limit, window_minutes): # Search backwards and forwards from the center time until shadow axis distance crosses radius limit. window = window_minutes / (24.0 * 60.0) @@ -6780,4 +6951,73 @@ def NextLunarEclipse(prevEclipseTime): LunarEclipseInfo """ startTime = prevEclipseTime.AddDays(10.0) - return SearchLunarEclipse(startTime) \ No newline at end of file + return SearchLunarEclipse(startTime) + + +def SearchGlobalSolarEclipse(startTime): + """Searches for a solar eclipse visible anywhere on the Earth's surface. + + This function finds the first solar eclipse that occurs after `startTime`. + A solar eclipse found may be partial, annular, or total. + See #GlobalSolarEclipseInfo for more information. + To find a series of solar eclipses, call this function once, + then keep calling #NextGlobalSolarEclipse as many times as desired, + passing in the `peak` value returned from the previous call. + + Parameters + ---------- + startTime : Time + The date and time for starting the search for a solar eclipse. + + Returns + ------- + GlobalSolarEclipseInfo + """ + PruneLatitude = 1.8 # Moon's ecliptic latitude beyond which eclipse is impossible + # Iterate through consecutive new moons until we find a solar eclipse visible somewhere on Earth. + nmtime = startTime + for nmcount in range(12): + # Search for the next new moon. Any eclipse will be near it. + newmoon = SearchMoonPhase(0.0, nmtime, 40.0) + if newmoon is None: + raise Error('Cannot find new moon') + + # Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible. + mp = _CalcMoon(newmoon) + if math.degrees(abs(mp.geo_eclip_lat)) < PruneLatitude: + # Search near the new moon for the time when the center of the Earth + # is closest to the line passing through the centers of the Sun and Moon. + shadow = _PeakMoonShadow(newmoon) + if shadow.r < shadow.p + _EARTH_MEAN_RADIUS_KM: + # This is at least a partial solar eclipse visible somewhere on Earth. + # Try to find an intersection between the shadow axis and the Earth's oblate geoid. + return _GeoidIntersect(shadow) + + # We didn't find an eclipse on this new moon, so search for the next one. + nmtime = newmoon.AddDays(10.0) + + # Safety valve to prevent infinite loop. + # This should never happen, because at least 2 solar eclipses happen per year. + raise Error('Failed to find solar eclipse within 12 full moons.') + + +def NextGlobalSolarEclipse(prevEclipseTime): + """Searches for the next global solar eclipse in a series. + + After using #SearchGlobalSolarEclipse to find the first solar eclipse + in a series, you can call this function to find the next consecutive solar eclipse. + Pass in the `peak` value from the #GlobalSolarEclipseInfo returned by the + previous call to `SearchGlobalSolarEclipse` or `NextGlobalSolarEclipse` + to find the next solar eclipse. + + Parameters + ---------- + prevEclipseTime : Time + A date and time near a new moon. Solar eclipse search will start at the next new moon. + + Returns + ------- + GlobalSolarEclipseInfo + """ + startTime = prevEclipseTime.AddDays(10.0) + return SearchGlobalSolarEclipse(startTime)