From d7e86fae52d48c09ae9d61b334bf682486d98e17 Mon Sep 17 00:00:00 2001
From: Don Cross <cosinekitty@gmail.com>
Date: Tue, 5 Apr 2022 12:30:13 -0400
Subject: [PATCH] C#: Search cleanup, add InternalError exception.

Added an InternalError class to explicitly indicate
that an exception occurs due to an internal assertion
failure inside Astronomy Engine. Any InternalError
should be considered a bug in Astronomy Engine, not
a bug in calling code.

Upon reviewing the code for searching moon phases,
I discovered that there was inconsistent behavior
in SearchMoonPhase. It was sometimes returning null,
other times throwing an exception. Because the caller
passes in `limitDays`, it makes sense to simply
return `null` in any case where the search fails.
This is to support callers that intentionally want
to find whether or not a moon phase occurs in a given
small window of time.

Updated internal callers of SearchMoonPhase to throw
an InternalError when they know they should always
find an event.

Internal function FindSeasonChange did not check to
make sure SearchSunLongitude succeeded. There is no
known case where this failure happens, but if it did,
a null AstroTime would have been stored in SeasonsInfo.
It is better to fail early with an explicit InternalError.

Other miscellaneous C# code cleanup.

In the Python code, I found a couple of `raise Error`
that needed to be changed to `raise InternalError`.
---
 demo/python/astronomy.py             |  10 ++-
 generate/template/astronomy.cs       | 130 ++++++++++++++-------------
 generate/template/astronomy.py       |  10 ++-
 source/csharp/README.md              |  11 ++-
 source/csharp/astronomy.cs           | 130 ++++++++++++++-------------
 source/python/astronomy/astronomy.py |  10 ++-
 6 files changed, 166 insertions(+), 135 deletions(-)

diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py
index 0a01fe76..616504e7 100644
--- a/demo/python/astronomy.py
+++ b/demo/python/astronomy.py
@@ -378,7 +378,7 @@ class InternalError(Error):
     Astronomy Engine for everyone! (Thank you in advance from the author.)
     """
     def __init__(self):
-        Error.__init__(self, 'Internal error - please report issue at https://github.com/cosinekitty/astronomy/issues')
+        Error.__init__(self, 'Internal error - please report issue, including stack trace, at https://github.com/cosinekitty/astronomy/issues')
 
 class NoConvergeError(Error):
     """A numeric solver did not converge.
@@ -8667,7 +8667,7 @@ def SearchLunarEclipse(startTime):
         # Search for the next full moon. Any eclipse will be near it.
         fullmoon = SearchMoonPhase(180, fmtime, 40)
         if fullmoon is None:
-            raise Error('Cannot find full moon.')
+            raise InternalError()   # should have always found the next full moon
 
         # Pruning: if the full Moon's ecliptic latitude is too large,
         # a lunar eclipse is not possible. Avoid needless work searching for
@@ -8751,7 +8751,7 @@ def SearchGlobalSolarEclipse(startTime):
         # 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')
+            raise InternalError()   # should always find the next new moon
 
         # Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible.
         eclip_lat = _MoonEclipticLatitudeDegrees(newmoon)
@@ -8826,6 +8826,8 @@ def SearchLocalSolarEclipse(startTime, observer):
     while True:
         # Search for the next new moon. Any eclipse will be near it.
         newmoon = SearchMoonPhase(0.0, nmtime, 40.0)
+        if newmoon is None:
+            raise InternalError()   # should always find the next new moon
 
         # Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible.
         eclip_lat = _MoonEclipticLatitudeDegrees(newmoon)
@@ -9080,7 +9082,7 @@ def SearchMoonNode(startTime):
             kind = NodeEventKind.Ascending if (eclip2.lat > eclip1.lat) else NodeEventKind.Descending
             result = Search(_MoonNodeSearchFunc, kind.value, time1, time2, 1.0)
             if result is None:
-                raise InternalError()
+                raise InternalError()   # should always find the next lunar node
             return NodeEventInfo(kind, result)
         time1 = time2
         eclip1 = eclip2
diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs
index cfdcb722..2b828724 100644
--- a/generate/template/astronomy.cs
+++ b/generate/template/astronomy.cs
@@ -51,7 +51,20 @@ namespace CosineKitty
     {
         /// <summary>Creates an exception indicating that the given body is not valid for this operation.</summary>
         public InvalidBodyException(Body body):
-            base(string.Format("Invalid body: {0}", body))
+            base("Invalid body: " + body)
+            {}
+    }
+
+    /// <summary>
+    /// This exception indicates an unexpected error occurred inside Astronomy Engine.
+    /// Please report any such errors by creating an issue at:
+    /// https://github.com/cosinekitty/astronomy/issues
+    /// </summary>
+    public class InternalError: Exception
+    {
+        /// <summary>Creates an exception indicating that an unexpected error ocurred.</summary>
+        public InternalError(string message):
+            base("Internal error. Please report an issue at: https://github.com/cosinekitty/astronomy/issues. Diagnostic: " + message)
             {}
     }
 
@@ -4001,7 +4014,7 @@ $ASTRO_IAU_DATA()
 
                     ltime = ltime2;
                 }
-                throw new Exception("Light travel time correction did not converge");
+                throw new InternalError("Light travel time correction did not converge");
             }
         }
 
@@ -4700,7 +4713,8 @@ $ASTRO_IAU_DATA()
         private static AstroTime FindSeasonChange(double targetLon, int year, int month, int day)
         {
             var startTime = new AstroTime(year, month, day, 0, 0, 0);
-            return SearchSunLongitude(targetLon, startTime, 4.0);
+            return SearchSunLongitude(targetLon, startTime, 4.0) ??
+                throw new InternalError($"Cannot find solution for Sun longitude {targetLon}");
         }
 
         /// <summary>
