---

worldline.C

Go to the documentation of this file.

       //=======================================================//    _\|/_
      //  __  _____           ___                    ___       //      /|\ ~
     //  /      |      ^     |   \  |         ^     |   \     //          _\|/_
    //   \__    |     / \    |___/  |        / \    |___/    //            /|\ ~
   //       \   |    /___\   |  \   |       /___\   |   \   // _\|/_
  //     ___/   |   /     \  |   \  |____  /     \  |___/  //   /|\ ~
 //                                                       //            _\|/_
//=======================================================//              /|\ ~

//
//----------------------------------------------------------------------
//
// Worldbundle member functions:
//
//      worldbundle::worldbundle(tdyn *b)
//      void worldbundle::print()
//      void worldbundle::dump()
//      int worldbundle::find_index(real id)
//      int worldbundle::find_index(char *name)
//      int worldbundle::find_index(_pdyn_ *b)
//      worldline *worldbundle::find_worldline(real id)
//      worldline *worldbundle::find_worldline(char *name)
//      worldline *worldbundle::find_worldline(_pdyn_ *b)
//      void worldbundle::attach(tdyn *bn, int verbose)
//      void worldbundle::print_worldline(char *name, real dt)
//
// Worldline member function:
//
//      void worldline::dump(int offset)
//
// Segment member function:
//
//      void segment::print(char *label)
//
// Other externally visible functions:
//
//      unique_id_t unique_id(char *name)
//      unique_id_t unique_id(node *b)
//      void print_id(void *id, char *label)
//      worldbundle *read_bundle(istream &s, int verbose)
//      tdyn *find_event(worldline *w, tdyn *bn, real t)
//      void print_event(worldline *w, tdyn *bn, real t)
//      vector interpolate_pos(tdyn *p, real t, tdyn *bn)
//      vector interpolate_pos(tdyn *p, real t, vector& pos, bool inc, tdyn *bn)
//      void set_interpolated_pos(tdyn *p, real t, vector& pos, tdyn *bn)
//      void set_interpolated_pos(tdyn *p, real t, pdyn *curr, tdyn *bn);
//      void inc_interpolated_pos(tdyn *p, real t, vector& pos, tdyn *bn)
//      vector get_pos(worldline *w, tdyn *b, tdyn *bn, real t)
//      pdyn *create_interpolated_tree(worldbundle *wb, real t)
//
//      real physical_mass_scale()
//
//----------------------------------------------------------------------

#include "worldline.h"
#include "inline.h"
#include <ctype.h>

#define NEW 0

#ifndef TOOLBOX

// Functions to define unique identification for particles.

#ifndef NEW_UNIQUE_ID_T

    // For use in old code:

    //#define ID_FAC    (1.0/8192)              // power of 2...
    //#define ID_FAC    (1.0/131072)            // power of 2...

#   define ID_FAC 0.000001                      // easier to read

    // To keep ids unique, need 1/FAC to exceed the maximum particle index.
    // Should possibly allow this to be varied at run time...

#else

//  For use in case NEW_UNIQUE_ID_T:

#   define ID_SHIFT 20

#endif

// Overloaded function...

unique_id_t unique_id(int index)
{
#ifndef NEW_UNIQUE_ID_T
    return (unique_id_t)(1 + ID_FAC*index);
#else
    return (unique_id_t)((1<<ID_SHIFT) + index);
#endif
}

unique_id_t unique_id(char *name)
{
    // Associate a unique identifier with the specified name.
    // This ID will be used to track the particle throughout.
    // Old style: a leaf has 1 < id < 2, a binary CM has 2 < id < 3, etc.
    // New style: a leaf has id>>20 = 1, binary CM has id>>20 = 2, etc.
    //
    // NOTE (Steve, 12/01): actually only necessary that the id be unique
    // at any given instant.  Not necessary that the id is unique for all
    // times (successive segments of a worldline may then refer to
    // different nodes, but that is OK...).

    // First see if name can be interpreted as an integer.

    int index = atoi(name);
    if (index > 0) return unique_id(index);

    // This is presumably a center of mass node.  Turn its
    // name into an identifying number.

    // Parse the name.

    unique_id_t id = 0;
    int n = 0;

#if 0

    // Old version:
    //
    // Look for commas and parentheses separating integers.
    // Combine ((a,b),c) --> a + ID_FAC*(b + ID_FAC*c), etc.
    // Factor ID_FAC is ~arbitrary (but see above).
    // Really need a better way of generating unique IDs...

    char tname[1024];
    strcpy(tname, name);        // work with a local copy of the name
    int l = strlen(tname);

    real fac = 1;
    int num = 0;

    for (int i = 0; i < l; i++) {
        if (tname[i] >= '0' && tname[i] <= '9') {
            if (num == 0) {
                fac *= ID_FAC;
                num = i+1;
            }
        } else {
            if (num > 0) {
                tname[i] = '\0';
                id += fac*atoi(tname+num-1);
                n++;
                num = 0;
            }
        }
    }

#else

    // Improved code:  no copy and uses strtol (Steve, 12/01).
    //
    // New (12/01) version encodes only the *first* number found
    // in name, relying on the value of n to keep id unique.
    //
    // New new version (12/01) uses an int to store the data, rather than
    // a real.  Versions distinguished by value of NEW_UNIQUE_ID_T.

    char *c = name, *end;
    bool num = false;

    while (1) {

        // Find the next digit in name.

        while (*c && !isdigit(*c)) c++;

        if (*c) {

            // New number starts at c.  Read and count it.  With the new code,
            // it may be quicker just to count commas rather than use strol
            // for numbers that won't be used.

            long int i = strtol(c, &end, 10);
            n++;

            if (n == 1) {

#ifndef NEW_UNIQUE_ID_T
                id = (unique_id_t)(ID_FAC*i);   // new (12/01)
#else
                id = (unique_id_t)(i);          // new new (12/01)
#endif
            }

            c = end;            // starting point for next search

        } else
            break;
    }

#endif

    // Portion of id before the decimal point or shifted by 20 indicates
    // the number of numbers found in the name.

#ifndef NEW_UNIQUE_ID_T
    id += n;
#else
    id += (n<<ID_SHIFT);                                // new new (12/01)
#endif

    return id;
}

