Description of Find_Orb's computation and caching of planetary positions:

   Early on in the history of Find_Orb,  I realized that the program spent a
_lot_ of time computing planetary positions.  I first dealt with this by
optimizing the JPL DE ephemeris code (jpleph.cpp).  That did help a fair bit,
but planetary positions still accounted for most of the CPU cycles.

   My next step was to try to re-use planetary positions by caching them:
compute 'em once,  then store them and recycle the computation the next time
the position was needed for the same time and object.  When you do a "full
step",  for example,  the program will usually integrate over a given time
range seven times.  With caching,  the first such integration would be pretty
slow,  with every position calculated from DE.  The subsequent integrations
of slightly "tweaked" particles (see 'full.txt') would almost entirely use
the same planetary positions,  with maybe a few new ones required due to
the "tweaked" object taking a slightly different path.

   I tried a series of techniques for doing this.  At first,  I stored the
planetary positions using a sort of balanced binary tree scheme,  which was
slow and overly complicated.  Then I used a hash table.  This was much
simpler and faster,  but too much time was still going into extracting the
positions from the hash table.  I think that was because the hash table
was usually tens or hundreds of MBytes in extent;  there was a lot of
thrashing around in memory going on.

   I did not really need anything very sophisticated.  The functions
required were "store this planet's position at this time" and "retrieve
that planet's position at that time".  No deletions,  in-order
traversals,  etc.  A good bit of dead space would be acceptable.

   What I ended up with was a tree with one root node and a theoretically
infinite number of leaves.  Planetary positions are stored in "leaf nodes"
of maybe a few hundred or a few thousand positions each. Each leaf node
contains the positions for a particular date range.  For example, by the
time we have enough positions cached to require three leaf nodes,  the data
might be separated as follows:

First node has a minimum JD of -infinity;
Second node has a minimum JD of 2448345.349;
Third node has a minimum JD of 2451545.018 (and implicitly,  a maximum JD
   of +infinity)

   In this case,  all positions with JDs less than 2448345.349 would be
cached in the first node.  Anything greater than that,  but before
2451545.018, would go into the second node.  Anything past that,  all the
way to +inf, would go into the third node.

   As new positions are computed,  they go into these nodes.  (We start out
with a single,  empty node with minimum JD -inf,  and implicit max JD of
+inf.)  If a node gets to being 80% full,  it's split in half (I will
explain shortly why we don't wait until it's completely full).  To split
it,  we sort the entries by JD;  we look at the JD of the middle entry of
the result, allocate a new node,  and split the entries between the old
and new nodes.

   (These bits are all pretty much standard B+ tree procedure. However,
having only one "root" index node with "leaves" directly below it really
simplifies thing tremendously.)

   It helps greatly that access is almost uniformly sequential;  i.e.,  we're
either integrating forward in time or backward in time,  accessing lots of
positions along the way.  So about 99% of the time,  we can just look within
whichever node was last accessed.  (Even when we don't,  binary searching
within the root node is fast.)

   A further wrinkle,  which I cannot believe is original to me,  but which
I've not seen discussed elsewhere,  lies in how items are stored within a
node.  Obviously,  once you have several hundred or thousand positions stored
within a node,  you don't want to have to say (to use the above example):
"We're trying to get a position for Jupiter for JD 2450000.0;  that's between
JD 2448345.349 and less than 2451545.018,  so we know we should look in the
second node,  but there are a few thousand items currently in that node;
let's look at every single one of those few thousand,  trying to find
the one we want."

   To get around this,  the items within a node are stored using a hash
table scheme.  JD=2450000.5 and planet_number = 5 are hashed,  and we look
at that point within the node,  doing a quadratic search.  (This is all
bog-standard hash table usage.)  In practice,  the nodes are usually about
half full,  so we don't have to do much searching.

   As hash tables fill,  collisions are more frequent,  which is why we
split nodes when they're 80% full.

   Note that deletions and in-order traversals would be hard to do.  But
we aren't doing them.  Insertion and finding are really,  really fast,  and
that's what we're doing here.  The method is very specific to the job
at hand.

   Yet another small wrinkle:  when a node is to be split,  we first check
its neighbors (or neighbor singular,  if it's the first or last node).  If
one of those neighbors is less than half full,  we don't actually split;  we
just re-distribute the "full" node so that both nodes end up with an equal
number of entries.  (This is referred to in the code as "spilling over" into
an adacent node.)  If,  for example,  the neighboring node is 30% full and
the node we're about to split is 80% full,  we re-distribute so that both
nodes are 55% full.  Doing this results in slightly fuller nodes,  and
therefore in fewer nodes being needed for a given number of entries.  (This
idea is borrowed from B-star trees,  though considerably simplified.  It
again helps greatly than all we do is add and search for nodes;  deletions
aren't an issue,  nor do we need to worry about in-order traversals.)