@@ -4823,7 +4837,7 @@ $ASTRO_IAU_DATA()
             for(;;)
             {
                 if (++iter > iter_limit)
-                    throw new Exception(string.Format("Search did not converge within {0} iterations.", iter_limit));
+                    throw new InternalError(string.Format("Search did not converge within {0} iterations.", iter_limit));
 
                 double dt = (t2.tt - t1.tt) / 2.0;
                 AstroTime tmid = t1.AddDays(dt);
@@ -5036,9 +5050,10 @@ $ASTRO_IAU_DATA()
         /// </returns>
         public static MoonQuarterInfo SearchMoonQuarter(AstroTime startTime)
         {
-            double angres = MoonPhase(startTime);
-            int quarter = (1 + (int)Math.Floor(angres / 90.0)) % 4;
-            AstroTime qtime = SearchMoonPhase(90.0 * quarter, startTime, 10.0);
+            double currentPhaseAngle = MoonPhase(startTime);
+            int quarter = (1 + (int)Math.Floor(currentPhaseAngle / 90.0)) % 4;
+            AstroTime qtime = SearchMoonPhase(90.0 * quarter, startTime, 10.0) ??
+                throw new InternalError($"Unable to find moon quarter {quarter} for startTime={startTime}.");
             return new MoonQuarterInfo(quarter, qtime);
         }
 
@@ -5055,15 +5070,15 @@ $ASTRO_IAU_DATA()
         /// <returns>The moon quarter that occurs next in time after the one passed in `mq`.</returns>
         public static MoonQuarterInfo NextMoonQuarter(MoonQuarterInfo mq)
         {
-            /* Skip 6 days past the previous found moon quarter to find the next one. */
-            /* This is less than the minimum possible increment. */
-            /* So far I have seen the interval well contained by the range (6.5, 8.3) days. */
+            // Skip 6 days past the previous found moon quarter to find the next one.
+            // This is less than the minimum possible increment.
+            // So far I have seen the interval well contained by the range (6.5, 8.3) days.
 
             AstroTime time = mq.time.AddDays(6.0);
             MoonQuarterInfo next_mq = SearchMoonQuarter(time);
             /* Verify that we found the expected moon quarter. */
             if (next_mq.quarter != (1 + mq.quarter) % 4)
-                throw new Exception("Internal error: found the wrong moon quarter.");
+                throw new InternalError("found the wrong moon quarter.");
             return next_mq;
         }
 
@@ -5096,7 +5111,7 @@ $ASTRO_IAU_DATA()
         /// </param>
         /// <returns>
         /// If successful, returns the date and time the moon reaches the phase specified by
-        /// `targetlon`. This function will return throw an exception if the phase does not
+        /// `targetlon`. This function will return `null` if the phase does not
         /// occur within `limitDays` of `startTime`; that is, if the search window is too small.
         /// </returns>
         public static AstroTime SearchMoonPhase(double targetLon, AstroTime startTime, double limitDays)
@@ -5111,7 +5126,6 @@ $ASTRO_IAU_DATA()
                 I have seen more than 0.9 days away from the simple prediction.
                 To be safe, we take the predicted time of the event and search
                 +/-1.5 days around it (a 3-day wide window).
-                Return null if the final result goes beyond limitDays after startTime.
             */
 
             const double uncertainty = 1.5;
@@ -5128,10 +5142,7 @@ $ASTRO_IAU_DATA()
                 dt2 = limitDays;
             AstroTime t1 = startTime.AddDays(dt1);
             AstroTime t2 = startTime.AddDays(dt2);
-            AstroTime time = Search(moon_offset, t1, t2, 1.0);
-            if (time == null)
-                throw new Exception(string.Format("Could not find moon longitude {0} within {1} days of {2}", targetLon, limitDays, startTime));
-            return time;
+            return Search(moon_offset, t1, t2, 1.0);
         }
 
         private static AstroTime InternalSearchAltitude(
@@ -5529,7 +5540,7 @@ $ASTRO_IAU_DATA()
                 }
             }
 
-            throw new Exception("Relative longitude search failed to converge.");
+            throw new InternalError("Relative longitude search failed to converge.");
         }
 
         private static double rlon_offset(Body body, AstroTime time, int direction, double targetRelLon)
@@ -5773,16 +5784,15 @@ $ASTRO_IAU_DATA()
                 /* Confirm the bracketing. */
                 double m1 = neg_elong_slope.Eval(t1);
                 if (m1 >= 0.0)
-                    throw new Exception("There is a bug in the bracketing algorithm! m1 = " + m1);
+                    throw new InternalError("There is a bug in the bracketing algorithm! m1 = " + m1);
 
                 double m2 = neg_elong_slope.Eval(t2);
                 if (m2 <= 0.0)
-                    throw new Exception("There is a bug in the bracketing algorithm! m2 = " + m2);
+                    throw new InternalError("There is a bug in the bracketing algorithm! m2 = " + m2);
 
                 /* Use the generic search algorithm to home in on where the slope crosses from negative to positive. */
-                AstroTime searchx = Search(neg_elong_slope, t1, t2, 10.0);
-                if (searchx == null)
-                    throw new Exception("Maximum elongation search failed.");
+                AstroTime searchx = Search(neg_elong_slope, t1, t2, 10.0) ??
+                    throw new InternalError("Maximum elongation search failed.");
 
                 if (searchx.tt >= startTime.tt)
                     return Elongation(body, searchx);
@@ -5793,7 +5803,7 @@ $ASTRO_IAU_DATA()
                 startTime = t2.AddDays(1.0);
             }
 
-            throw new Exception("Maximum elongation search iterated too many times.");
+            throw new InternalError("Maximum elongation search iterated too many times.");
         }
 
         ///