unique_id_t unique_id(node *b)
{
    // Return a unique non-negative real number identifying this node.
    // Use the index if available (assume unique!), otherwise, construct
    // an integer from the node name.  If the node has no name or index,
    // ignore it (id = -1).

    // The root node should be called "root", and will therefore
    // be assigned an ID of 0.

    if (b->get_index() > 0)
        return unique_id(b->get_index());

    else if (b->get_name())
        return unique_id(b->format_label());

    else
        return -1;
}

void print_id(void *id, char *s = NULL, int offset = 0)
{
    real *r = (real *)id;
    unsigned long long *u = (unsigned long long *)id;

    real rr = *r;
    unsigned long long uu = *u;

    if (offset > 0) PRI(offset);
    if (s)
        cerr << s << ": ";
    else
        cerr << "    ";

    fprintf(stderr, "%.16lf = %qu = %qx\n", rr, uu, uu);
}

void worldline::dump(int offset,        // default = 0
                     real t1,           // default = 0
                     real t2)           // default = VERY_LARGE_NUMBER
{
    // Print a worldline, from time t1 to time t2.

    PRI(offset+4); PRC(start_esc_flag); PRC(end_esc_flag); PRL(t_esc);

    segment *s = get_first_segment();
    int is = 0;
    while (s) {
        real t_start = s->get_t_start();
        real t_end = s->get_t_end();
        if (t_start <= t2 && t_end >= t1) { 
            PRI(offset+4); cerr << "segment " << is << " (" << s << "): ";
            PRC(t_start); PRL(t_end);
            print_id(&id, "id", offset+8);
        }
        tdyn *b = s->get_first_event();
        int ie = 0;
        while (b) {
            real t = b->get_time();
            if (t >= t1 && t <= t2) {
                PRI(offset+8);
                cerr << "event " << ie << ": name = "
                     << b->format_label() << ", ";
                PRL(t);
            }
            b = b->get_next();
            ie++;
        }
        if (t_start <= t2 && t_end >= t1) { 
            if (ie == 1) {
                PRI(offset+8); cerr << "*** single event ***" << endl;
            }
        }
        s = s->get_next();
        is++;
    }
}

local void check_print_worldline_header(bool& print, int i)
{
    cerr << "worldline " << i << ":" << endl;
    print = false;
}

void worldline::check(int i)    // default = -1
{
    // Check a worldline, start to finish.

    bool print = true;

    segment *s = get_first_segment();
    int is = 0;
    while (s) {
        tdyn *b0 = s->get_first_event();
        tdyn *b = b0;
        int ie = 0;
        while (b) {
            real t = b->get_time();
            if (b != b0) {
                tdyn *p = b->get_prev();
                if (!p) {
                    if (print) check_print_worldline_header(print, i);
                    cerr << "    segment " << is << ", event "
                         << ie << " has prev = NULL" << endl;
                } else if (p->get_next() != b) {
                    if (print) check_print_worldline_header(print, i);
                    cerr << "    segment " << is << ", event "
                         << ie << " has prev->next != this" << endl;
                }
            }
            b = b->get_next();
            ie++;
        }
        if (ie == 1) {
            if (print) check_print_worldline_header(print, i);
            cerr << "    *** single event in segment " << is << " ***" << endl;
        }

        s = s->get_next();
        is++;
    }
}

void segment::print(char *label)
{
    if (label == NULL) label = "    ";
    cerr << label << "name " << first_event->format_label()
         << ", id = " << id
         << ", t_start = " << t_start
         << ", t_end = " << t_end
         << endl;
}

//======================================================================

local int wlcompare(const void *a, const void *b)       // a and b are actually
                                                        // of type worldlineptr*
{
    worldlineptr wa = *((worldlineptr*)a);              // ugly!
    worldlineptr wb = *((worldlineptr*)b);

    int result = 0;
//    cerr << "comparing " << wa->get_id() << " and " << wb->get_id() << endl;

    if (wa->get_id() < wb->get_id())
        result = -1;
    else if (wa->get_id() > wb->get_id())
        result = 1;

//    PRL(result);
    return result;
}

worldbundle::worldbundle(tdyn *b)
{
    // Create a new worldbundle from the nodes below b.

    nw_max = 0;
    for_all_nodes(tdyn, b, bb) nw_max++;

    nw_max *= 2;                                // conservative

    bundle = new worldlineptr[nw_max];
    nw = 0;

    t_min = VERY_LARGE_NUMBER;
    t_max = -VERY_LARGE_NUMBER;

    t_int = -VERY_LARGE_NUMBER;
    root = NULL;

    // Add nodes to the list.

    for_all_nodes(tdyn, b, bb) {

        // *** Flag recomputation of acc and jerk. ***

        bb->clear_jerk();

        real t = bb->get_time();
        worldline *w = new worldline(bb);

        bundle[nw++] = w;

        t_min = min(t_min, t);
        t_max = max(t_max, t);
    }

    // Sort the list by ID.

    qsort((void*)bundle, (size_t)nw, sizeof(worldlineptr), wlcompare);

//    for (int i = 0; i < nw; i++) {
//      PRC(i); PRC(bundle[i]->get_id());
//      cerr << bundle[i]->get_first_segment()
//                       ->get_first_event()->format_label() << endl;
//    }
}

