Change main computation routine to a generator

This is kind of magic.  Because it is now a generator, you have the
option of asking for position, then only proceeding to the velocity if
you need it, without having to repeat the work of getting to that point
again.

bin/benchmark.py

@@ -9,9 +9,9 @@ from jplephem.spk import SPK
 def main():
     for size in 10, 1000, 100000:
         jd = np.linspace(2414992.5, 2471184.50, size)
-        spk = SPK.open('de421.bsp')
+        kernel = SPK.open('de421.bsp')
         ephem = Ephemeris(de421)
-        mars = spk.targets[4]
+        mars = kernel[0,4]
 
         print(size)
         print('-- old code (2 successive runs):')

jplephem/__init__.py

@@ -104,7 +104,7 @@ of a 3-vector.  For an ephemeris that does not, you can ask for the
 Chebyshev polynomial to be differentiated to produce a velocity, which
 is delivered as a second return value:
 
->>> position, velocity = kernel[0,4].compute(2457061.5, differentiate=True)
+>>> position, velocity = kernel[0,4].compute_and_differentiate(2457061.5)
 >>> print(position)
 [  2.05700211e+08   4.25141646e+07   1.39379183e+07]
 >>> print(velocity)
@@ -114,8 +114,8 @@ is delivered as a second return value:
 Details of the API
 ------------------
 
-Here are a few details for people wanting to go beyond the high-level
-API provided above and read through the code to learn about the details.
+Here are a few details for people ready to go beyond the high-level API
+provided above and read through the code to learn more.
 
 * Instead of reading an entire ephemeris into memory, ``jplephem``
   memory-maps the underlying file so that the operating system can

jplephem/jpltest.py

@@ -59,21 +59,20 @@ def run_testpo(spk, testpo_file):
 def _position(spk, jed, target):
     """Compute position given a JPL test file target integer identifier."""
 
-
     if target == 3:
-        p1, v1 = spk[0,3].compute(jed, differentiate=True)
-        p2, v2 = spk[3,399].compute(jed, differentiate=True)
+        p1, v1 = spk[0,3].compute_and_differentiate(jed)
+        p2, v2 = spk[3,399].compute_and_differentiate(jed)
         p = p1 + p2
         v = v1 + v2
     elif target == 10:
-        p1, v1 = spk[0,3].compute(jed, differentiate=True)
-        p2, v2 = spk[3,301].compute(jed, differentiate=True)
+        p1, v1 = spk[0,3].compute_and_differentiate(jed)
+        p2, v2 = spk[3,301].compute_and_differentiate(jed)
         p = p1 + p2
         v = v1 + v2
     elif target == 12:
         return np.zeros((6, 1))  # solar system barycenter is the origin
     else:
-        p, v = spk[0,target].compute(jed, differentiate=True)
+        p, v = spk[0,target].compute_and_differentiate(jed)
 
     return np.concatenate((p, v))
 

jplephem/spk.py

@@ -103,6 +103,15 @@ class Segment(object):
                      .format(self, self.source.decode('ascii')))
         return text
 
+    def compute(self, tdb, tdb2=0.0):
+        """Compute the component values for the time `tdb` plus `tdb2`."""
+        for position in self.generate(tdb, tdb2):
+            return position
+
+    def compute_and_differentiate(self, tdb, tdb2=0.0):
+        """Compute components and differentials for time `tdb` plus `tdb2`."""
+        return list(self.generate(tdb, tdb2))
+
     def _load(self):
         """Map the coefficients into memory using a NumPy array."""
 
@@ -125,15 +134,15 @@ class Segment(object):
         coefficients = rollaxis(coefficients, 1)
         return initial_epoch, interval_length, coefficients
 
-    def compute(self, tdb, tdb2=0.0, differentiate=False):
-        """Compute the component values for the time `tdb` plus `tdb2`.
+    def generate(self, tdb, tdb2):
+        """Generate components and differentials for time `tdb` plus `tdb2`.
 
-        If `differentiate` is false, then an array of components is
-        returned.
-
-        If `differentiate` is true, then a tuple is returned whose first
-        element is the array of components and whose second element is
-        an array of rates at which the components are changing.
+        Most uses will simply want to call the `compute()` method or the
+        `compute_differentials()` method, for convenience.  But in those
+        cases (see Skyfield) where you want to compute a position and
+        examine it before deciding whether to proceed with the velocity,
+        but without losing all of the work that it took to get to that
+        point, this generator lets you get them as two separate steps.
 
         """
         scalar = not getattr(tdb, 'shape', 0) and not getattr(tdb2, 'shape', 0)
@@ -175,8 +184,8 @@ class Segment(object):
         components = (T.T * coefficients).sum(axis=2)
         if scalar:
             components = components[:,0]
-        if not differentiate:
-            return components
+
+        yield components
 
         # Chebyshev differentiation.
 
@@ -192,7 +201,8 @@ class Segment(object):
         rates = (dT.T * coefficients).sum(axis=2)
         if scalar:
             rates = rates[:,0]
-        return components, rates
+
+        yield rates
 
 
 def titlecase(name):

jplephem/test.py

@@ -153,7 +153,7 @@ class SPKTests(_CommonTests, TestCase):
 
     def position_and_velocity(self, name, tdb, tdb2=0.0):
         segment = self.spk[0, target_names[name]]
-        return segment.compute(tdb, tdb2, differentiate=True)
+        return segment.compute_and_differentiate(tdb, tdb2)
 
     def test_str(self):
         str(self.spk)  # just to confirm it does not raise an exception