@@ -5846,7 +5856,7 @@ $ASTRO_IAU_DATA()
         {
             double r = a.Length() * b.Length();
             if (r < 1.0e-8)
-                throw new Exception("Cannot find angle between vectors because they are too short.");
+                throw new ArgumentException("Cannot find angle between vectors because they are too short.");
 
             double dot = (a.x*b.x + a.y*b.y + a.z*b.z) / r;
 
@@ -5928,11 +5938,11 @@ $ASTRO_IAU_DATA()
                     else
                     {
                         /* This should never happen. It should not be possible for both slopes to be zero. */
-                        throw new Exception("Internal error with slopes in SearchLunarApsis");
+                        throw new InternalError("both slopes are zero in SearchLunarApsis.");
                     }
 
                     if (search == null)
-                        throw new Exception("Failed to find slope transition in lunar apsis search.");
+                        throw new InternalError("Failed to find slope transition in lunar apsis search.");
 
                     double dist_au = SearchContext_MoonDistanceSlope.MoonDistance(search);
                     return new ApsisInfo(search, kind, dist_au);
@@ -5943,7 +5953,7 @@ $ASTRO_IAU_DATA()
             }
 
             /* It should not be possible to fail to find an apsis within 2 synodic months. */
-            throw new Exception("Internal error: should have found lunar apsis within 2 synodic months.");
+            throw new InternalError("should have found lunar apsis within 2 synodic months.");
         }
 
         /// <summary>
@@ -5973,7 +5983,7 @@ $ASTRO_IAU_DATA()
             AstroTime time = apsis.time.AddDays(skip);
             ApsisInfo next =  SearchLunarApsis(time);
             if ((int)next.kind + (int)apsis.kind != 1)
-                throw new Exception(string.Format("Internal error: previous apsis was {0}, but found {1} for next apsis.", apsis.kind, next.kind));
+                throw new InternalError($"Internal error: previous apsis was {apsis.kind}, but found {next.kind} for next apsis.");
             return next;
         }
 
@@ -6096,7 +6106,7 @@ $ASTRO_IAU_DATA()
             if (aphelion.time.tt >= startTime.tt)
                 return aphelion;
 
-            throw new Exception("Internal error: failed to find planet apsis.");
+            throw new InternalError("failed to find planet apsis.");
         }
 
 
@@ -6171,12 +6181,11 @@ $ASTRO_IAU_DATA()
                     else
                     {
                         /* This should never happen. It should not be possible for both slopes to be zero. */
-                        throw new Exception("Internal error with slopes in SearchPlanetApsis");
+                        throw new InternalError("Both slopes were zero in SearchPlanetApsis");
                     }
 
-                    AstroTime search = Search(slope_func, t1, t2, 1.0);
-                    if (search == null)
-                        throw new Exception("Failed to find slope transition in planetary apsis search.");
+                    AstroTime search = Search(slope_func, t1, t2, 1.0) ??
+                        throw new InternalError("Failed to find slope transition in planetary apsis search.");
 
                     double dist = HelioDistance(body, search);
                     return new ApsisInfo(search, kind, dist);
@@ -6186,7 +6195,7 @@ $ASTRO_IAU_DATA()
                 m1 = m2;
             }
             /* It should not be possible to fail to find an apsis within 2 planet orbits. */
-            throw new Exception("Internal error: should have found planetary apsis within 2 orbital periods.");
+            throw new InternalError("should have found planetary apsis within 2 orbital periods.");
         }
 
         /// <summary>
@@ -6224,7 +6233,7 @@ $ASTRO_IAU_DATA()
 
             /* Verify that we found the opposite apsis from the previous one. */
             if ((int)next.kind + (int)apsis.kind != 1)
-                throw new Exception(string.Format("Internal error: previous apsis was {0}, but found {1} for next apsis.", apsis.kind, next.kind));
+                throw new InternalError($"previous apsis was {apsis.kind}, but found {next.kind} for next apsis.");
 
             return next;
         }
@@ -6268,7 +6277,8 @@ $ASTRO_IAU_DATA()
             for (int fmcount=0; fmcount < 12; ++fmcount)
             {
                 // Search for the next full moon. Any eclipse will be near it.
-                AstroTime fullmoon = SearchMoonPhase(180.0, fmtime, 40.0);
+                AstroTime fullmoon = SearchMoonPhase(180.0, fmtime, 40.0) ??
+                    throw new InternalError("Failed to find next full moon.");
 
                 /*
                     Pruning: if the full Moon's ecliptic latitude is too large,
@@ -6313,7 +6323,7 @@ $ASTRO_IAU_DATA()
             }
 
             // This should never happen, because there should be at least 2 lunar eclipses per year.
-            throw new Exception("Internal error: failed to find lunar eclipse within 12 full moons.");
+            throw new InternalError("failed to find lunar eclipse within 12 full moons.");
         }
 
 
@@ -6373,7 +6383,8 @@ $ASTRO_IAU_DATA()
             for (int nmcount=0; nmcount < 12; ++nmcount)
             {
                 /* Search for the next new moon. Any eclipse will be near it. */
-                AstroTime newmoon = SearchMoonPhase(0.0, nmtime, 40.0);
+                AstroTime newmoon = SearchMoonPhase(0.0, nmtime, 40.0) ??
+                    throw new InternalError("Failed to find next new moon.");
 
                 /* Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible. */
                 double eclip_lat = MoonEclipticLatitudeDegrees(newmoon);
@@ -6396,7 +6407,7 @@ $ASTRO_IAU_DATA()
 
             /* Safety valve to prevent infinite loop. */
             /* This should never happen, because at least 2 solar eclipses happen per year. */
-            throw new Exception("Failure to find global solar eclipse.");
+            throw new InternalError("Failure to find global solar eclipse.");
         }
 
 
@@ -6512,7 +6523,7 @@ $ASTRO_IAU_DATA()
                 /* 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 || surface.r < 0.0)
-                    throw new Exception("Invalid surface distance from intersection.");
+                    throw new InternalError("Invalid surface distance from intersection.");
 
                 eclipse.kind = EclipseKindFromUmbra(surface.k);
             }
@@ -6675,7 +6686,8 @@ $ASTRO_IAU_DATA()
             for(;;)
             {
                 /* Search for the next new moon. Any eclipse will be near it. */