void worldbundle::print()
{
    // Print out basic information about this worldbundle.

    for (int i = 0; i < nw; i++) {

        worldline *w = bundle[i];
        segment *s = w->get_first_segment();
        int segnum = 0;

        cerr << i << " (" << s->get_first_event()->format_label();
        int p = cerr.precision(HIGH_PRECISION);
        cerr << ", id = " << w->get_id() << ")";
        cerr.precision(p);
        cerr << "  t_start = " << w->get_t_start()
             << ", t_end = " << w->get_t_end()
             << endl;

        while (s) {

            segnum++;
            int nn = 1;

            tdyn *b = s->get_first_event();
            while (b->get_next()) {
                b = b->get_next();
                nn++;
            }

            cerr << "        segment " << segnum << ":  "
                 << s->get_t_start() << " (" << nn << ") " << s->get_t_end()
                 << endl;

            s = s->get_next();
        }
    }
}

void worldbundle::dump(int offset)      // default = 0
{
    for (int i = 0; i < get_nw(); i++) {
        PRI(offset); cerr << "worldline " << i << ":" << endl;
        get_worldline(i)->dump(offset);
    }
}

void worldbundle::check()
{
    cerr << "checking worldbundle " << this << endl;
    for (int i = 0; i < get_nw(); i++)
        get_worldline(i)->check(i);
}

//======================================================================

int worldbundle::find_index(unique_id_t id)
{
    // Find the index of the worldline corresponding to id.

    // Return values:   0 - nw-1        index of id
    //                  -1              id out of range below
    //                  -nw-1           id out of range above
    //                  -nw -- -2       -(first index above id) - 1
    //
    // (so -find_index - 1 is where a new worldline should go in the list).

    // The list is ordered in id; use bisection to search it.

    int loc;

    if (nw <= 0 || id < bundle[0]->get_id()) {
        PRC(id); PRL(bundle[0]->get_id());
        loc = -1;
    }

    else if (id > bundle[nw-1]->get_id())
        loc = -nw-1;

    else {

        int low = 0, high = nw-1;
        loc = (low+high)/2;

        if (bundle[low]->get_id() == id)
            loc = low;

        else if (bundle[high]->get_id() == id)
            loc = high;

        else if (bundle[loc]->get_id() != id) {

            bool found = false;

            while (low < high-1) {

                if (bundle[loc]->get_id() < id)
                    low = loc;
                else if (bundle[loc]->get_id() > id)
                    high = loc;
                else {
                    found = true;
                    break;
                }

                loc = (low+high)/2;
            }

            if (!found) loc = -high - 1;
        }
    }

    return loc;
}

int worldbundle::find_index(char *name) {return find_index(unique_id(name));}
int worldbundle::find_index(_pdyn_ *b)  {return find_index(unique_id(b));}

worldline *worldbundle::find_worldline(unique_id_t id)
{
    int i = find_index(id);
    if (i >= 0 && i < nw)
        return bundle[i];
    else
        return NULL;
}

worldline *worldbundle::find_worldline(char *name)
{return find_worldline(unique_id(name));}

worldline *worldbundle::find_worldline(_pdyn_ *b)
{return find_worldline(unique_id(b));}

//======================================================================

// Update the world bundle.

void worldbundle::attach(tdyn *bn,
                         int verbose)           // default = 0
{
    // Locate and attach all components of bn in the 4tree hierarchy.

    if (bn) {

        for_all_nodes(tdyn, bn, b) {

            // *** Flag recomputation of acc and jerk. ***

            b->clear_jerk();

            // Find the start of b's current worldline segment, or
            // add b to the list.  There are three possibilities:
            //
            //  (1) b's ID is on the list and b is valid
            //                          ==> extend the current segment
            //                              or begin a new one if the
            //                              previous b was defunct
            //
            //  (2) b's ID is on the list and b is defunct
            //                          ==> end the current segment
            //
            //  (3) b's ID is not on the list
            //                          ==> extend the list and open
            //                              a new segment
            //
            // Cases (1) and (2) are handled at the same time.

            // Find b on the worldline list or add it to the list
            // (maintaining the list ordering).

            unique_id_t id = unique_id(b);

            if (id >= 0) {

                real t = b->get_time();
                int loc = find_index(id);

//              int p;
//              bool deb = false;
//              if (t >= 13.275 && t <= 13.313) {
//                  if (node_contains(b, "(12,25)")) {
//                      deb = true;
//                      p = cerr.precision(HIGH_PRECISION);
//                      PRC(b->format_label()); PRL(b->is_defunct());
//                  }
//              }

                if (loc >= 0) {

                    // Extend an existing worldline.

                    worldline *w = bundle[loc];
                    segment *s = w->get_last_segment(); // "last" = "current"

//                  if (deb) {
//                      PRC(s); PRL(t);
//                      PRC(b->format_label()); PRL(id);
//                      w->dump(0, 13.25);
//                  }

                    // Check case (2) for previous instance of b.

                    if (s->get_last_event()->is_defunct()) {

                        if (verbose)
                            cerr << "starting new segment for "
                                 << b->format_label()
                                 << " (id = " << id
                                 << ") at t = " << t
                                 << endl;

                        // Create a new segment starting at b.

                        s = new segment(b);
                        w->add_segment(s, true);  // "true" ==> skip id check
                                                  // (see comments below...)

                        if (verbose)
                            cerr << "new segment = " << s << endl;

                    } else {

                        // Case (1): extend the current segment.

                        s->add_event(b, true);    // "true" ==> skip id check

                        // Note from Steve (8/01):  For unknown reasons,
                        // the id checks in add_event() and add_segment()
                        // may cause errors when optimization is turned on.
                        // Problems arise in multiple systems (so we are
                        // working near the last significant digit), where
                        // the ids may return unequal even though the bits
                        // are identical.  Didn't occur in the July 29,
                        // 2001 version of the code.  Does occur in the
                        // August 2 version!
                        //
                        // Workaround for now: skip the check, which is
                        // redundant here anyway.  Longer-term: improve
                        // the handling of unique_id.  Comparing reals is
                        // probably always a bad idea...
                    }

                    // if (deb) cerr.precision(p);

                    w->set_t_end(t);

                } else {

                    // Case (3): create a new worldline starting at b.
                    // Maintain the ordering of the array by placing the
                    // new pointer at location -loc - 1.

                    if (nw >= nw_max) {

                        // Extend the storage array -- use new and
                        // delete to mimic realloc().

                        nw_max *= 5;
                        nw_max /= 4;

                        cerr << "extending bundle array for worldbundle "
                             << this << ":  nw_max = "
                             << nw_max << endl;

                        worldlineptr *tmp = new worldlineptr[nw_max];
                        for (int i = 0; i < nw; i++)
                            tmp[i] = bundle[i];
                        delete [] bundle;
                        bundle = tmp;
                    }

                    // Where to insert the new worldline:

                    loc = -loc - 1;

                    for (int i = nw-1; i >= loc; i--)
                        bundle[i+1] = bundle[i];

                    bundle[loc] = new worldline(b);
                    nw++;

                    if (verbose)
                        cerr << "adding " << b->format_label()
                             << " (id = " << id
                             << ") at t = " << b->get_time()
                             << endl;
                }

                t_max = max(t_max, t);
            }
        }
    }
}

