Calculate barycentric state of Pluto.

The BaryState function did not support Pluto before.
Refactored the code so that the internal CalcPluto function
returns both the position and velocity, and its caller
can select from heliocentric or barycentric coordinates.
HelioVector asks for heliocentric coordinates and keeps
only the position vector. BaryState asks for barycentric
coordinates and returns both position and velocity.

I added test data for Pluto generated by JPL Horizons.
It turns out the Pluto system barycenter is the best fit
for TOP2013, presumably because Charon causes Pluto to
wobble quite a bit.

I also generated JPL Horizons test data for the Moon
and the Earth/Moon barycenter, anticipating that I will
support calculating their barycentric state vectors soon.

I had to increase the enforced size limit for minified
JavaScript from 100000 bytes to 120000 bytes.
I guess this is like raising the "debt ceiling".

Fixed a bug in Python unit tests: if "-v" verbose option
was specified, it was printing a summary line for every
single line of input, instead of a single summary after
processing the whole file, as was intended. This is one
of those Python whitespace indentation bugs!

source/python/README.md

@@ -1031,7 +1031,7 @@ The vectors are expressed in equatorial J2000 coordinates (EQJ).
 
 | Type | Parameter | Description |
 | --- | --- | --- |
-| [`Body`](#Body) | `body` | The celestial body whose barycentric state vector is to be calculated. Supported values are `Body.Sun`, `Body.SSB`, and all planets except Pluto: `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, `Body.Jupiter`, `Body.Saturn`, `Body.Uranus`, `Body.Neptune`. |
+| [`Body`](#Body) | `body` | The celestial body whose barycentric state vector is to be calculated. Supported values are `Body.Sun`, `Body.SSB`, and all planets: `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, `Body.Jupiter`, `Body.Saturn`, `Body.Uranus`, `Body.Neptune`, `Body.Pluto`. |
 | [`Time`](#Time) | `time` | The date and time for which to calculate position and velocity. |
 
 ### Returns: [`StateVector`](#StateVector)

source/python/astronomy.py

@@ -3472,6 +3472,13 @@ def _UpdatePosition(dt, r, v, a):
         r.z + dt*(v.z + dt*a.z/2.0)
     )
 
+def _UpdateVelocity(dt, v, a):
+    return _TerseVector(
+        v.x + dt*a.x,
+        v.y + dt*a.y,
+        v.z + dt*a.z
+    )
+
 
 def _GravSim(tt2, calc1):
     dt = tt2 - calc1.tt
@@ -3567,7 +3574,8 @@ def _CalcPlutoOneWay(entry, target_tt, dt):
     return sim
 
 
-def _CalcPluto(time):
+def _CalcPluto(time, helio):
+    bary = None
     seg = _GetSegment(_pluto_cache, time.tt)
     if seg is None:
         # The target time is outside the year range 0000..4000.
@@ -3578,6 +3586,7 @@ def _CalcPluto(time):
         else:
             sim = _CalcPlutoOneWay(_PlutoStateTable[_PLUTO_NUM_STATES-1], time.tt, +_PLUTO_DT)
         r = sim.grav.r
+        v = sim.grav.v
         bary = sim.bary
     else:
         left = _ClampIndex((time.tt - seg[0].tt) / _PLUTO_DT, _PLUTO_NSTEPS-1)
@@ -3589,15 +3598,25 @@ def _CalcPluto(time):
 
         # Use Newtonian mechanics to extrapolate away from t1 in the positive time direction.
         ra = _UpdatePosition(time.tt - s1.tt, s1.r, s1.v, acc)
+        va = _UpdateVelocity(time.tt - s1.tt, s1.v, acc)
 
         # Use Newtonian mechanics to extrapolate away from t2 in the negative time direction.
         rb = _UpdatePosition(time.tt - s2.tt, s2.r, s2.v, acc)
+        vb = _UpdateVelocity(time.tt - s2.tt, s2.v, acc)
 
         # Use fade in/out idea to blend the two position estimates.
         ramp = (time.tt - s1.tt)/_PLUTO_DT
         r = ra*(1 - ramp) + rb*ramp
-        bary = _major_bodies_t(time.tt)
-    return (r - bary.Sun.r).ToAstroVector(time)
+        v = va*(1 - ramp) + vb*ramp
+
+    if helio:
+        # Convert barycentric vectors to heliocentric vectors.
+        if bary is None:
+            bary = _major_bodies_t(time.tt)
+        r -= bary.Sun.r
+        v -= bary.Sun.v
+
+    return StateVector(r.x, r.y, r.z, v.x, v.y, v.z, time)
 
 
 # END Pluto Integrator
@@ -4087,7 +4106,8 @@ def HelioVector(body, time):
         at the given time.
     """
     if body == Body.Pluto:
-        return _CalcPluto(time)
+        planet = _CalcPluto(time, True)
+        return Vector(planet.x, planet.y, planet.z, time)
 
     if 0 <= body.value < len(_vsop):
         return _CalcVsop(_vsop[body.value], time)
@@ -4234,9 +4254,9 @@ def BaryState(body, time):
     ----------
     body : Body
         The celestial body whose barycentric state vector is to be calculated.
-        Supported values are `Body.Sun`, `Body.SSB`, and all planets except Pluto:
+        Supported values are `Body.Sun`, `Body.SSB`, and all planets:
         `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, `Body.Jupiter`,
-        `Body.Saturn`, `Body.Uranus`, `Body.Neptune`.
+        `Body.Saturn`, `Body.Uranus`, `Body.Neptune`, `Body.Pluto`.
     time : Time
         The date and time for which to calculate position and velocity.
 
@@ -4249,6 +4269,9 @@ def BaryState(body, time):
     if body == Body.SSB:
         return StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time)
 
+    if body == Body.Pluto:
+        return _CalcPluto(time, False)
+
     # Find the barycentric positions and velocities for the 5 major bodies.
     bary = _major_bodies_t(time.tt)
 
@@ -4286,7 +4309,7 @@ def BaryState(body, time):
             time
         )
 
-    # FIXFIXFIX: later, we can add support for Pluto, Moon, EMB, etc.
+    # FIXFIXFIX: later, we can add support for Moon, EMB, etc.
     raise InvalidBodyError()