-                AstroTime newmoon = SearchMoonPhase(0.0, nmtime, 40.0);
+                AstroTime newmoon = SearchMoonPhase(0.0, nmtime, 40.0) ??
+                    throw new InternalError("Failed to find next new moon.");
 
                 /* Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible. */
                 double eclip_lat = MoonEclipticLatitudeDegrees(newmoon);
@@ -6777,9 +6789,8 @@ $ASTRO_IAU_DATA()
             AstroTime t2)
         {
             var context = new SearchContext_LocalEclipseTransition(func, direction, observer);
-            AstroTime search = Search(context, t1, t2, 1.0);
-            if (search == null)
-                throw new Exception("Local eclipse transition search failed.");
+            AstroTime search = Search(context, t1, t2, 1.0) ??
+                throw new InternalError("Local eclipse transition search failed.");
             return CalcEvent(observer, search);
         }
 
@@ -6808,9 +6819,8 @@ $ASTRO_IAU_DATA()
         {
             /* Search for the time the planet's penumbra begins/ends making contact with the center of the Earth. */
             var context = new SearchContext_PlanetShadowBoundary(body, planet_radius_km, direction);
-            AstroTime time = Search(context, t1, t2, 1.0);
-            if (time == null)
-                throw new Exception("Planet transit boundary search failed");
+            AstroTime time = Search(context, t1, t2, 1.0) ??
+                throw new InternalError("Planet transit boundary search failed");
             return time;
         }
 
@@ -6963,9 +6973,8 @@ $ASTRO_IAU_DATA()
                         context.Direction = -1;
                         kind = NodeEventKind.Descending;
                     }
-                    AstroTime result = Search(context, time1, time2, 1.0);
-                    if (result == null)
-                        throw new Exception("Could not find Moon node.");   // should never happen
+                    AstroTime result = Search(context, time1, time2, 1.0) ??
+                        throw new InternalError("Could not find Moon node.");
 
                     return new NodeEventInfo { time = result, kind = kind };
                 }
@@ -6992,16 +7001,16 @@ $ASTRO_IAU_DATA()
             {
                 case NodeEventKind.Ascending:
                     if (node.kind != NodeEventKind.Descending)
-                        throw new Exception("Internal error: previous node was ascending, but this node was: " + node.kind);
+                        throw new InternalError("previous node was ascending, but this node was: " + node.kind);
                     break;
 
                 case NodeEventKind.Descending:
                     if (node.kind != NodeEventKind.Ascending)
-                        throw new Exception("Internal error: previous node was descending, but this node was: " + node.kind);
+                        throw new InternalError("previous node was descending, but this node was: " + node.kind);
                     break;
 
                 default:
-                    throw new Exception("Previous node has an invalid node kind.");
+                    throw new ArgumentException("Previous node has an invalid node kind.");
             }
             return node;
         }
@@ -7272,16 +7281,15 @@ $ASTRO_IAU_DATA()
                 /* Confirm the bracketing. */
                 double m1 = mag_slope.Eval(t1);
                 if (m1 >= 0.0)
-                    throw new Exception("Internal error: m1 >= 0");    /* should never happen! */
+                    throw new InternalError("m1 >= 0");    /* should never happen! */
 
                 double m2 = mag_slope.Eval(t2);
                 if (m2 <= 0.0)
-                    throw new Exception("Internal error: m2 <= 0");    /* should never happen! */
+                    throw new InternalError("m2 <= 0");    /* should never happen! */
 
                 /* Use the generic search algorithm to home in on where the slope crosses from negative to positive. */
-                AstroTime tx = Search(mag_slope, t1, t2, 10.0);
-                if (tx == null)
-                    throw new Exception("Failed to find magnitude slope transition.");
+                AstroTime tx = Search(mag_slope, t1, t2, 10.0) ??
+                    throw new InternalError("Failed to find magnitude slope transition.");
 
                 if (tx.tt >= startTime.tt)
                     return Illumination(body, tx);
@@ -7292,7 +7300,7 @@ $ASTRO_IAU_DATA()
                 startTime = t2.AddDays(1.0);
             }
             // This should never happen. If it does, please report as a bug in Astronomy Engine.
-            throw new Exception("Peak magnitude search failed.");
+            throw new InternalError("Peak magnitude search failed.");
         }
 
         /// <summary>
@@ -8742,7 +8750,7 @@ $ASTRO_IAU_DATA()
                     return new ConstellationInfo(ConstelNames[b.index].symbol, ConstelNames[b.index].name, equ1875.ra, equ1875.dec);
 
             // This should never happen!
-            throw new Exception($"Unable to find constellation for coordinates: RA={ra}, DEC={dec}");
+            throw new InternalError($"Unable to find constellation for coordinates: RA={ra}, DEC={dec}");
         }
 
 $ASTRO_CONSTEL()
diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py
index 8f9486a1..40713d65 100644
--- a/generate/template/astronomy.py
+++ b/generate/template/astronomy.py
@@ -378,7 +378,7 @@ class InternalError(Error):
     Astronomy Engine for everyone! (Thank you in advance from the author.)
     """
     def __init__(self):
-        Error.__init__(self, 'Internal error - please report issue at https://github.com/cosinekitty/astronomy/issues')
+        Error.__init__(self, 'Internal error - please report issue, including stack trace, at https://github.com/cosinekitty/astronomy/issues')
 
 class NoConvergeError(Error):
     """A numeric solver did not converge.
@@ -6174,7 +6174,7 @@ def SearchLunarEclipse(startTime):
         # Search for the next full moon. Any eclipse will be near it.
         fullmoon = SearchMoonPhase(180, fmtime, 40)
         if fullmoon is None:
-            raise Error('Cannot find full moon.')
+            raise InternalError()   # should have always found the next full moon
 
         # Pruning: if the full Moon's ecliptic latitude is too large,
         # a lunar eclipse is not possible. Avoid needless work searching for
@@ -6258,7 +6258,7 @@ def SearchGlobalSolarEclipse(startTime):
         # 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')
+            raise InternalError()   # should always find the next new moon
 
         # Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible.
         eclip_lat = _MoonEclipticLatitudeDegrees(newmoon)
@@ -6333,6 +6333,8 @@ def SearchLocalSolarEclipse(startTime, observer):
     while True:
         # Search for the next new moon. Any eclipse will be near it.
         newmoon = SearchMoonPhase(0.0, nmtime, 40.0)
+        if newmoon is None:
+            raise InternalError()   # should always find the next new moon
 
         # Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible.
         eclip_lat = _MoonEclipticLatitudeDegrees(newmoon)
@@ -6587,7 +6589,7 @@ def SearchMoonNode(startTime):
             kind = NodeEventKind.Ascending if (eclip2.lat > eclip1.lat) else NodeEventKind.Descending
             result = Search(_MoonNodeSearchFunc, kind.value, time1, time2, 1.0)
             if result is None:
-                raise InternalError()
+                raise InternalError()   # should always find the next lunar node
             return NodeEventInfo(kind, result)
         time1 = time2
         eclip1 = eclip2
diff --git a/source/csharp/README.md b/source/csharp/README.md
index 0cbe96c7..edda7ed9 100644
--- a/source/csharp/README.md
+++ b/source/csharp/README.md
@@ -1740,7 +1740,7 @@ This function is useful for finding general phase angles outside those four quar
 | [`AstroTime`](#AstroTime) | `startTime` | The beginning of the time window in which to search for the Moon reaching the specified phase. |
 | `double` | `limitDays` | The number of days after `startTime` that limits the time window for the search. |
 
-**Returns:** If successful, returns the date and time the moon reaches the phase specified by `targetlon`. This function will return throw an exception if the phase does not occur within `limitDays` of `startTime`; that is, if the search window is too small.
+**Returns:** If successful, returns the date and time the moon reaches the phase specified by `targetlon`. This function will return `null` if the phase does not occur within `limitDays` of `startTime`; that is, if the search window is too small.
 
 <a name="Astronomy.SearchMoonQuarter"></a>
 ### Astronomy.SearchMoonQuarter(startTime) &#8658; [`MoonQuarterInfo`](#MoonQuarterInfo)
@@ -2528,6 +2528,15 @@ to report the visual magnitude and illuminated fraction of a celestial body at a
 
 ---
 
+<a name="InternalError"></a>
+## `class InternalError`
+
+**This exception indicates an unexpected error occurred inside Astronomy Engine.
+Please report any such errors by creating an issue at:
+https://github.com/cosinekitty/astronomy/issues**
+
+---
+
 <a name="InvalidBodyException"></a>
 ## `class InvalidBodyException`
 
diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs
index dac6e2d2..1afdb2b1 100644
--- a/source/csharp/astronomy.cs
+++ b/source/csharp/astronomy.cs
@@ -51,7 +51,20 @@ namespace CosineKitty
     {
         /// <summary>Creates an exception indicating that the given body is not valid for this operation.</summary>
         public InvalidBodyException(Body body):
-            base(string.Format("Invalid body: {0}", body))
+            base("Invalid body: " + body)
+            {}
+    }
+
+    /// <summary>
+    /// This exception indicates an unexpected error occurred inside Astronomy Engine.
+    /// Please report any such errors by creating an issue at:
+    /// https://github.com/cosinekitty/astronomy/issues
+    /// </summary>
+    public class InternalError: Exception
+    {
+        /// <summary>Creates an exception indicating that an unexpected error ocurred.</summary>
+        public InternalError(string message):
+            base("Internal error. Please report an issue at: https://github.com/cosinekitty/astronomy/issues. Diagnostic: " + message)
             {}
     }
 
@@ -5213,7 +5226,7 @@ namespace CosineKitty
 
                     ltime = ltime2;
                 }
-                throw new Exception("Light travel time correction did not converge");
+                throw new InternalError("Light travel time correction did not converge");
             }
         }
 
@@ -5912,7 +5925,8 @@ namespace CosineKitty
         private static AstroTime FindSeasonChange(double targetLon, int year, int month, int day)
         {
             var startTime = new AstroTime(year, month, day, 0, 0, 0);
-            return SearchSunLongitude(targetLon, startTime, 4.0);
+            return SearchSunLongitude(targetLon, startTime, 4.0) ??
+                throw new InternalError($"Cannot find solution for Sun longitude {targetLon}");
         }
 
         /// <summary>