//======================================================================

static real physical_mass = 0;

real mass_scale_factor()
{
    return physical_mass;
}

local void check_final_escapers(worldbundle *wb, tdyn *b)
{
    for_all_nodes(tdyn, b, bb) {
        worldline *w = NULL;
        if (bb->get_prev()) {                   // esc flag is set
            bb->set_prev(NULL);
            w = wb->find_worldline(bb);
            if (w) w->set_end_esc_flag(1);
        }
        if (bb->get_next()) {                   // t_esc is specified
            real *t_esc = (real*)bb->get_next();
            bb->set_next(NULL);
            if (!w) w = wb->find_worldline(bb);
            if (w) w->set_t_esc(*t_esc);
            delete t_esc;
        }
    }
}

// From Steve (10/00):
//
// Root dumps are done by kira twice per synchronization interval, so each
// worldbundle is self-contained (this is the *same* as the way in which
// tree changes are handled...).

worldbundle *read_bundle(istream &s,
                         int verbose)                   // default = 0
{
    // Read a bundle of worldlines, from root dump to root dump.

    // First read the base root node.  Assume without checking that
    // the next portion of input data starts with a complete system.

    tdyn *b = get_tdyn(s, NULL, NULL, false);           // "stripped tdyn"
    if (!b || !streq(b->format_label(), "root")) return NULL;

    // Need to keep track of the physical initial mass, but this only
    // appears in the first Log story and is not saved in the usual
    // story structure, to save space.  Grab it here and make it
    // globally accessible.

    if (b->get_log_story())
        physical_mass = getrq(b->get_log_story(), "physical_initial_mass");

    // Create the initial list of node IDs, esc_flags, etc.  Note that
    // we use the prev pointer in tdyn to carry temporary information
    // about cluster membership.  These should be new worldline(b),
    // which is called by new worldbundle(b).  We do *not* check here!

    worldbundle *wb = new worldbundle(b);

    if (verbose)
        cerr << endl
             << "new worldbundle: created initial list, nw = "
             << wb->get_nw() << endl;

    // Now read in individual nodes and attach them to the 4tree.
    // Stop when another tree is read in.

    while (b = get_tdyn(s, NULL, NULL, false)) {        // "stripped tdyn"

        // A new tree fragment will be read in every time the tree
        // changes, so pointers should be correctly maintained.
        // Links are found at the start of each worldline segment.

        // If this is the terminating full snapshot, we must check
        // the prev pointers for updated information on escapers
        // before we attach to the 4tree.

        bool stop = (b->get_oldest_daughter()
                     && streq(b->format_label(), "root"));
        if (stop)
            check_final_escapers(wb, b);

        // Connect b to the 4tree.

        wb->attach(b, verbose);

        // Stop once we have read in another complete system.

        if (stop) {
            if (verbose)
                cerr << "break at t = " << b->get_time() << endl;
            break;
        }
    }

    if (verbose) {

        cerr << "after last input:  nw = " << wb->get_nw()
             << endl;

        // Print some global diagnostics:

        for (int i = 0; i < wb->get_nw(); i++) {
            worldline *w = wb->get_bundle()[i];

            // First check segments.

            int seg = 0, nseg = 0;
            segment *s = w->get_first_segment();
            tdyn *b = s->get_first_event();
            real t = s->get_t_start();

            while (s) {
                nseg++;
                s = s->get_next();
            }

            s = w->get_first_segment();
            int n_jump = 0;

            while (s) {

                if (s->get_t_start() >= s->get_t_end()) {
                    int p = cerr.precision(HIGH_PRECISION);
                    cerr << "worldline " << i << " (" << b->format_label()
                         << "), id = " << w->get_id() << endl;
                    cerr.precision(p);
                    PRI(10); cerr << "zero-length segment at t = "
                                  << s->get_t_start()
                                  << ";  segment " << seg << " of " << nseg
                                  << endl;
                }

                if (s->get_t_start() > t) {
                    n_jump++;
                    if (verbose > 1) {
                        int p = cerr.precision(HIGH_PRECISION);
                        cerr << "worldline " << i << " (" << b->format_label()
                            << "), id = " << w->get_id() << endl;
                        cerr.precision(p);
                        PRI(10); cerr << "jump from " << t
                                      << " to " << s->get_t_start()
                                      << " after segment " << seg
                                      << " of " << nseg
                                      << endl;
                    }
                }

                t = s->get_t_end();
                seg++;
                s = s->get_next();
            }

            if (verbose == 1 && n_jump > 0) {
                cerr << "worldline " << i << " (" << b->format_label()
                     << "), id = " << w->get_id()
                     << " has " << n_jump << " jump";
                if (n_jump > 1) cerr << "s";
                cerr << endl;
            }

            // Then check events within each segment.

            seg = 0;
            s = w->get_first_segment();
            while (s) {

                b = s->get_first_event();

                if (b->get_time() != s->get_t_start()) {
                    int p = cerr.precision(HIGH_PRECISION);
                    cerr << "worldline " << i << " (" << b->format_label()
                         << "), id = " << w->get_id() << endl;
                    cerr.precision(p);
                    PRI(10); cerr << "t_start =  " << s->get_t_start()
                                  << ", t = "
                                  << b->get_time() << " in segment " << seg
                                  << " of " << nseg << endl;
                }

                while (b && b->get_next()) b = b->get_next();

                if (b->get_time() != s->get_t_end()) {
                    int p = cerr.precision(HIGH_PRECISION);
                    cerr << "worldline " << i << " (" << b->format_label()
                         << "), id = " << w->get_id() << endl;
                    cerr.precision(p);
                    PRI(10); cerr << "t_end =  " << s->get_t_end() << ", t = "
                                  << b->get_time() << " in segment " << seg
                                  << " of " << nseg << endl;
                }

                seg++;
                s = s->get_next();
            }
        }
    }

    return wb;
}

