diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs
index ce28eb7d..91051e48 100644
--- a/generate/template/astronomy.cs
+++ b/generate/template/astronomy.cs
@@ -3741,8 +3741,7 @@ $ASTRO_CSHARP_CHEBYSHEV(8);
 
             for (i=0; i < npoints; ++i)
             {
-                double ut = t1.ut + (i * interval);
-                AstroTime time = new AstroTime(ut);
+                AstroTime time = t1.AddDays(i * interval);
                 double dist = NeptuneHelioDistance(time);
                 if (i == 0)
                 {
@@ -3837,7 +3836,7 @@ $ASTRO_CSHARP_CHEBYSHEV(8);
                 {
                     /* There is a change of slope polarity within the time range [t1, t2]. */
                     /* Therefore this time range contains an apsis. */
-                    /* Figure out whether it is perigee or apogee. */
+                    /* Figure out whether it is perihelion or aphelion. */
 
                     SearchContext_PlanetDistanceSlope slope_func;
                     ApsisKind kind;
diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py
index 0d6f958e..57cba722 100644
--- a/generate/template/astronomy.py
+++ b/generate/template/astronomy.py
@@ -53,6 +53,7 @@ _SECONDS_PER_DAY = 24.0 * 3600.0
 _SOLAR_DAYS_PER_SIDEREAL_DAY = 0.9972695717592592
 _MEAN_SYNODIC_MONTH = 29.530588
 _EARTH_ORBITAL_PERIOD = 365.256
+_NEPTUNE_ORBITAL_PERIOD = 60189.0
 _REFRACTION_NEAR_HORIZON = 34.0 / 60.0
 _SUN_RADIUS_AU  = 4.6505e-3
 _MOON_RADIUS_AU = 1.15717e-5
@@ -171,7 +172,7 @@ _PlanetOrbitalPeriod = [
     4332.589,
     10759.22,
     30685.4,
-    60189.0,
+    _NEPTUNE_ORBITAL_PERIOD,
     90560.0
 ]
 
@@ -200,6 +201,11 @@ class BadVectorError(Error):
     def __init__(self):
         Error.__init__(self, 'Vector is too small to have a direction.')
 
+class BadTimeError(Error):
+    """Cannot calculate Pluto position for this date/time."""
+    def __init__(self, time):
+        Error.__init__(self, 'Cannot calculate Pluto position for time: {}'.format(time))
+
 class InternalError(Error):
     """An internal error occured that should be reported as a bug.
 
@@ -966,16 +972,17 @@ _vsop = [
     $ASTRO_LIST_VSOP(Neptune),
 ]
 
+def _VsopFormula(formula, t):
+    tpower = 1.0
+    coord = 0.0
+    for series in formula:
+        coord += tpower * sum(A * math.cos(B + C*t) for (A, B, C) in series)
+        tpower *= t
+    return coord
+
 def _CalcVsop(model, time):
-    spher = []
     t = time.tt / 365250.0
-    for formula in model:
-        tpower = 1.0
-        coord = 0.0
-        for series in formula:
-            coord += tpower * sum(A * math.cos(B + C*t) for (A, B, C) in series)
-            tpower *= t
-        spher.append(coord)
+    spher = [_VsopFormula(formula, t) for formula in model]
 
     # Convert spherical coordinates to ecliptic cartesian coordinates.
     r_coslat = spher[2] * math.cos(spher[1])
@@ -989,6 +996,12 @@ def _CalcVsop(model, time):
     vz = 0.397776982902*ey + 0.917482137087*ez
     return Vector(vx, vy, vz, time)
 
+def _VsopHelioDistance(model, time):
+    # The caller only wants to know the distance between the planet and the Sun.
+    # So we only need to calculate the radial component of the spherical coordinates.
+    # There is no need to translate coordinates.
+    return _VsopFormula(model[2], time.tt / 365250.0)
+
 def _CalcEarth(time):
     return _CalcVsop(_vsop[Body.Earth], time)
 
@@ -1020,7 +1033,7 @@ def _CalcChebyshev(model, time):
                     p1 = p2
                 pos.append(sum - coeff[0][d]/2)
             return Vector(pos[0], pos[1], pos[2], time)
-    raise Error('Cannot extrapolate Chebyshev model for given Terrestrial Time: {}'.format(time.tt))
+    raise BadTimeError(time)
 
 # END CHEBYSHEV
 #----------------------------------------------------------------------------
@@ -1244,7 +1257,7 @@ def HelioVector(body, time):
     if body == Body.Pluto:
         return _CalcChebyshev(_pluto, time)
 
-    if 0 <= body <= len(_vsop):
+    if 0 <= body < len(_vsop):
         return _CalcVsop(_vsop[body], time)
 
     if body == Body.Sun:
@@ -1258,6 +1271,37 @@ def HelioVector(body, time):
     raise InvalidBodyError()
 
 
+def HelioDistance(body, time):
+    """Calculates the distance between a body and the Sun at a given time.
+
+    Given a date and time, this function calculates the distance between
+    the center of `body` and the center of the Sun.
+    For the planets Mercury through Neptune, this function is significantly
+    more efficient than calling #HelioVector followed by taking the length
+    of the resulting vector.
+
+    Parameters
+    ----------
+    body : Body
+        A body for which to calculate a heliocentric distance:
+        the Sun, Moon, or any of the planets.
+    time : Time
+        The date and time for which to calculate the heliocentric distance.
+
+    Returns
+    -------
+    float
+        The heliocentric distance in AU.
+    """
+    if body == Body.Sun:
+        return 0.0
+
+    if 0 <= body < len(_vsop):
+        return _VsopHelioDistance(_vsop[body], time)
+
+    return HelioVector(body, time).Length()
+
+
 def GeoVector(body, time, aberration):
     """Calculates geocentric Cartesian coordinates of a body in the J2000 equatorial system.
 
@@ -2921,7 +2965,7 @@ def Seasons(year):
 def _MoonDistance(time):
     return _CalcMoon(time).distance_au
 
-def _distance_slope(direction, time):
+def _moon_distance_slope(direction, time):
     dt = 0.001
     t1 = time.AddDays(-dt/2.0)
     t2 = time.AddDays(+dt/2.0)
@@ -3008,11 +3052,11 @@ def SearchLunarApsis(startTime):
     """
     increment = 5.0     # number of days to skip on each iteration
     t1 = startTime
-    m1 = _distance_slope(+1, t1)
+    m1 = _moon_distance_slope(+1, t1)
     iter = 0
     while iter * increment < 2.0 * _MEAN_SYNODIC_MONTH:
         t2 = t1.AddDays(increment)
-        m2 = _distance_slope(+1, t2)
+        m2 = _moon_distance_slope(+1, t2)
         if m1 * m2 <= 0.0:
             # There is a change of slope polarity within the time range [t1, t2].
             # Therefore this time range contains an apsis.
@@ -3020,12 +3064,12 @@ def SearchLunarApsis(startTime):
             if m1 < 0.0 or m2 > 0.0:
                 # We found a minimum-distance event: perigee.
                 # Search the time range for the time when the slope goes from negative to positive.
-                apsis_time = Search(_distance_slope, +1, t1, t2, 1.0)
+                apsis_time = Search(_moon_distance_slope, +1, t1, t2, 1.0)
                 kind = ApsisKind.Pericenter
             elif m1 > 0.0 or m2 < 0.0:
                 # We found a maximum-distance event: apogee.
                 # Search the time range for the time when the slope goes from positive to negative.
-                apsis_time = Search(_distance_slope, -1, t1, t2, 1.0)
+                apsis_time = Search(_moon_distance_slope, -1, t1, t2, 1.0)
                 kind = ApsisKind.Apocenter
             else:
                 # This should never happen. It should not be possible for both slopes to be zero.
@@ -3072,10 +3116,189 @@ def NextLunarApsis(apsis):
     if apsis.kind not in [ApsisKind.Apocenter, ApsisKind.Pericenter]:
         raise Error('Parameter "apsis" contains an invalid "kind" value.')
     if next.kind + apsis.kind != 1:
-        raise InternalError()
+        raise InternalError()   # should have found opposite apsis kind
     return next
 
 
+def _planet_distance_slope(context, time):
+    (direction, body) = context
+    dt = 0.001
+    t1 = time.AddDays(-dt/2.0)
+    t2 = time.AddDays(+dt/2.0)
+    dist1 = HelioDistance(body, t1)
+    dist2 = HelioDistance(body, t2)
+    return direction * (dist2 - dist1) / dt
+
+
+def SearchPlanetApsis(body, startTime):
+    """Finds the next planet perihelion or aphelion, after a given time.
+
+    Given a date and time to start the search in `startTime`, this function finds the
+    next date and time that the center of the specified planet reaches the closest or farthest point
+    in its orbit with respect to the center of the Sun, whichever comes first after `startTime`.
+
+    The closest point is called *perihelion* and the farthest point is called *aphelion*.
+    The word *apsis* refers to either event.
+
+    To iterate through consecutive alternating perihelion and aphelion events,
+    call `SearchPlanetApsis` once, then use the return value to call #NextPlanetApsis.
+    After that, keep feeding the previous return value from `NextPlanetApsis`
+    into another call of `NextPlanetApsis` as many times as desired.
+
+    Parameters
+    ----------
+    body : Body
+        The planet for which to find the next perihelion/aphelion event.
+        Not allowed to be `Body.Sun` or `Body.Moon`.
+    startTime : Time
+        The date and time at which to start searching for the next perihelion or aphelion.
+
+    Returns
+    -------
+    Apsis
+    """
+    if body == Body.Neptune:
+        return _SearchNeptuneApsis(startTime)
+    positive_slope = (+1.0, body)
+    negative_slope = (-1.0, body)
+    orbit_period_days = _PlanetOrbitalPeriod[body]
+    increment = orbit_period_days / 6.0
+    t1 = startTime
+    m1 = _planet_distance_slope(positive_slope, t1)
+    iter = 0
+    while iter * increment < 2 * orbit_period_days:
+        t2 = t1.AddDays(increment)
+        m2 = _planet_distance_slope(positive_slope, t2)
+        if m1 * m2 <= 0.0:
+            # There is a change of slope polarity within the time range [t1, t2].
+            # Therefore this time range contains an apsis.
+            # Figure out whether it is perihelion or aphelion.
+            if m1 < 0.0 or m2 > 0.0:
+                # We found a minimum-distance event: perihelion.
+                # Search the time range for the time when the slope goes from negative to positive.
+                context = positive_slope
+                kind = ApsisKind.Pericenter
+            elif m1 > 0.0 or m2 < 0.0:
+                # We found a maximum-distance event: aphelion.
+                # Search the time range for the time when the slope goes from positive to negative.
+                context = negative_slope
+                kind = ApsisKind.Apocenter
+            else:
+                raise InternalError()   # at least one of the planet distance slopes should have been nonzero
+            search = Search(_planet_distance_slope, context, t1, t2, 1.0)
+            if search is None:
+                raise InternalError()   # failed to find where planet distance slope passed through zero
+
+            dist = HelioDistance(body, search)
+            return Apsis(search, kind, dist)
+        t1 = t2
+        m1 = m2
+        iter += 1
+    raise InternalError()   # should have found planet apsis within 2 planet orbits
+
+
+def NextPlanetApsis(body, apsis):
+    """Finds the next planetary perihelion or aphelion event in a series.
+
+    This function requires an #Apsis value obtained from a call
+    to #SearchPlanetApsis or `NextPlanetApsis`.
+    Given an aphelion event, this function finds the next perihelion event, and vice versa.
+    See #SearchPlanetApsis for more details.
+
+    Parameters
+    ----------
+    body : Body
+        The planet for which to find the next perihelion/aphelion event.
+        Not allowed to be `Body.Sun` or `Body.Moon`.
+        Must match the body passed into the call that produced the `apsis` parameter.
+    apsis : Apsis
+        An apsis event obtained from a call to #SearchPlanetApsis or `NextPlanetApsis`.
+
+    Returns
+    -------
+    Apsis
+    """
+    if apsis.kind not in [ApsisKind.Apocenter, ApsisKind.Pericenter]:
+        raise Error('Parameter "apsis" contains an invalid "kind" value.')
+    # Skip 1/4 of an orbit before starting search again.
+    skip = 0.25 * _PlanetOrbitalPeriod[body]
+    time = apsis.time.AddDays(skip)
+    next = SearchPlanetApsis(body, time)
+    # Verify that we found the opposite apsis from the previous one.
+    if next.kind + apsis.kind != 1:
+        raise InternalError()   # should have found opposite planetary apsis type
+    return next
+
+def _NeptuneHelioDistance(time):
+    return _VsopHelioDistance(_vsop[Body.Neptune], time)
+
+def _NeptuneExtreme(kind, start_time, dayspan):
+    direction = +1.0 if (kind == ApsisKind.Apocenter) else -1.0
+    npoints = 10
+    while True:
+        interval = dayspan / (npoints - 1)
+        # Iterate until uncertainty is less than one minute.
+        if interval < 1/1440:
+            apsis_time = start_time.AddDays(interval/2)
+            dist_au = _NeptuneHelioDistance(apsis_time)
+            return Apsis(apsis_time, kind, dist_au)
+        for i in range(npoints):
+            time = start_time.AddDays(i * interval)
+            dist = direction * _NeptuneHelioDistance(time)
+            if i==0 or dist > best_dist:
+                best_i = i
+                best_dist = dist
+        # Narrow in on the extreme point.
+        start_time = start_time.AddDays((best_i - 1) * interval)
+        dayspan = 2 * interval
+
+def _SearchNeptuneApsis(startTime):
+    # Neptune is a special case for two reasons:
+    # 1. Its orbit is nearly circular (low orbital eccentricity).
+    # 2. It is so distant from the Sun that the orbital period is very long.
+    # Put together, this causes wobbling of the Sun around the Solar System Barycenter (SSB)
+    # to be so significant that there are 3 local minima in the distance-vs-time curve
+    # near each apsis. Therefore, unlike for other planets, we can't use an optimized
+    # algorithm for finding dr/dt = 0.
+    # Instead, we use a dumb, brute-force algorithm of sampling and finding min/max
+    # heliocentric distance.
+    #
+    # Rewind approximately 30 degrees in the orbit,
+    # then search forward for 270 degrees.
+    # This is a very cautious way to prevent missing an apsis.
+    # Typically we will find two apsides, and we pick whichever
+    # apsis is ealier, but after startTime.
+    # Sample points around this orbital arc and find when the distance
+    # is greatest and smallest.
+    t1 = startTime.AddDays(_NEPTUNE_ORBITAL_PERIOD * ( -30 / 360))
+    t2 = startTime.AddDays(_NEPTUNE_ORBITAL_PERIOD * (+270 / 360))
+    t_min = t_max = t1
+    npoints = 1000
+    interval = (t2.ut - t1.ut) / (npoints - 1)
+
+    for i in range(npoints):
+        time = t1.AddDays(i * interval)
+        dist = _NeptuneHelioDistance(time)
+        if i == 0:
+            min_dist = max_dist = dist
+        else:
+            if dist > max_dist:
+                max_dist = dist
+                t_max = time
+            if dist < min_dist:
+                min_dist = dist
+                t_min = time
+
+    perihelion = _NeptuneExtreme(ApsisKind.Pericenter, t_min.AddDays(-2*interval), 4*interval)
+    aphelion = _NeptuneExtreme(ApsisKind.Apocenter, t_max.AddDays(-2*interval), 4*interval)
+    if perihelion.time.tt >= startTime.tt:
+        if startTime.tt <= aphelion.time.tt < perihelion.time.tt:
+            return aphelion
+        return perihelion
+    if aphelion.time.tt >= startTime.tt:
+        return aphelion
+    raise InternalError()   # failed to find Neptune apsis
+
 def VectorFromSphere(sphere, time):
     """Converts spherical coordinates to Cartesian coordinates.
 
diff --git a/generate/test.py b/generate/test.py
index bd5d00a2..07b974b9 100755
--- a/generate/test.py
+++ b/generate/test.py
@@ -2,6 +2,7 @@
 import sys
 import math
 import re
+import os
 sys.path.append('../source/python')
 import astronomy
 
@@ -432,22 +433,10 @@ ElongTestData = [
     ( astronomy.Body.Venus,   "2018-08-07T17:02Z", "2018-08-17T17:02Z", 45.90, 'evening' )
 ]
 
-def ParseDate(text):
-    m = re.match(r'^(\d+)-(\d+)-(\d+)T(\d+):(\d+)Z$', text)
-    if not m:
-        print('ParseDate: invalid date text "{}"'.format(text))
-        raise Exception('Bad elongation test data')
-    year = int(m.group(1))
-    month = int(m.group(2))
-    day = int(m.group(3))
-    hour = int(m.group(4))
-    minute = int(m.group(5))
-    return astronomy.Time.Make(year, month, day, hour, minute, 0)
-
 def TestMaxElong(body, searchText, eventText, angle, visibility):
     name = body.name
-    searchTime = ParseDate(searchText)
-    eventTime = ParseDate(eventText)
+    searchTime = astronomy.Time.Parse(searchText)
+    eventTime = astronomy.Time.Parse(eventText)
     evt = astronomy.SearchMaxElongation(body, searchTime)
     if evt is None:
         print('TestMaxElong({} {}): SearchMaxElongation failed.'.format(name, searchText))
@@ -551,7 +540,7 @@ def CheckSaturn():
     ]
     error = 0
     for (dtext, mag, tilt) in data:
-        time = ParseDate(dtext)
+        time = astronomy.Time.Parse(dtext)
         illum = astronomy.Illumination(astronomy.Body.Saturn, time)
         print('Saturn: date={}  calc mag={:12.8f}  ring_tilt={:12.8f}'.format(dtext, illum.mag, illum.ring_tilt))
         mag_diff = abs(illum.mag - mag)
@@ -580,8 +569,8 @@ def TestMaxMag(body, filename):
             line = line.strip()
             tokenlist = line.split()
             if len(tokenlist) == 5:
-                time1 = ParseDate(tokenlist[0])
-                time2 = ParseDate(tokenlist[1])
+                time1 = astronomy.Time.Parse(tokenlist[0])
+                time2 = astronomy.Time.Parse(tokenlist[1])
                 if time1 and time2:
                     center_time = time1.AddDays(0.5*(time2.ut - time1.ut))
                     correct_mag = float(tokenlist[4])
@@ -727,7 +716,7 @@ def LunarApsis(filename):
             if len(tokenlist) != 3:
                 print('LunarApsis({} line {}): invalid data format'.format(filename, lnum))
                 return 1
-            correct_time = ParseDate(tokenlist[1])
+            correct_time = astronomy.Time.Parse(tokenlist[1])
             if not correct_time:
                 print('LunarApsis({} line {}): invalid time'.format(filename, lnum))
                 return 1
@@ -749,12 +738,6 @@ def LunarApsis(filename):
     print('LunarApsis: found {} events, max time error = {:0.3f} minutes, max distance error = {:0.3f} km.'.format(lnum, max_minutes, max_km))
     return 0
 
-
-def Test_Apsis():
-    if 0 != LunarApsis('apsides/moon.txt'):
-        return 1
-    return 0
-
 #-----------------------------------------------------------------------------------------------------------
 
 def CompareMatrices(caller, a, b, tolerance):
@@ -1051,6 +1034,78 @@ def Test_Refraction():
 
 #-----------------------------------------------------------------------------------------------------------
 
+def Test_PlanetApsis():
+    degree_threshold = 0.1
+    start_time = astronomy.Time.Make(1700, 1, 1, 0, 0, 0)
+    found_bad_planet = False
+    body = astronomy.Body.Mercury
+    while body <= astronomy.Body.Pluto:
+        count = 1
+        period = astronomy._PlanetOrbitalPeriod[body]
+        filename = os.path.join('apsides', 'apsis_{}.txt'.format(int(body)))
+        min_interval = -1.0
+        max_diff_days = 0.0
+        max_dist_ratio = 0.0
+        apsis = astronomy.SearchPlanetApsis(body, start_time)
+        with open(filename, 'rt') as infile:
+            for line in infile:
+                token = line.split()
+                if len(token) != 3:
+                    print('Test_PlanetApsis({} line {}): Invalid data format: {} tokens'.format(filename, count, len(token)))
+                    return 1
+                expected_kind = int(token[0])
+                expected_time = astronomy.Time.Parse(token[1])
+                expected_distance = float(token[2])
+                if apsis.kind != expected_kind:
+                    print('Test_PlanetApsis({} line {}): WRONG APSIS KIND: expected {}, found {}'.format(filename, count, expected_kind, apsis.kind))
+                    return 1
+                diff_days = abs(expected_time.tt - apsis.time.tt)
+                max_diff_days = max(max_diff_days, diff_days)
+                diff_degrees = (diff_days / period) * 360
+                if diff_degrees > degree_threshold:
+                    found_bad_planet = True
+                diff_dist_ratio = abs(expected_distance - apsis.dist_au) / expected_distance
+                max_dist_ratio = max(max_dist_ratio, diff_dist_ratio)
+                if diff_dist_ratio > 1.0e-4:
+                    print('Test_PlanetApsis({} line {}): distance ratio {} is too large.'.format(filename, count, diff_dist_ratio))
+                    return 1
+
+                # Calculate the next apsis.
+                prev_time = apsis.time
+                try:
+                    apsis = astronomy.NextPlanetApsis(body, apsis)
+                except astronomy.BadTimeError:
+                    if body != astronomy.Body.Pluto:
+                        raise
+                    break   # It is OK for us to go beyond Pluto calculation's time domain.
+                count += 1
+                interval = apsis.time.tt - prev_time.tt
+                if min_interval < 0.0:
+                    min_interval = max_interval = interval
+                else:
+                    min_interval = min(min_interval, interval)
+                    max_interval = max(max_interval, interval)
+            if count < 2:
+                print('Test_PlanetApsis: FAILED to find apsides for {}'.format(body))
+                return 1
+            print('Test_PlanetApsis: {:4d} apsides for {:<9s} -- intervals: min={:0.2f}, max={:0.2f}, ratio={:0.6f}; max day={:0.3f}, degrees={:0.3f}, dist ratio={:0.6f}'.format(
+                count,
+                body.name,
+                min_interval, max_interval, max_interval / min_interval,
+                max_diff_days,
+                (max_diff_days / period) * 360.0,
+                max_dist_ratio
+            ))
+        body = astronomy.Body(body + 1)
+
+    if found_bad_planet:
+        print('Test_PlanetApsis: FAIL - planet(s) exceeded angular threshold ({} degrees)'.format(degree_threshold))
+        return 1
+    print('Test_PlanetApsis: PASS')
+    return 0
+
+#-----------------------------------------------------------------------------------------------------------
+
 if __name__ == '__main__':
     if len(sys.argv) == 2:
         if sys.argv[1] == 'time':
@@ -1066,8 +1121,10 @@ if __name__ == '__main__':
             sys.exit(Test_Elongation())
         if sys.argv[1] == 'magnitude':
             sys.exit(Test_Magnitude())
-        if sys.argv[1] == 'apsis':
-            sys.exit(Test_Apsis())
+        if sys.argv[1] == 'lunar_apsis':
+            sys.exit(LunarApsis('apsides/moon.txt'))
+        if sys.argv[1] == 'planet_apsis':
+            sys.exit(Test_PlanetApsis())
         if sys.argv[1] == 'issue46':
             sys.exit(Test_Issue46())
         if sys.argv[1] == 'rotation':
diff --git a/generate/unit_test_python b/generate/unit_test_python
index 4533b957..ba6acf0f 100755
--- a/generate/unit_test_python
+++ b/generate/unit_test_python
@@ -9,7 +9,8 @@ python3 --version || Fail "Cannot print python version"
 python3 test.py rotation || Fail "Failed Python rotation tests."
 python3 test.py refraction || Fail "Failed Python refraction tests."
 python3 test.py time || Fail "Failure reported by test.py (time)"
-python3 test.py apsis || Fail "Failed Python apsis tests."
+python3 test.py lunar_apsis || Fail "Failed Python lunar apsis tests."
+python3 test.py planet_apsis || Fail "Failed Python planet apsis tests."
 python3 test.py magnitude || Fail "Failed Python magnitude tests."
 python3 test.py moon || Fail "Failure reported by test.py (moon)"
 python3 test.py seasons seasons/seasons.txt || Fail "Failed Python seasons test."
@@ -27,4 +28,5 @@ time python3 test.py astro_check > temp/py_check.txt || Fail "Failure in Python
 ./generate check temp/py_check.txt || Fail "Verification failure for Python unit test output."
 ./ctest diff temp/{py,c}_check.txt || Fail "Diff(py,js) failure."
 
+echo "$0: PASS"
 exit 0
diff --git a/pydown/py_prefix.md b/pydown/py_prefix.md
index a5120427..30d9a24f 100644
--- a/pydown/py_prefix.md
+++ b/pydown/py_prefix.md
@@ -59,6 +59,13 @@ To get started quickly, here are some [examples](../../demo/python/).
 | [SearchLunarApsis](#SearchLunarApsis) | Finds the next perigee or apogee of the Moon after a specified date. |
 | [NextLunarApsis](#NextLunarApsis) | Given an already-found apsis, finds the next perigee or apogee of the Moon. |
 
+### Planet perihelion and aphelion
+
+| Function | Description |
+| -------- | ----------- |
+| [SearchPlanetApsis](#SearchPlanetApsis) | Finds the next perihelion or aphelion of a planet after a specified date. |
+| [NextPlanetApsis](#NextPlanetApsis) | Given an already-found apsis, finds the next perihelion or aphelion of a planet. |
+
 ### Visual magnitude and elongation
 
 | Function | Description |
diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs
index 823788d3..ffd4278b 100644
--- a/source/csharp/astronomy.cs
+++ b/source/csharp/astronomy.cs
@@ -5037,8 +5037,7 @@ namespace CosineKitty
 
             for (i=0; i < npoints; ++i)
             {
-                double ut = t1.ut + (i * interval);
-                AstroTime time = new AstroTime(ut);
+                AstroTime time = t1.AddDays(i * interval);
                 double dist = NeptuneHelioDistance(time);
                 if (i == 0)
                 {
@@ -5133,7 +5132,7 @@ namespace CosineKitty
                 {
                     /* There is a change of slope polarity within the time range [t1, t2]. */
                     /* Therefore this time range contains an apsis. */
-                    /* Figure out whether it is perigee or apogee. */
+                    /* Figure out whether it is perihelion or aphelion. */
 
                     SearchContext_PlanetDistanceSlope slope_func;
                     ApsisKind kind;
diff --git a/source/python/README.md b/source/python/README.md
index 87d0a2a7..fe4f19ea 100644
--- a/source/python/README.md
+++ b/source/python/README.md
@@ -59,6 +59,13 @@ To get started quickly, here are some [examples](../../demo/python/).
 | [SearchLunarApsis](#SearchLunarApsis) | Finds the next perigee or apogee of the Moon after a specified date. |
 | [NextLunarApsis](#NextLunarApsis) | Given an already-found apsis, finds the next perigee or apogee of the Moon. |
 
+### Planet perihelion and aphelion
+
+| Function | Description |
+| -------- | ----------- |
+| [SearchPlanetApsis](#SearchPlanetApsis) | Finds the next perihelion or aphelion of a planet after a specified date. |
+| [NextPlanetApsis](#NextPlanetApsis) | Given an already-found apsis, finds the next perihelion or aphelion of a planet. |
+
 ### Visual magnitude and elongation
 
 | Function | Description |
@@ -550,6 +557,13 @@ as seen by an observer on the surface of the Earth.
 
 ---
 
+<a name="BadTimeError"></a>
+### BadTimeError
+
+Cannot calculate Pluto position for this date/time.
+
+---
+
 <a name="BadVectorError"></a>
 ### BadVectorError
 
@@ -842,6 +856,27 @@ A geocentric position vector of the center of the given body.
 
 ---
 
+<a name="HelioDistance"></a>
+### HelioDistance(body, time)
+
+**Calculates the distance between a body and the Sun at a given time.**
+
+Given a date and time, this function calculates the distance between
+the center of `body` and the center of the Sun.
+For the planets Mercury through Neptune, this function is significantly
+more efficient than calling [`HelioVector`](#HelioVector) followed by taking the length
+of the resulting vector.
+
+| Type | Parameter | Description |
+| --- | --- | --- |
+| [`Body`](#Body) | `body` | A body for which to calculate a heliocentric distance: the Sun, Moon, or any of the planets. |
+| [`Time`](#Time) | `time` | The date and time for which to calculate the heliocentric distance. |
+
+### Returns: `float`
+The heliocentric distance in AU.
+
+---
+
 <a name="HelioVector"></a>
 ### HelioVector(body, time)
 
@@ -1096,6 +1131,25 @@ the one passed in as the parameter `mq`.
 
 ---
 
+<a name="NextPlanetApsis"></a>
+### NextPlanetApsis(body, apsis)
+
+**Finds the next planetary perihelion or aphelion event in a series.**
+
+This function requires an [`Apsis`](#Apsis) value obtained from a call
+to [`SearchPlanetApsis`](#SearchPlanetApsis) or `NextPlanetApsis`.
+Given an aphelion event, this function finds the next perihelion event, and vice versa.
+See [`SearchPlanetApsis`](#SearchPlanetApsis) for more details.
+
+| Type | Parameter | Description |
+| --- | --- | --- |
+| [`Body`](#Body) | `body` | The planet for which to find the next perihelion/aphelion event. Not allowed to be `Body.Sun` or `Body.Moon`. Must match the body passed into the call that produced the `apsis` parameter. |
+| [`Apsis`](#Apsis) | `apsis` | An apsis event obtained from a call to [`SearchPlanetApsis`](#SearchPlanetApsis) or `NextPlanetApsis`. |
+
+### Returns: [`Apsis`](#Apsis)
+
+---
+
 <a name="RefractionAngle"></a>
 ### RefractionAngle(refraction, altitude)
 
@@ -1597,6 +1651,30 @@ However, the difference is minor and has little practical value.
 
 ---
 
+<a name="SearchPlanetApsis"></a>
+### SearchPlanetApsis(body, startTime)
+
+**Finds the next planet perihelion or aphelion, after a given time.**
+
+Given a date and time to start the search in `startTime`, this function finds the
+next date and time that the center of the specified planet reaches the closest or farthest point
+in its orbit with respect to the center of the Sun, whichever comes first after `startTime`.
+The closest point is called *perihelion* and the farthest point is called *aphelion*.
+The word *apsis* refers to either event.
+To iterate through consecutive alternating perihelion and aphelion events,
+call `SearchPlanetApsis` once, then use the return value to call [`NextPlanetApsis`](#NextPlanetApsis).
+After that, keep feeding the previous return value from `NextPlanetApsis`
+into another call of `NextPlanetApsis` as many times as desired.
+
+| Type | Parameter | Description |
+| --- | --- | --- |
+| [`Body`](#Body) | `body` | The planet for which to find the next perihelion/aphelion event. Not allowed to be `Body.Sun` or `Body.Moon`. |
+| [`Time`](#Time) | `startTime` | The date and time at which to start searching for the next perihelion or aphelion. |
+
+### Returns: [`Apsis`](#Apsis)
+
+---
+
 <a name="SearchRelativeLongitude"></a>
 ### SearchRelativeLongitude(body, targetRelLon, startTime)
 
diff --git a/source/python/astronomy.py b/source/python/astronomy.py
index c341a09a..64422608 100644
--- a/source/python/astronomy.py
+++ b/source/python/astronomy.py
@@ -53,6 +53,7 @@ _SECONDS_PER_DAY = 24.0 * 3600.0
 _SOLAR_DAYS_PER_SIDEREAL_DAY = 0.9972695717592592
 _MEAN_SYNODIC_MONTH = 29.530588
 _EARTH_ORBITAL_PERIOD = 365.256
+_NEPTUNE_ORBITAL_PERIOD = 60189.0
 _REFRACTION_NEAR_HORIZON = 34.0 / 60.0
 _SUN_RADIUS_AU  = 4.6505e-3
 _MOON_RADIUS_AU = 1.15717e-5
@@ -171,7 +172,7 @@ _PlanetOrbitalPeriod = [
     4332.589,
     10759.22,
     30685.4,
-    60189.0,
+    _NEPTUNE_ORBITAL_PERIOD,
     90560.0
 ]
 
@@ -200,6 +201,11 @@ class BadVectorError(Error):
     def __init__(self):
         Error.__init__(self, 'Vector is too small to have a direction.')
 
+class BadTimeError(Error):
+    """Cannot calculate Pluto position for this date/time."""
+    def __init__(self, time):
+        Error.__init__(self, 'Cannot calculate Pluto position for time: {}'.format(time))
+
 class InternalError(Error):
     """An internal error occured that should be reported as a bug.
 
@@ -2908,16 +2914,17 @@ _vsop = [
 ],
 ]
 
+def _VsopFormula(formula, t):
+    tpower = 1.0
+    coord = 0.0
+    for series in formula:
+        coord += tpower * sum(A * math.cos(B + C*t) for (A, B, C) in series)
+        tpower *= t
+    return coord
+
 def _CalcVsop(model, time):
-    spher = []
     t = time.tt / 365250.0
-    for formula in model:
-        tpower = 1.0
-        coord = 0.0
-        for series in formula:
-            coord += tpower * sum(A * math.cos(B + C*t) for (A, B, C) in series)
-            tpower *= t
-        spher.append(coord)
+    spher = [_VsopFormula(formula, t) for formula in model]
 
     # Convert spherical coordinates to ecliptic cartesian coordinates.
     r_coslat = spher[2] * math.cos(spher[1])
@@ -2931,6 +2938,12 @@ def _CalcVsop(model, time):
     vz = 0.397776982902*ey + 0.917482137087*ez
     return Vector(vx, vy, vz, time)
 
+def _VsopHelioDistance(model, time):
+    # The caller only wants to know the distance between the planet and the Sun.
+    # So we only need to calculate the radial component of the spherical coordinates.
+    # There is no need to translate coordinates.
+    return _VsopFormula(model[2], time.tt / 365250.0)
+
 def _CalcEarth(time):
     return _CalcVsop(_vsop[Body.Earth], time)
 
@@ -3130,7 +3143,7 @@ def _CalcChebyshev(model, time):
                     p1 = p2
                 pos.append(sum - coeff[0][d]/2)
             return Vector(pos[0], pos[1], pos[2], time)
-    raise Error('Cannot extrapolate Chebyshev model for given Terrestrial Time: {}'.format(time.tt))
+    raise BadTimeError(time)
 
 # END CHEBYSHEV
 #----------------------------------------------------------------------------
@@ -3354,7 +3367,7 @@ def HelioVector(body, time):
     if body == Body.Pluto:
         return _CalcChebyshev(_pluto, time)
 
-    if 0 <= body <= len(_vsop):
+    if 0 <= body < len(_vsop):
         return _CalcVsop(_vsop[body], time)
 
     if body == Body.Sun:
@@ -3368,6 +3381,37 @@ def HelioVector(body, time):
     raise InvalidBodyError()
 
 
+def HelioDistance(body, time):
+    """Calculates the distance between a body and the Sun at a given time.
+
+    Given a date and time, this function calculates the distance between
+    the center of `body` and the center of the Sun.
+    For the planets Mercury through Neptune, this function is significantly
+    more efficient than calling #HelioVector followed by taking the length
+    of the resulting vector.
+
+    Parameters
+    ----------
+    body : Body
+        A body for which to calculate a heliocentric distance:
+        the Sun, Moon, or any of the planets.
+    time : Time
+        The date and time for which to calculate the heliocentric distance.
+
+    Returns
+    -------
+    float
+        The heliocentric distance in AU.
+    """
+    if body == Body.Sun:
+        return 0.0
+
+    if 0 <= body < len(_vsop):
+        return _VsopHelioDistance(_vsop[body], time)
+
+    return HelioVector(body, time).Length()
+
+
 def GeoVector(body, time, aberration):
     """Calculates geocentric Cartesian coordinates of a body in the J2000 equatorial system.
 
@@ -5031,7 +5075,7 @@ def Seasons(year):
 def _MoonDistance(time):
     return _CalcMoon(time).distance_au
 
-def _distance_slope(direction, time):
+def _moon_distance_slope(direction, time):
     dt = 0.001
     t1 = time.AddDays(-dt/2.0)
     t2 = time.AddDays(+dt/2.0)
@@ -5118,11 +5162,11 @@ def SearchLunarApsis(startTime):
     """
     increment = 5.0     # number of days to skip on each iteration
     t1 = startTime
-    m1 = _distance_slope(+1, t1)
+    m1 = _moon_distance_slope(+1, t1)
     iter = 0
     while iter * increment < 2.0 * _MEAN_SYNODIC_MONTH:
         t2 = t1.AddDays(increment)
-        m2 = _distance_slope(+1, t2)
+        m2 = _moon_distance_slope(+1, t2)
         if m1 * m2 <= 0.0:
             # There is a change of slope polarity within the time range [t1, t2].
             # Therefore this time range contains an apsis.
@@ -5130,12 +5174,12 @@ def SearchLunarApsis(startTime):
             if m1 < 0.0 or m2 > 0.0:
                 # We found a minimum-distance event: perigee.
                 # Search the time range for the time when the slope goes from negative to positive.
-                apsis_time = Search(_distance_slope, +1, t1, t2, 1.0)
+                apsis_time = Search(_moon_distance_slope, +1, t1, t2, 1.0)
                 kind = ApsisKind.Pericenter
             elif m1 > 0.0 or m2 < 0.0:
                 # We found a maximum-distance event: apogee.
                 # Search the time range for the time when the slope goes from positive to negative.
-                apsis_time = Search(_distance_slope, -1, t1, t2, 1.0)
+                apsis_time = Search(_moon_distance_slope, -1, t1, t2, 1.0)
                 kind = ApsisKind.Apocenter
             else:
                 # This should never happen. It should not be possible for both slopes to be zero.
@@ -5182,10 +5226,189 @@ def NextLunarApsis(apsis):
     if apsis.kind not in [ApsisKind.Apocenter, ApsisKind.Pericenter]:
         raise Error('Parameter "apsis" contains an invalid "kind" value.')
     if next.kind + apsis.kind != 1:
-        raise InternalError()
+        raise InternalError()   # should have found opposite apsis kind
     return next
 
 
+def _planet_distance_slope(context, time):
+    (direction, body) = context
+    dt = 0.001
+    t1 = time.AddDays(-dt/2.0)
+    t2 = time.AddDays(+dt/2.0)
+    dist1 = HelioDistance(body, t1)
+    dist2 = HelioDistance(body, t2)
+    return direction * (dist2 - dist1) / dt
+
+
+def SearchPlanetApsis(body, startTime):
+    """Finds the next planet perihelion or aphelion, after a given time.
+
+    Given a date and time to start the search in `startTime`, this function finds the
+    next date and time that the center of the specified planet reaches the closest or farthest point
+    in its orbit with respect to the center of the Sun, whichever comes first after `startTime`.
+
+    The closest point is called *perihelion* and the farthest point is called *aphelion*.
+    The word *apsis* refers to either event.
+
+    To iterate through consecutive alternating perihelion and aphelion events,
+    call `SearchPlanetApsis` once, then use the return value to call #NextPlanetApsis.
+    After that, keep feeding the previous return value from `NextPlanetApsis`
+    into another call of `NextPlanetApsis` as many times as desired.
+
+    Parameters
+    ----------
+    body : Body
+        The planet for which to find the next perihelion/aphelion event.
+        Not allowed to be `Body.Sun` or `Body.Moon`.
+    startTime : Time
+        The date and time at which to start searching for the next perihelion or aphelion.
+
+    Returns
+    -------
+    Apsis
+    """
+    if body == Body.Neptune:
+        return _SearchNeptuneApsis(startTime)
+    positive_slope = (+1.0, body)
+    negative_slope = (-1.0, body)
+    orbit_period_days = _PlanetOrbitalPeriod[body]
+    increment = orbit_period_days / 6.0
+    t1 = startTime
+    m1 = _planet_distance_slope(positive_slope, t1)
+    iter = 0
+    while iter * increment < 2 * orbit_period_days:
+        t2 = t1.AddDays(increment)
+        m2 = _planet_distance_slope(positive_slope, t2)
+        if m1 * m2 <= 0.0:
+            # There is a change of slope polarity within the time range [t1, t2].
+            # Therefore this time range contains an apsis.
+            # Figure out whether it is perihelion or aphelion.
+            if m1 < 0.0 or m2 > 0.0:
+                # We found a minimum-distance event: perihelion.
+                # Search the time range for the time when the slope goes from negative to positive.
+                context = positive_slope
+                kind = ApsisKind.Pericenter
+            elif m1 > 0.0 or m2 < 0.0:
+                # We found a maximum-distance event: aphelion.
+                # Search the time range for the time when the slope goes from positive to negative.
+                context = negative_slope
+                kind = ApsisKind.Apocenter
+            else:
+                raise InternalError()   # at least one of the planet distance slopes should have been nonzero
+            search = Search(_planet_distance_slope, context, t1, t2, 1.0)
+            if search is None:
+                raise InternalError()   # failed to find where planet distance slope passed through zero
+
+            dist = HelioDistance(body, search)
+            return Apsis(search, kind, dist)
+        t1 = t2
+        m1 = m2
+        iter += 1
+    raise InternalError()   # should have found planet apsis within 2 planet orbits
+
+
+def NextPlanetApsis(body, apsis):
+    """Finds the next planetary perihelion or aphelion event in a series.
+
+    This function requires an #Apsis value obtained from a call
+    to #SearchPlanetApsis or `NextPlanetApsis`.
+    Given an aphelion event, this function finds the next perihelion event, and vice versa.
+    See #SearchPlanetApsis for more details.
+
+    Parameters
+    ----------
+    body : Body
+        The planet for which to find the next perihelion/aphelion event.
+        Not allowed to be `Body.Sun` or `Body.Moon`.
+        Must match the body passed into the call that produced the `apsis` parameter.
+    apsis : Apsis
+        An apsis event obtained from a call to #SearchPlanetApsis or `NextPlanetApsis`.
+
+    Returns
+    -------
+    Apsis
+    """
+    if apsis.kind not in [ApsisKind.Apocenter, ApsisKind.Pericenter]:
+        raise Error('Parameter "apsis" contains an invalid "kind" value.')
+    # Skip 1/4 of an orbit before starting search again.
+    skip = 0.25 * _PlanetOrbitalPeriod[body]
+    time = apsis.time.AddDays(skip)
+    next = SearchPlanetApsis(body, time)
+    # Verify that we found the opposite apsis from the previous one.
+    if next.kind + apsis.kind != 1:
+        raise InternalError()   # should have found opposite planetary apsis type
+    return next
+
+def _NeptuneHelioDistance(time):
+    return _VsopHelioDistance(_vsop[Body.Neptune], time)
+
+def _NeptuneExtreme(kind, start_time, dayspan):
+    direction = +1.0 if (kind == ApsisKind.Apocenter) else -1.0
+    npoints = 10
+    while True:
+        interval = dayspan / (npoints - 1)
+        # Iterate until uncertainty is less than one minute.
+        if interval < 1/1440:
+            apsis_time = start_time.AddDays(interval/2)
+            dist_au = _NeptuneHelioDistance(apsis_time)
+            return Apsis(apsis_time, kind, dist_au)
+        for i in range(npoints):
+            time = start_time.AddDays(i * interval)
+            dist = direction * _NeptuneHelioDistance(time)
+            if i==0 or dist > best_dist:
+                best_i = i
+                best_dist = dist
+        # Narrow in on the extreme point.
+        start_time = start_time.AddDays((best_i - 1) * interval)
+        dayspan = 2 * interval
+
+def _SearchNeptuneApsis(startTime):
+    # Neptune is a special case for two reasons:
+    # 1. Its orbit is nearly circular (low orbital eccentricity).
+    # 2. It is so distant from the Sun that the orbital period is very long.
+    # Put together, this causes wobbling of the Sun around the Solar System Barycenter (SSB)
+    # to be so significant that there are 3 local minima in the distance-vs-time curve
+    # near each apsis. Therefore, unlike for other planets, we can't use an optimized
+    # algorithm for finding dr/dt = 0.
+    # Instead, we use a dumb, brute-force algorithm of sampling and finding min/max
+    # heliocentric distance.
+    #
+    # Rewind approximately 30 degrees in the orbit,
+    # then search forward for 270 degrees.
+    # This is a very cautious way to prevent missing an apsis.
+    # Typically we will find two apsides, and we pick whichever
+    # apsis is ealier, but after startTime.
+    # Sample points around this orbital arc and find when the distance
+    # is greatest and smallest.
+    t1 = startTime.AddDays(_NEPTUNE_ORBITAL_PERIOD * ( -30 / 360))
+    t2 = startTime.AddDays(_NEPTUNE_ORBITAL_PERIOD * (+270 / 360))
+    t_min = t_max = t1
+    npoints = 1000
+    interval = (t2.ut - t1.ut) / (npoints - 1)
+
+    for i in range(npoints):
+        time = t1.AddDays(i * interval)
+        dist = _NeptuneHelioDistance(time)
+        if i == 0:
+            min_dist = max_dist = dist
+        else:
+            if dist > max_dist:
+                max_dist = dist
+                t_max = time
+            if dist < min_dist:
+                min_dist = dist
+                t_min = time
+
+    perihelion = _NeptuneExtreme(ApsisKind.Pericenter, t_min.AddDays(-2*interval), 4*interval)
+    aphelion = _NeptuneExtreme(ApsisKind.Apocenter, t_max.AddDays(-2*interval), 4*interval)
+    if perihelion.time.tt >= startTime.tt:
+        if startTime.tt <= aphelion.time.tt < perihelion.time.tt:
+            return aphelion
+        return perihelion
+    if aphelion.time.tt >= startTime.tt:
+        return aphelion
+    raise InternalError()   # failed to find Neptune apsis
+
 def VectorFromSphere(sphere, time):
     """Converts spherical coordinates to Cartesian coordinates.