@@ -6035,7 +6049,7 @@ namespace CosineKitty
             for(;;)
             {
                 if (++iter > iter_limit)
-                    throw new Exception(string.Format("Search did not converge within {0} iterations.", iter_limit));
+                    throw new InternalError(string.Format("Search did not converge within {0} iterations.", iter_limit));
 
                 double dt = (t2.tt - t1.tt) / 2.0;
                 AstroTime tmid = t1.AddDays(dt);
@@ -6248,9 +6262,10 @@ namespace CosineKitty
         /// </returns>
         public static MoonQuarterInfo SearchMoonQuarter(AstroTime startTime)
         {
-            double angres = MoonPhase(startTime);
-            int quarter = (1 + (int)Math.Floor(angres / 90.0)) % 4;
-            AstroTime qtime = SearchMoonPhase(90.0 * quarter, startTime, 10.0);
+            double currentPhaseAngle = MoonPhase(startTime);
+            int quarter = (1 + (int)Math.Floor(currentPhaseAngle / 90.0)) % 4;
+            AstroTime qtime = SearchMoonPhase(90.0 * quarter, startTime, 10.0) ??
+                throw new InternalError($"Unable to find moon quarter {quarter} for startTime={startTime}.");
             return new MoonQuarterInfo(quarter, qtime);
         }
 
@@ -6267,15 +6282,15 @@ namespace CosineKitty
         /// <returns>The moon quarter that occurs next in time after the one passed in `mq`.</returns>
         public static MoonQuarterInfo NextMoonQuarter(MoonQuarterInfo mq)
         {
-            /* Skip 6 days past the previous found moon quarter to find the next one. */
-            /* This is less than the minimum possible increment. */
-            /* So far I have seen the interval well contained by the range (6.5, 8.3) days. */
+            // Skip 6 days past the previous found moon quarter to find the next one.
+            // This is less than the minimum possible increment.
+            // So far I have seen the interval well contained by the range (6.5, 8.3) days.
 
             AstroTime time = mq.time.AddDays(6.0);
             MoonQuarterInfo next_mq = SearchMoonQuarter(time);
             /* Verify that we found the expected moon quarter. */
             if (next_mq.quarter != (1 + mq.quarter) % 4)
-                throw new Exception("Internal error: found the wrong moon quarter.");
+                throw new InternalError("found the wrong moon quarter.");
             return next_mq;
         }
 
@@ -6308,7 +6323,7 @@ namespace CosineKitty
         /// </param>
         /// <returns>
         /// If successful, returns the date and time the moon reaches the phase specified by
-        /// `targetlon`. This function will return throw an exception if the phase does not
+        /// `targetlon`. This function will return `null` if the phase does not
         /// occur within `limitDays` of `startTime`; that is, if the search window is too small.
         /// </returns>
         public static AstroTime SearchMoonPhase(double targetLon, AstroTime startTime, double limitDays)
@@ -6323,7 +6338,6 @@ namespace CosineKitty
                 I have seen more than 0.9 days away from the simple prediction.
                 To be safe, we take the predicted time of the event and search
                 +/-1.5 days around it (a 3-day wide window).
-                Return null if the final result goes beyond limitDays after startTime.
             */
 
             const double uncertainty = 1.5;
@@ -6340,10 +6354,7 @@ namespace CosineKitty
                 dt2 = limitDays;
             AstroTime t1 = startTime.AddDays(dt1);
             AstroTime t2 = startTime.AddDays(dt2);
-            AstroTime time = Search(moon_offset, t1, t2, 1.0);
-            if (time == null)
-                throw new Exception(string.Format("Could not find moon longitude {0} within {1} days of {2}", targetLon, limitDays, startTime));
-            return time;
+            return Search(moon_offset, t1, t2, 1.0);
         }
 
         private static AstroTime InternalSearchAltitude(
@@ -6741,7 +6752,7 @@ namespace CosineKitty
                 }
             }
 
-            throw new Exception("Relative longitude search failed to converge.");
+            throw new InternalError("Relative longitude search failed to converge.");
         }
 
         private static double rlon_offset(Body body, AstroTime time, int direction, double targetRelLon)
@@ -6985,16 +6996,15 @@ namespace CosineKitty
                 /* Confirm the bracketing. */
                 double m1 = neg_elong_slope.Eval(t1);
                 if (m1 >= 0.0)
-                    throw new Exception("There is a bug in the bracketing algorithm! m1 = " + m1);
+                    throw new InternalError("There is a bug in the bracketing algorithm! m1 = " + m1);
 
                 double m2 = neg_elong_slope.Eval(t2);
                 if (m2 <= 0.0)
-                    throw new Exception("There is a bug in the bracketing algorithm! m2 = " + m2);
+                    throw new InternalError("There is a bug in the bracketing algorithm! m2 = " + m2);
 
                 /* Use the generic search algorithm to home in on where the slope crosses from negative to positive. */
-                AstroTime searchx = Search(neg_elong_slope, t1, t2, 10.0);
-                if (searchx == null)
-                    throw new Exception("Maximum elongation search failed.");
+                AstroTime searchx = Search(neg_elong_slope, t1, t2, 10.0) ??
+                    throw new InternalError("Maximum elongation search failed.");
 
                 if (searchx.tt >= startTime.tt)
                     return Elongation(body, searchx);
@@ -7005,7 +7015,7 @@ namespace CosineKitty
                 startTime = t2.AddDays(1.0);
             }
 
-            throw new Exception("Maximum elongation search iterated too many times.");
+            throw new InternalError("Maximum elongation search iterated too many times.");
         }
 
         ///
@@ -7058,7 +7068,7 @@ namespace CosineKitty
         {
             double r = a.Length() * b.Length();
             if (r < 1.0e-8)
-                throw new Exception("Cannot find angle between vectors because they are too short.");
+                throw new ArgumentException("Cannot find angle between vectors because they are too short.");
 
             double dot = (a.x*b.x + a.y*b.y + a.z*b.z) / r;
 
@@ -7140,11 +7150,11 @@ namespace CosineKitty
                     else
                     {
                         /* This should never happen. It should not be possible for both slopes to be zero. */
-                        throw new Exception("Internal error with slopes in SearchLunarApsis");
+                        throw new InternalError("both slopes are zero in SearchLunarApsis.");
                     }
 
                     if (search == null)
-                        throw new Exception("Failed to find slope transition in lunar apsis search.");
+                        throw new InternalError("Failed to find slope transition in lunar apsis search.");
 
                     double dist_au = SearchContext_MoonDistanceSlope.MoonDistance(search);
                     return new ApsisInfo(search, kind, dist_au);
@@ -7155,7 +7165,7 @@ namespace CosineKitty
             }
 
             /* It should not be possible to fail to find an apsis within 2 synodic months. */
-            throw new Exception("Internal error: should have found lunar apsis within 2 synodic months.");
+            throw new InternalError("should have found lunar apsis within 2 synodic months.");
         }
 
         /// <summary>