int count_segments(worldbundle *wb)
{
    int ns = 0;
    for (int i = 0; i < wb->get_nw(); i++) {
        worldline *w = wb->get_worldline(i);
        segment *s = w->get_first_segment();
        while (s) {
            ns++;
            s = s->get_next();
        }
    }
    return ns;
}

int count_events(worldbundle *wb)
{
    int ne = 0;
    for (int i = 0; i < wb->get_nw(); i++) {
        worldline *w = wb->get_worldline(i);
        segment *s = w->get_first_segment();
        while (s) {
            tdyn *b = s->get_first_event();
            while (b) {
                ne++;
                b = b->get_next();
            }
            s = s->get_next();
        }
    }
    return ne;
}

//======================================================================

// Navigation of the 4tree data and diagnostic output.

tdyn *find_event(worldline *w, tdyn *bn, real t)
{
    // Find time t along the segment of worldline w starting at base
    // node bn.  Return a pointer to the event immediately preceding t.
    //
    // If worldline w is set, first see if any current_event pointer
    // exists; if so, start our search there.  The logic in the calling
    // function should ensure that time t lies in the current segment.

    // Default is forward search starting at bn.  This may be modified
    // if current_event information exists.

    tdyn *b = bn;

    if (w) {

        tdyn *curr = w->get_current_event();

        if (curr) {

            real t_curr = curr->get_time();
            real t_bn = bn->get_time();         // evaluating bn->get_time() in
                                                // if (..) below causes strange
                                                // error on margaux (DecUNIX)...

            if (t_curr > t) {

                if (t > 0.5*(t_bn + t_curr)) {

                    // Do a backward search from curr.

                    b = curr;
#if 0
                    if (b != bn && !b->get_prev()) {
                        cerr << endl << "prev = NULL for node " << b << endl;
                        PRC(b); PRL(bn);
                        cerr << "worldline dump:" << endl;
                        w->dump();
                        cerr << "checking worldline..." << endl;
                        w->check(-1);
                    }
#endif

                    tdyn *p;
                    while (b && (p = b->get_prev()) && b->get_time() > t)
                        b = p;

                    return b;
                }

            } else {

                // Explicitly test the most obvious possibility first.
                // Case next = NULL probably means an error in the data,
                // but do something reasonable anyway (return curr).

                tdyn *n = curr->get_next();
                if (!n || n->get_time() >= t) return curr;

                // Do a forward search from n.

                b = n;

            }
        }
    }

    // Do a forward search from b.  Note that repeated "get_next" calls
    // may be quite time consuming if we keep going outside the cache.

    tdyn *n;
    while (b && (n = b->get_next()) && n->get_time() < t)
        b = n;

    // The portion of the trajectory spanning t runs from b to b->next.

    return b;
}

void print_event(worldline *w, tdyn *bn, real t)
{
    // Print info on the portion of the worldline portion
    // starting at bn that spans t.

    tdyn *b = find_event(w, bn, t);

    if (b) {
        PRC(t); PRC(b->get_time()); PRL(b->get_next());
        if (b->get_next()) PRL(b->get_next()->get_time());
    }
}

//----------------------------------------------------------------------
#include "print_local.C"        // avoid repeating local print functions
//----------------------------------------------------------------------

local inline bool check_and_initialize(tdyn *p,
                                       real tp,
                                       real t,
                                       vector& pos,
                                       bool inc,
                                       tdyn *bn)
{
    if (tp < t) {

        tdyn *n = p->get_next();

        if (!n) {

            cerr << "interpolate_pos: error 2: next node not found: ";
            PRC(p->format_label()); PRL(t);
            int pr = cerr.precision(HIGH_PRECISION);
            PRL(unique_id(p)); cerr.precision(pr);
            PRI(26); PRC(bn); PRL(bn->get_time());
            PRI(26); PRC(p); PRL(p->get_time());
            PRI(26); PRL(p->get_time()-t);
            if (bn) {
                PRI(26); PRL(bn->format_label());
                pr = cerr.precision(HIGH_PRECISION);
                PRL(unique_id(bn)); cerr.precision(pr);
            }
            // print_details(wb, p, t);

            if (inc)
                pos += p->get_pos();
            else
                pos = p->get_pos();
            return false;

        } else if (n->get_time() < t) {

            cerr << "interpolate_pos: error 3: specified time above range: ";
            PRC(p->format_label()); PRC(t); PRL(p->get_time());
            // print_details(wb, p, t);

            if (inc)
                pos += p->get_pos();
            else
                pos = p->get_pos();
            return false;
        }

        // We now have tp < t <= n->time.  We will interpolate using
        // a fit to pos and vel.  Check to see if the interpolating
        // coefficients must be set.

        if (p->get_jerk()[0] == 0) {

            // Recompute acc/2 and jerk/6 to fit pos and vel.

            // Note that we overwrite acc and jerk by equivalent
            // vectors that guarantee continuity of pos and vel.

            // *** Flag this to prevent recalculation by setting ***
            // *** jerk = 0 as the tree is constructed.          ***

            // Set up the interpolating acc and jerk.

            real dt = n->get_time() - tp;
            real dti = 1/dt;

            if (streq(p->get_name(), "root")) {         // infrequent...

                // Use linear interpolation for the cluster center.
                // Discard vel information to make pos continuous.

                p->set_vel(dti*(n->get_pos() - p->get_pos()));
                p->set_acc(0);
                p->set_jerk(vector(VERY_SMALL_NUMBER, 0, 0));

            } else {

                // Standard fourth-order interpolation otherwise.

                p->set_acc((3*(n->get_pos()-p->get_pos())
                            - dt*(2*p->get_vel()+n->get_vel()))*dti*dti);
                p->set_jerk((p->get_vel()+n->get_vel()
                             - 2*dti*(n->get_pos()-p->get_pos()))*dti*dti);
            }
        }
        return true;

    } else if (tp > t) {

        cerr << "interpolate_pos: error 1: specified time below range: ";
        PRC(p->format_label()); PRC(t); PRL(p->get_time());
        // print_details(wb, p, t); // wb not available here...
    }

    // (Special case tp == t is handled here.)

    if (inc)
        pos += p->get_pos();
    else
        pos = p->get_pos();

    return false;
}

// Interpolation of part of the 4tree to a specific time.

// Retain numerous "interpolate_pos" functions for timing and testing...
// Apparently very little time difference between them.  New code may be
// marginally faster, by perhaps 5 percent, but its use unaccountably
// increases the time spent in other portions of the code, for a net
// improvement of only 1 percent or so.  See profile_*_interp.txt.  The
// function node::next_node() seems to be used 4 times more frequently
// in the new code, offsetting most of the gain in the other functions.

// Old version:

vector interpolate_pos(tdyn *p, real t,
                       tdyn *bn)                // default = NULL; specifies
                                                // actual base node
{
    // The range (p to p->next) includes time t.
    // Interpolate and return the pos of this particle.

    // Check...

    if (p->get_time() > t) {
        cerr << "interpolate_pos: error 1: specified time below range: ";
        PRC(p->format_label()); PRC(t); PRL(p->get_time());
        // print_details(wb, p, t); // wb not available here...

        return p->get_pos();
    }

    // Special case:

    if (p->get_time() == t) return p->get_pos();

    tdyn *n = p->get_next();

    if (!n) {
        cerr << "interpolate_pos: error 2: next node not found: ";
        PRC(p->format_label()); PRL(t);
        int pr = cerr.precision(HIGH_PRECISION);
        PRL(unique_id(p)); cerr.precision(pr);
        PRI(26); PRC(bn); PRL(bn->get_time());
        PRI(26); PRC(p); PRL(p->get_time());
        PRI(26); PRL(p->get_time()-t);
        if (bn) {
            PRI(26); PRL(bn->format_label());
            pr = cerr.precision(HIGH_PRECISION);
            PRL(unique_id(bn)); cerr.precision(pr);
        }
        // print_details(wb, p, t);

        return p->get_pos();
    }

    if (n->get_time() < t) {
        cerr << "interpolate_pos: error 3: specified time above range: ";
        PRC(p->format_label()); PRC(t); PRL(p->get_time());
        // print_details(wb, p, t);

        return p->get_pos();
    }

    // Time t is included in the range.

    // Interpolate using pos and vel for now...
    // Note that we overwrite acc and jerk by equivalent
    // vectors that guarantee continuity of pos and vel.

    real tp = p->get_time();
    real dt = n->get_time() - tp;

    if (dt > 0) {

        // Recompute acc/2 and jerk/6 to fit pos and vel.

        // *** Flag this to prevent recalculation by setting ***
        // *** jerk = 0 as the tree is constructed.          ***

        if (p->get_jerk()[0] == 0) {

            // Set up the interpolating acc and jerk.

            real dti = 1/dt;

            if (streq(p->get_name(), "root")) {         // infrequent...

                // Use linear interpolation for the cluster center.
                // Discard vel information to make pos continuous.

                p->set_vel(dti*(n->get_pos() - p->get_pos()));
                p->set_acc(0);
                p->set_jerk(vector(VERY_SMALL_NUMBER, 0, 0));

            } else {

                // Standard fourth-order interpolation otherwise.

                p->set_acc((3*(n->get_pos()-p->get_pos())
                            - dt*(2*p->get_vel()+n->get_vel()))*dti*dti);
                p->set_jerk((p->get_vel()+n->get_vel()
                             - 2*dti*(n->get_pos()-p->get_pos()))*dti*dti);
            }
        }

        dt = t - tp;

        return p->get_pos() + dt * (p->get_vel()
                                    + dt * (p->get_acc()
                                            + dt * p->get_jerk()));
    } else

        return p->get_pos();
}

// Overloaded (and not used?):

void interpolate_pos(tdyn *p,
                     real t,
                     vector& pos,
                     bool inc,                  // default = false
                     tdyn *bn)                  // default = NULL; specifies
                                                // the actual base node
{
    // The range (p to p->next) includes time t.
    // Interpolate and return the pos of this particle.

    real tp = p->get_time();

    if (check_and_initialize(p, tp, t, pos, inc, bn)) {

        real dt = t - tp;

        if (inc)
            pos += p->get_pos() + dt * (p->get_vel()
                                           + dt * (p->get_acc()
                                                      + dt * p->get_jerk()));
        else
            pos = p->get_pos() + dt * (p->get_vel()
                                          + dt * (p->get_acc()
                                                     + dt * p->get_jerk()));
    }
}

// New versions (NEW_INTERP):

void set_interpolated_pos(tdyn *p,
                          real t,
                          vector& pos,
                          tdyn *bn)             // default = NULL; specifies
                                                // the actual base node
{
    // The range (p to p->next) includes time t.
    // Interpolate and set the value of pos.

    real tp = p->get_time();

    if (check_and_initialize(p, tp, t, pos, false, bn)) {

        real dt = t - tp;
        pos = p->get_pos() + dt * (p->get_vel()
                                      + dt * (p->get_acc()
                                                 + dt * p->get_jerk()));
    }
}