@@ -7185,7 +7195,7 @@ namespace CosineKitty
             AstroTime time = apsis.time.AddDays(skip);
             ApsisInfo next =  SearchLunarApsis(time);
             if ((int)next.kind + (int)apsis.kind != 1)
-                throw new Exception(string.Format("Internal error: previous apsis was {0}, but found {1} for next apsis.", apsis.kind, next.kind));
+                throw new InternalError($"Internal error: previous apsis was {apsis.kind}, but found {next.kind} for next apsis.");
             return next;
         }
 
@@ -7308,7 +7318,7 @@ namespace CosineKitty
             if (aphelion.time.tt >= startTime.tt)
                 return aphelion;
 
-            throw new Exception("Internal error: failed to find planet apsis.");
+            throw new InternalError("failed to find planet apsis.");
         }
 
 
@@ -7383,12 +7393,11 @@ namespace CosineKitty
                     else
                     {
                         /* This should never happen. It should not be possible for both slopes to be zero. */
-                        throw new Exception("Internal error with slopes in SearchPlanetApsis");
+                        throw new InternalError("Both slopes were zero in SearchPlanetApsis");
                     }
 
-                    AstroTime search = Search(slope_func, t1, t2, 1.0);
-                    if (search == null)
-                        throw new Exception("Failed to find slope transition in planetary apsis search.");
+                    AstroTime search = Search(slope_func, t1, t2, 1.0) ??
+                        throw new InternalError("Failed to find slope transition in planetary apsis search.");
 
                     double dist = HelioDistance(body, search);
                     return new ApsisInfo(search, kind, dist);
@@ -7398,7 +7407,7 @@ namespace CosineKitty
                 m1 = m2;
             }
             /* It should not be possible to fail to find an apsis within 2 planet orbits. */
-            throw new Exception("Internal error: should have found planetary apsis within 2 orbital periods.");
+            throw new InternalError("should have found planetary apsis within 2 orbital periods.");
         }
 
         /// <summary>
@@ -7436,7 +7445,7 @@ namespace CosineKitty
 
             /* Verify that we found the opposite apsis from the previous one. */
             if ((int)next.kind + (int)apsis.kind != 1)
-                throw new Exception(string.Format("Internal error: previous apsis was {0}, but found {1} for next apsis.", apsis.kind, next.kind));
+                throw new InternalError($"previous apsis was {apsis.kind}, but found {next.kind} for next apsis.");
 
             return next;
         }
@@ -7480,7 +7489,8 @@ namespace CosineKitty
             for (int fmcount=0; fmcount < 12; ++fmcount)
             {
                 // Search for the next full moon. Any eclipse will be near it.
-                AstroTime fullmoon = SearchMoonPhase(180.0, fmtime, 40.0);
+                AstroTime fullmoon = SearchMoonPhase(180.0, fmtime, 40.0) ??
+                    throw new InternalError("Failed to find next full moon.");
 
                 /*
                     Pruning: if the full Moon's ecliptic latitude is too large,
@@ -7525,7 +7535,7 @@ namespace CosineKitty
             }
 
             // This should never happen, because there should be at least 2 lunar eclipses per year.
-            throw new Exception("Internal error: failed to find lunar eclipse within 12 full moons.");
+            throw new InternalError("failed to find lunar eclipse within 12 full moons.");
         }
 
 
@@ -7585,7 +7595,8 @@ namespace CosineKitty
             for (int nmcount=0; nmcount < 12; ++nmcount)
             {
                 /* Search for the next new moon. Any eclipse will be near it. */
-                AstroTime newmoon = SearchMoonPhase(0.0, nmtime, 40.0);
+                AstroTime newmoon = SearchMoonPhase(0.0, nmtime, 40.0) ??
+                    throw new InternalError("Failed to find next new moon.");
 
                 /* Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible. */
                 double eclip_lat = MoonEclipticLatitudeDegrees(newmoon);
@@ -7608,7 +7619,7 @@ namespace CosineKitty
 
             /* Safety valve to prevent infinite loop. */
             /* This should never happen, because at least 2 solar eclipses happen per year. */
-            throw new Exception("Failure to find global solar eclipse.");
+            throw new InternalError("Failure to find global solar eclipse.");
         }
 
 
@@ -7724,7 +7735,7 @@ namespace CosineKitty
                 /* 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 || surface.r < 0.0)
-                    throw new Exception("Invalid surface distance from intersection.");
+                    throw new InternalError("Invalid surface distance from intersection.");
 
                 eclipse.kind = EclipseKindFromUmbra(surface.k);
             }
@@ -7887,7 +7898,8 @@ namespace CosineKitty
             for(;;)
             {
                 /* Search for the next new moon. Any eclipse will be near it. */
-                AstroTime newmoon = SearchMoonPhase(0.0, nmtime, 40.0);
+                AstroTime newmoon = SearchMoonPhase(0.0, nmtime, 40.0) ??
+                    throw new InternalError("Failed to find next new moon.");
 
                 /* Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible. */
                 double eclip_lat = MoonEclipticLatitudeDegrees(newmoon);
@@ -7989,9 +8001,8 @@ namespace CosineKitty
             AstroTime t2)
         {
             var context = new SearchContext_LocalEclipseTransition(func, direction, observer);
-            AstroTime search = Search(context, t1, t2, 1.0);
-            if (search == null)
-                throw new Exception("Local eclipse transition search failed.");
+            AstroTime search = Search(context, t1, t2, 1.0) ??
+                throw new InternalError("Local eclipse transition search failed.");
             return CalcEvent(observer, search);
         }
 
@@ -8020,9 +8031,8 @@ namespace CosineKitty
         {
             /* Search for the time the planet's penumbra begins/ends making contact with the center of the Earth. */
             var context = new SearchContext_PlanetShadowBoundary(body, planet_radius_km, direction);
-            AstroTime time = Search(context, t1, t2, 1.0);
-            if (time == null)
-                throw new Exception("Planet transit boundary search failed");
+            AstroTime time = Search(context, t1, t2, 1.0) ??
+                throw new InternalError("Planet transit boundary search failed");
             return time;
         }
 
@@ -8175,9 +8185,8 @@ namespace CosineKitty
                         context.Direction = -1;
                         kind = NodeEventKind.Descending;
                     }
-                    AstroTime result = Search(context, time1, time2, 1.0);
-                    if (result == null)
-                        throw new Exception("Could not find Moon node.");   // should never happen
+                    AstroTime result = Search(context, time1, time2, 1.0) ??
+                        throw new InternalError("Could not find Moon node.");
 
                     return new NodeEventInfo { time = result, kind = kind };
                 }
@@ -8204,16 +8213,16 @@ namespace CosineKitty
             {
                 case NodeEventKind.Ascending:
                     if (node.kind != NodeEventKind.Descending)
-                        throw new Exception("Internal error: previous node was ascending, but this node was: " + node.kind);
+                        throw new InternalError("previous node was ascending, but this node was: " + node.kind);
                     break;
 
                 case NodeEventKind.Descending:
                     if (node.kind != NodeEventKind.Ascending)
-                        throw new Exception("Internal error: previous node was descending, but this node was: " + node.kind);
+                        throw new InternalError("previous node was descending, but this node was: " + node.kind);
                     break;
 
                 default:
-                    throw new Exception("Previous node has an invalid node kind.");
+                    throw new ArgumentException("Previous node has an invalid node kind.");
             }
             return node;
         }
@@ -8484,16 +8493,15 @@ namespace CosineKitty
                 /* Confirm the bracketing. */
                 double m1 = mag_slope.Eval(t1);
                 if (m1 >= 0.0)
-                    throw new Exception("Internal error: m1 >= 0");    /* should never happen! */
+                    throw new InternalError("m1 >= 0");    /* should never happen! */
 
                 double m2 = mag_slope.Eval(t2);
                 if (m2 <= 0.0)
-                    throw new Exception("Internal error: m2 <= 0");    /* should never happen! */
+                    throw new InternalError("m2 <= 0");    /* should never happen! */
 
                 /* Use the generic search algorithm to home in on where the slope crosses from negative to positive. */
-                AstroTime tx = Search(mag_slope, t1, t2, 10.0);
-                if (tx == null)
-                    throw new Exception("Failed to find magnitude slope transition.");
+                AstroTime tx = Search(mag_slope, t1, t2, 10.0) ??
+                    throw new InternalError("Failed to find magnitude slope transition.");
 
                 if (tx.tt >= startTime.tt)
                     return Illumination(body, tx);
@@ -8504,7 +8512,7 @@ namespace CosineKitty
                 startTime = t2.AddDays(1.0);
             }
             // This should never happen. If it does, please report as a bug in Astronomy Engine.
-            throw new Exception("Peak magnitude search failed.");
+            throw new InternalError("Peak magnitude search failed.");
         }
 
         /// <summary>
@@ -9954,7 +9962,7 @@ namespace CosineKitty
                     return new ConstellationInfo(ConstelNames[b.index].symbol, ConstelNames[b.index].name, equ1875.ra, equ1875.dec);
 
             // This should never happen!
-            throw new Exception($"Unable to find constellation for coordinates: RA={ra}, DEC={dec}");
+            throw new InternalError($"Unable to find constellation for coordinates: RA={ra}, DEC={dec}");
         }
 
         private static readonly constel_info_t[] ConstelNames = new constel_info_t[]
diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py
index 0a01fe76..616504e7 100644
--- a/source/python/astronomy/astronomy.py
+++ b/source/python/astronomy/astronomy.py
@@ -378,7 +378,7 @@ class InternalError(Error):
     Astronomy Engine for everyone! (Thank you in advance from the author.)
     """
     def __init__(self):
-        Error.__init__(self, 'Internal error - please report issue at https://github.com/cosinekitty/astronomy/issues')
+        Error.__init__(self, 'Internal error - please report issue, including stack trace, at https://github.com/cosinekitty/astronomy/issues')
 
 class NoConvergeError(Error):
     """A numeric solver did not converge.
@@ -8667,7 +8667,7 @@ def SearchLunarEclipse(startTime):
         # Search for the next full moon. Any eclipse will be near it.
         fullmoon = SearchMoonPhase(180, fmtime, 40)
         if fullmoon is None:
-            raise Error('Cannot find full moon.')
+            raise InternalError()   # should have always found the next full moon
 
         # Pruning: if the full Moon's ecliptic latitude is too large,
         # a lunar eclipse is not possible. Avoid needless work searching for
@@ -8751,7 +8751,7 @@ def SearchGlobalSolarEclipse(startTime):
         # 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')
+            raise InternalError()   # should always find the next new moon
 
         # Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible.
         eclip_lat = _MoonEclipticLatitudeDegrees(newmoon)
@@ -8826,6 +8826,8 @@ def SearchLocalSolarEclipse(startTime, observer):
     while True:
         # Search for the next new moon. Any eclipse will be near it.
         newmoon = SearchMoonPhase(0.0, nmtime, 40.0)
+        if newmoon is None:
+            raise InternalError()   # should always find the next new moon
 
         # Pruning: if the new moon's ecliptic latitude is too large, a solar eclipse is not possible.
         eclip_lat = _MoonEclipticLatitudeDegrees(newmoon)
@@ -9080,7 +9082,7 @@ def SearchMoonNode(startTime):
             kind = NodeEventKind.Ascending if (eclip2.lat > eclip1.lat) else NodeEventKind.Descending
             result = Search(_MoonNodeSearchFunc, kind.value, time1, time2, 1.0)
             if result is None:
-                raise InternalError()
+                raise InternalError()   # should always find the next lunar node
             return NodeEventInfo(kind, result)
         time1 = time2
         eclip1 = eclip2