// Overloaded:

void set_interpolated_pos(tdyn *p,
                          real t,
                          pdyn *curr,
                          tdyn *bn)             // default = NULL; specifies
                                                // the actual base node
{
    // The range (p to p->next) includes time t.
    // Interpolate and set the value of pos.

    real tp = p->get_time();
    vector pos;

    if (check_and_initialize(p, tp, t, pos, false, bn)) {

        real dt = t - tp;
        curr->set_pos(p->get_pos() + dt * (p->get_vel()
                                            + dt * (p->get_acc()
                                                     + dt * p->get_jerk())));
    } else

        curr->set_pos(pos);
}

void inc_interpolated_pos(tdyn *p,
                          real t,
                          vector& pos,
                          tdyn *bn)             // default = NULL; specifies
                                                // the actual base node
{
    // The range (p to p->next) includes time t.
    // Interpolate and increment the value of pos.

    real tp = p->get_time();

    if (check_and_initialize(p, tp, t, pos, true, bn)) {

        real dt = t - tp;
        pos += p->get_pos() + dt * (p->get_vel()
                                       + dt * (p->get_acc()
                                                  + dt * p->get_jerk()));
    }
}

vector interpolate_vel(tdyn *p, real t,
                       tdyn *bn)                // default = NULL; specifies
                                                // actual base node
{
    // Interpolate and return the vel of this particle.

    // Only called after a call to interpolate_pos, so skip checks
    // and recomputation of acc and jerk.

    // Special case:

    if (p->get_time() == t) return p->get_vel();

    tdyn *n = p->get_next();
    real tp = p->get_time();
    real dt = n->get_time() - tp;

    if (dt > 0) {

        dt = t - tp;
        return p->get_vel() + dt * (2*p->get_acc()
                                    + dt * 3*p->get_jerk());
    } else

        return p->get_vel();
}

// Older versions of interpolate_pos() and interpolate_vel() are
// (for now) saved in old_interpolate.C

local vector get_pos(worldline *w,
                     tdyn *b, tdyn *bn,
                     real t = VERY_LARGE_NUMBER)
{
    // Return the current absolute position of node b, properly
    // corrected for its location in the tree.

    // The node of interest is b; the base node for the segment of
    // worldline w (containing all parental information) is bn.

    if (t == -VERY_LARGE_NUMBER) t = b->get_time();

    vector pos;
    set_interpolated_pos(b, t, pos, bn);

    // The tree structure associated with the base node should extend
    // to the top level, but will not contain the root node unless bn
    // is part of a full dump.

    // Ascend the base tree, including parent contributions (at time t)
    // as needed.  Exclude root, but note that is_root() will fail, as
    // it simply checks for a null parent.  For now, check the name
    // explicitly...

    tdyn *bp = bn->get_parent();

    while (bp) {

        // Find the portion of the parent worldline spanning t.

        tdyn *p = find_event(w, bp, t);

        inc_interpolated_pos(p, t, pos, bn);
        bp = bp->get_parent();
    }

    return pos;
}

void worldbundle::print_worldline(char *name,
                                  real dt)      // default = 0
{
    // Print out the entire worldline (all segments) of particle 'name'.
    // Take steps of length dt (0 ==> just use the stored times).

    int loc = find_index(name);                 // find the worldline
                                                // for particle 'name'
    if (loc >= 0) {

        worldline *w = bundle[loc];
        segment *s = w->get_first_segment();

        if (dt == 0) {

            while (s) {                             // loop over segments

                tdyn *bn = s->get_first_event();    // starting point
                tdyn *b = bn;
                while (b) {                         // loop over events
                    cout << "    " << b->get_time()
                         << " " << get_pos(w, b, bn) << endl << flush;
                    b = b->get_next();
                }

                s = s->get_next();
            }

        } else {

            real t = get_t_min();
            while (t < get_t_max() - 0.5*dt) {

                // Code here assumes that we are moving sequentially
                // through the data, as we start with the previous s.

                while (s && s->get_t_end() < t) s = s->get_next();

                if (s) {
                    tdyn *bn = s->get_first_event();
                    tdyn *b = find_event(w, bn, t);
                    cout << "    " << t
                         << " " << get_pos(w, b, bn) << endl << flush;
                }

                t += dt;
            }
        }

    } else
        cerr << name << " not found" << endl;
}

//======================================================================

// Creation of an entire interpolated tree at some specific time.

//----------------------------------------------------------------------
#include "print_local2.C"       // avoid repeating local print functions
//----------------------------------------------------------------------

//----------------------------------------------------------------------
#include "attach_new_node.C"    // avoid repeating local print functions
//----------------------------------------------------------------------

//----------------------------------------------------------------------
#include "update_node.C"        // avoid repeating local functions
//----------------------------------------------------------------------

local INLINE void add_to_interpolated_tree(worldbundle *wb,
                                           worldline *w, segment *s,
                                           pdyn *root, real t, bool vel,
                                           bool debug)
{
    // Current worldline is w, segment is s.  Find the base node
    // (containing all relevant tree information) for this particle,
    // and the current event.

    // *** By construction, on entry, w and s represent a leaf. ***

    tdyn *bn = s->get_first_event();

    if (debug) {
        cerr << "s = " << s << endl
             << "first event = " << bn << " "
             << bn->format_label() << " "
             << bn->get_time() << endl;

        print_details(wb, bn, t);
    }

    // Create the entire portion of the tree corresponding to the
    // top-level node of bn.  Note that top will probably be bn most
    // of the time.

    // Need to be careful with top-level nodes, as they  generally
    // won't have parent nodes (no root), so standard functions like
    // "is_top_level_node()" and "get_top_level_node()" may fail.
    // (Choice of root in calling function should correct this...)

    tdyn *top = bn;
    while (top->get_parent()
           && unique_id(top->get_parent()) > 0)
        top = top->get_parent();

    if (debug) {
        PRC(top); PRL(top->format_label());
        // put_node(cerr, *top, false);
    }

    // Copy the entire tree below top to the new tree.

    for_all_nodes(tdyn, top, bb) {

        // Find the worldline of bb.

        worldline *ww;

        if (bb == bn)
            ww = w;
        else {
            ww = wb->find_worldline(bb);
            if (!ww)
                cerr << "create_interpolated_tree: error:"
                     << "can't find worldline of subtree member." << endl;
        }

        if (ww) {

            // Create and attach a new node.

            attach_new_node(wb, ww, root, top, bb, debug);

            // Update the new tree entry: copy or interpolate all
            // relevant quantities from bb to curr.

            update_node(wb, ww, t, bb, top, vel, debug);
        }
    }
}

int id_n_clump(unique_id_t id)                  // number of particles in id
{
#ifndef NEW_UNIQUE_ID_T
    return (int) id;
#else
    return (id>>ID_SHIFT);                      // new new (12/01)
#endif
}

local inline bool id_is_leaf(unique_id_t id)
{
    return (id_n_clump(id) == 1);
}

#define EPS 1.e-12

pdyn *create_interpolated_tree(worldbundle *wb, real t,
                               bool vel)                // default = false
{
    // All leaves on the bundle list are still current, by construction.
    // Find the base segment corresponding to each and use it to build
    // a new pdyn tree interpolated to the current time.

    bool debug = false;

    // Try to take care of rounding error in t.  Management of worldbundles
    // should be the responsibility of the calling program.  (This would not
    // be an issue if timesteps were constrained to be powers of 2...)

    real dt = t - wb->get_t_max();
    if (dt > EPS)
        return NULL;
    else if (dt > 0)
        t = wb->get_t_max();

    dt = t - wb->get_t_min();
    if (dt < -EPS)
        return NULL;
    else if (dt < 0)
        t = wb->get_t_min();

    // Find the array of worldlines (worldbundle), and initialize
    // the tree_node pointers.

    worldlineptr *bundle = wb->get_bundle();

    for (int i = 0; i < wb->get_nw(); i++)
        bundle[i]->clear_tree_node();

    // Create a root node.

    pdyn *root = new pdyn(NULL, NULL, false);

    // Replace by alloc_pdyn() if we do our own memory management.

    root->set_system_time(t);
    root->set_pos(0);

    // Establish root as the global root node (should clean up problems
    // with "top_level_node" functions...).  Stored in location 0 of the
    // worldline array.

    root->set_root(root);
    bundle[0]->set_tree_node(root);

    // Loop through remaining worldlines and take action for leaves only.
    // Logic: we are really dealing with top-level nodes, but we don't know
    // which component we will encounter first.  When we find a component
    // of a clump of particles, we update the entire clump, and flag all
    // components accordingly.

    for (int i = 1; i < wb->get_nw(); i++) {

        worldline *w = bundle[i];

        if (!w->get_tree_node()) {

            // Tree entry for this node has not yet been created/updated.

            unique_id_t id = w->get_id();

//          if (id >= 1 && id < 2) {
            if (id_is_leaf(id)) {

                // Worldline w represents a leaf.  Locate its current segment.

                if (debug) {
                    cerr.precision(16);
                    cerr << endl
                         << "following worldline " << i << ", id = "
                         << w->get_id() << endl;
                }

                segment *s = w->get_first_segment();
                while (s && s->get_t_end() < t) s = s->get_next();

                if (debug) {
                    PRC(w); cerr << "segment "; PRC(s);
                    print_events(s, t);
                }

                // If s is NULL, it should mean that w refers to a CM node
                // that does not exist at time t.  Particle worldlines
                // should be continuous -- could check in debugging mode...

                if (s) add_to_interpolated_tree(wb, w, s,
                                                root, t, vel, debug);
            }
        }
    }

    return root;
}

#else

main(int argc, char *argv[])
{
    // Fixed format on the command line:  filename  time.

    if (argc <= 2) exit(1);

    ifstream s(argv[1]);
    if (!s) exit(2);

    real t = atof(argv[2]);

    worldbundle *wb = read_bundle(s);

    pdyn *root = create_interpolated_tree(wb, t, true);
    put_node(cout, *root, false);
    rmtree(root);

#if 0
    // Some statistics:

    wb->print();

    cerr << wb->get_nw() << " worldlines, "
         << count_segments(wb) << " segments, "
         << count_events(wb) << " events, t = "
         << wb->get_t_min() << " to " << wb->get_t_max()
         << endl;
#endif

#if 0
    // Locate a specific time along the worldline of some particle.

    char name[128];
    real t;
    while (1) {
        cerr << endl << "time, name: "; cin >> t >> name;
        if (cin.eof()) {
            cerr << endl;
            break;
        }
        int loc = wb->find_index(name);
        if (loc >= 0) {
            PRC(unique_id(name)); PRL(loc);
            worldline *w = wb->get_bundle()[loc];
            segment *s = w->get_first_segment();
            while (s && s->get_t_start() > t) s = s->get_next();
            print_event(w, s->get_first_event(), t);
        } else
            cerr << "not found" << endl;
    }
#endif

#if 0
    // Print out the entire worldline of a specified particle,
    // taking steps of the specified length.

    real dt;
    char name[128];
    while (1) {

        cerr << endl << "dt, name: " << flush; cin >> dt >> name;
        if (cin.eof()) {
            cerr << endl;
            break;
        }

        wb->print_worldline(name, dt);
    }
#endif

#if 0
    real t = 0.8;
    while (t < 0.85) {
        pdyn *root = create_interpolated_tree(wb, t);
        put_node(cout, *root, false);
        rmtree(root);
        t += 0.01;
    }
#endif

}

#endif

---