hdyn_tt.C Source File

00001 
00002        //=======================================================//    _\|/_
00003       //  __  _____           ___                    ___       //      /|\ ~
00004      //  /      |      ^     |   \  |         ^     |   \     //          _\|/_
00005     //   \__    |     / \    |___/  |        / \    |___/    //            /|\ ~
00006    //       \   |    /___\   |  \   |       /___\   |   \   // _\|/_
00007   //     ___/   |   /     \  |   \  |____  /     \  |___/  //   /|\ ~
00008  //                                                       //            _\|/_
00009 //=======================================================//              /|\ ~
00010 
00011 // hdyn_tt:  Starlab hdyn-specific node-handling functions.
00012 
00013 // Externally visible functions:
00014 //
00015 //      hdyn* hdyn::next_node
00016 //      hdyn* hdyn::find_perturber_node
00017 //      bool hdyn::nn_stats
00018 //      real hdyn::print_pert
00019 //      void hdyn::setup_binary_node
00020 //
00021 //      void print_nn
00022 //      void print_coll
00023 
00024 #include "hdyn.h"
00025 
00026 #ifndef TOOLBOX
00027 
00028 void hdyn::null_pointers()
00029 {
00030     // Clear all pointers (don't touch what they point to)
00031     // -- for use in cleaning up temporary nodes...  Careful!!
00032 
00033     nn = coll = NULL;
00034     n_perturbers = 0;
00035     perturber_list = NULL;
00036     valid_perturbers = false;
00037 
00038     _dyn_::null_pointers();
00039 }
00040 
00041 // Note from Steve 8/20/98.  This is based on the node version, but is
00042 // modified to include the possiblilty of treating unperturbed centers
00043 // of mass as leaves.
00044 
00045 hdyn* hdyn::next_node(hdyn* base,
00046                       bool resolve_unperturbed) // default = true
00047 {
00048     // If resolve_unperturbed is false, we only descend to the next level
00049     // in the tree if the oldest daughter is perturbed (Steve, 8/20/98).
00050 
00051     if (oldest_daughter
00052         && (resolve_unperturbed
00053             || get_oldest_daughter()->kep == NULL
00054             || !get_oldest_daughter()->fully_unperturbed)) {
00055 
00056         return get_oldest_daughter();
00057 
00058     } else if (this == base) {          // in case base is a leaf...
00059 
00060         return NULL;
00061 
00062     } else if (younger_sister) {
00063 
00064         return get_younger_sister();
00065 
00066     } else {
00067 
00068         // Can't go down or right.  See if we can go up.
00069 
00070         if (this == base)               // in case base is a leaf...
00071             return NULL;
00072 
00073         if (parent == NULL)             // in case b is atomic...
00074             return NULL;
00075 
00076         hdyn* tmp = get_parent();
00077         while (tmp->get_younger_sister() == NULL) {
00078             if (tmp == base)
00079                 return NULL;
00080             else
00081                 tmp = tmp->get_parent();
00082         }
00083 
00084         // Now tmp is the lowest-level ancestor with a younger sister.
00085 
00086         if (tmp == base)
00087             return NULL;
00088         else
00089             return tmp->get_younger_sister();
00090     }
00091 
00092     return NULL;        // To keep some compilers happy...
00093 }
00094 
00095 // find_perturber_node:  return a pointer to the first ancestor node
00096 //                       found that contains a valid perturber list,
00097 //                       or NULL otherwise.
00098 
00099 hdyn* hdyn::find_perturber_node()
00100 {
00101     if (is_root() || is_top_level_node()) return NULL;
00102 
00103     hdyn* pnode;
00104     hdyn* top_level = get_top_level_node();
00105 
00106     if (!ALLOW_LOW_LEVEL_PERTURBERS) {
00107 
00108         // Same result, so avoid the extra work in this case.
00109 
00110         pnode = top_level;
00111 
00112     } else {
00113 
00114         pnode = get_parent();
00115         while (pnode != top_level && !pnode->valid_perturbers)
00116             pnode = pnode->get_parent();
00117 
00118     }
00119 
00120     if (!pnode->valid_perturbers) pnode = NULL;
00121     return pnode;
00122 }
00123 
00124 // print_nn: print out the ID of b's nearest neighbor.
00125 
00126 void print_nn(hdyn* b,
00127               int level,        // default = 0
00128               ostream& s)       // default = cerr
00129 {
00130     if (level > 0)
00131         s << "nn of " << b->format_label() << " is ";
00132 
00133     if (b->get_nn()) {
00134         s << "(" << b->get_nn() << ") ";
00135         if (b->get_nn()->is_valid()) {
00136             s << b->get_nn()->format_label();
00137             if (level > 1)
00138                 s << " (d = " << sqrt(b->get_d_nn_sq()) << ")";
00139         } else
00140             s << "invalid";
00141     } else
00142         s << "(NULL)";
00143 
00144     if (level > 0)
00145         s << endl;
00146 }
00147 
00148 // print_coll: print out the ID of b's coll particle.
00149 
00150 void print_coll(hdyn* b,
00151                 int level,      // default = 0
00152                 ostream& s)     // default = cerr
00153 {
00154     if (level > 0)
00155         s << "coll of " << b->format_label() << " is ";
00156 
00157     if (b->get_coll()) {
00158         s << "(" << b->get_coll() << ") ";
00159         if (b->get_coll()->is_valid()) {
00160             s << b->get_coll()->format_label();
00161             if (level > 1)
00162                 s << " (d = " << sqrt(b->get_d_coll_sq())
00163                   << ", R = " << b->get_coll()->get_radius()
00164                   << ")";
00165         } else
00166             s << "invalid";
00167     } else
00168         s << "(NULL)";
00169 
00170     if (level > 0)
00171         s << endl;
00172 }
00173 
00174 bool hdyn::nn_stats(real energy_cutoff, real kT,
00175                     vector center, bool verbose,
00176                     bool long_binary_output,            // default = true
00177                     int which)                          // default = 0
00178 {
00179     if (which == 0 && verbose)
00180         cerr << "\n  Bound nn pairs:\n";
00181 
00182     if (!nn) return false;
00183 
00184     // Avoid printing mutual nearest neighbors twice, and print
00185     // the data from the primary (where possible):
00186 
00187     if (nn->nn == this) {
00188         if (mass < nn->mass)                        // fails for equal masses
00189             return false;
00190         else if (mass == nn->mass) {
00191 
00192             // The following *can* occur in the GRAPE-4 version:
00193 
00194 #if 0
00195             if (nn == this)
00196                 cerr << "    warning: " << format_label()
00197                      << " is its own nearest neighbor..." << endl;
00198 #endif
00199 
00200             if (this >= nn)                         // always legal?
00201                 return false;                       // compare labels if not...
00202         }
00203     }
00204 
00205     real M = mass + nn->mass;
00206     if (M <= 0) return false;
00207 
00208     vector dx, dv;
00209 
00210     if (parent == nn->parent) {
00211 
00212         dx = pos - nn->pos;
00213         dv = vel - nn->vel;
00214 
00215     } else {
00216 
00217         // Won't happen if nn search is limited to top-level nodes
00218         // (and probably unnecessary otherwise, as binary should
00219         // already have been picked up)...
00220 
00221         dx = hdyn_something_relative_to_root(this, &hdyn::get_pos)
00222            - hdyn_something_relative_to_root(nn, &hdyn::get_pos);
00223         dv = hdyn_something_relative_to_root(this, &hdyn::get_vel)
00224            - hdyn_something_relative_to_root(nn, &hdyn::get_vel);
00225 
00226     }
00227 
00228     real mu = mass * nn->mass / M;
00229     real E = mu*(0.5*dv*dv - M / abs(dx));
00230 
00231     // Convention: energy_cutoff is in terms of kT if set,
00232     //             in terms of E/mu if kT = 0.
00233 
00234     bool found = false;
00235 
00236     if ((kT > 0 && E < -energy_cutoff*kT)
00237         || (kT == 0 && E <  -energy_cutoff*mu)) {
00238 
00239         // cerr << " nn of " << format_label();
00240         // cerr << " is " << nn->format_label() << endl << flush;
00241 
00242         print_binary_from_dyn_pair(this, nn, kT, center,
00243                                    verbose, long_binary_output);
00244         cerr << endl;
00245 
00246         found = true;
00247     }
00248 
00249     return found;
00250 }
00251 
00252 real hdyn::print_pert(bool long_binary_output,          // default = true
00253                       int indent)                       // default = BIN_INDENT
00254 {
00255     // Print out information on perturbers and nn of this low-level node.
00256     // Formatting is to match sys_stats output.
00257     // Return zero or the energy of this binary, if unperturned.
00258 
00259     // NOTE: 'this' is a binary component.
00260 
00261     if (long_binary_output) {
00262 
00263         hdyn* parent = get_parent();
00264         hdyn* nn = parent->get_nn();
00265 
00266         // cerr << endl;                // omitted by print_binary_params...
00267 
00268         if (timestep > 0 && parent->timestep > 0) {
00269             PRI(indent);
00270             cerr << "cpt timestep = " << timestep
00271                  << "  CM timestep = " << parent->timestep << endl;
00272         }
00273 
00274         PRI(indent);
00275         real pert_sq = 0;
00276         if(get_perturbation_squared()>=0)
00277           pert_sq = sqrt(get_perturbation_squared());
00278         cerr << "pert = " << pert_sq;
00279 
00280         if (slow) cerr << " [" << get_kappa() << "]";
00281 
00282         cerr << " (";
00283         hdyn* pnode = find_perturber_node();
00284         if (pnode && pnode->get_valid_perturbers())
00285             cerr << pnode->get_n_perturbers()
00286                  << " perturbers)";
00287         else
00288             cerr << "no perturber list)";
00289 
00290         cerr << "  nn is ";
00291         if (nn == NULL || nn == parent)
00292             cerr << "unknown" << endl;
00293         else {
00294             cerr << nn->format_label() << endl;
00295             real rnn = abs(parent->get_pos() - nn->get_pos());
00296             if (rnn > 0) {
00297                 real energy =
00298                     -(parent->get_mass() + nn->get_mass()) / rnn
00299                         + 0.5 * square(parent->get_vel()
00300                                        - nn->get_vel());
00301 
00302                 // Some of this information may already have been
00303                 // printed out for bound nn pairs, but repeat it
00304                 // here for uniformity.
00305 
00306                 // if (energy >= 0) {
00307                 PRI(indent);
00308                 cerr << "nn mass = " << nn->get_mass()
00309                      << "  dist = " << rnn
00310                      << "  E/mu = " << energy << endl;
00311                 // }
00312             }
00313         }
00314 
00315     } else {
00316 
00317         if (slow) cerr << " [" << get_kappa() << "]";
00318         cerr << endl;
00319 
00320     }
00321 
00322     real e_unp = 0;
00323 
00324     if (kep)
00325         e_unp -= mass * get_binary_sister()->mass
00326                       / (2 * kep->get_semi_major_axis());
00327 
00328     return e_unp;
00329 }
00330 
00331 #define TOLERANCE (1e-12)
00332 
00333 inline int within_tolerance(real x, real scale)
00334 {
00335     if (scale != 0.0) {
00336         return abs(x) / abs(scale) < TOLERANCE;
00337     } else {
00338         return abs(x) < TOLERANCE;
00339     }
00340 }
00341 
00342 void check_consistency_of_node(hdyn * node,
00343                                hdyn_VMF_ptr get_something,
00344                                char *id)
00345 {
00346     hdyn *first_child = node->get_oldest_daughter();
00347     hdyn *second_child = first_child->get_younger_sister();
00348 
00349     if (second_child->get_younger_sister() != NULL) {
00350         node->pretty_print_node(cerr);
00351         cerr << " is not a binary node \n";
00352         return;
00353     }
00354     vector d = (first_child->*get_something) ()
00355                 - (second_child->*get_something) ();
00356     vector cm = (first_child->*get_something) () * first_child->get_mass()
00357                 + (second_child->*get_something) () * second_child->get_mass();
00358     real mass_of_node = first_child->get_mass() + second_child->get_mass();
00359 
00360     if (!within_tolerance(mass_of_node - node->get_mass(), node->get_mass())) {
00361 
00362         node->pretty_print_node(cerr);
00363         cerr << " has incorrect total mass "
00364              << mass_of_node << " " << node->get_mass();
00365 
00366     } else {
00367 
00368         cm /= mass_of_node;
00369         if (!within_tolerance(abs(cm), abs(d))) {
00370             node->pretty_print_node(cerr);
00371             cerr << " has incorrect center of mass " << id
00372                  << ": " << cm << endl;
00373         }
00374     }
00375 }
00376 
00377 void check_consistency_of_nodes(hdyn * node)
00378 {
00379     if (node->get_parent() != NULL) {
00380 
00381         // node is not root
00382 
00383         if (node->get_oldest_daughter() != NULL) {
00384           check_consistency_of_node(node, &hdyn::get_pos, "pos");
00385           check_consistency_of_node(node, &hdyn::get_vel, "vel");
00386           check_consistency_of_node(node, &hdyn::get_pred_pos, "ppos");
00387           check_consistency_of_node(node, &hdyn::get_pred_vel, "pvel");
00388           check_consistency_of_node(node, &hdyn::get_acc, "acc");
00389           check_consistency_of_node(node, &hdyn::get_jerk, "jerk");
00390           check_consistency_of_node(node, &hdyn::get_old_acc, "old_acc");
00391           check_consistency_of_node(node, &hdyn::get_old_jerk, "old_jerk");
00392         }
00393     }
00394 
00395     if (node->get_oldest_daughter() != NULL) {
00396         for (hdyn * bb = node->get_oldest_daughter(); bb != NULL;
00397              bb = bb->get_younger_sister()) {
00398             check_consistency_of_nodes(bb);
00399         }
00400     }
00401 }
00402 
00403 void hdyn::setup_binary_node()
00404 {
00405     hdyn *younger_daughter = (hdyn *) oldest_daughter->get_younger_sister();
00406     hdyn *older_daughter = (hdyn *) oldest_daughter;
00407 
00408     // Check if the times of daughters are exactly the same.
00409 
00410     if (older_daughter->get_time() != younger_daughter->get_time()) {
00411         cerr << "setup_binary_node, times of daughters are different:"
00412              << older_daughter->get_time() << " "
00413              << younger_daughter->get_time() << "\n";
00414         exit(1);
00415     }
00416 
00417     // Set up physical quantities of the parent node.
00418     // Note that pot is not set -- must be recalculated.
00419 
00420     time = older_daughter->time;
00421     timestep = min(older_daughter->timestep, younger_daughter->timestep);
00422 
00423     mass = older_daughter->mass + younger_daughter->mass;
00424 
00425     real f1 = older_daughter->mass / mass;
00426     real f2 = younger_daughter->mass / mass;
00427     pos = f1 * older_daughter->pos + f2 * younger_daughter->pos;
00428     vel = f1 * older_daughter->vel + f2 * younger_daughter->vel;
00429     acc = f1 * older_daughter->acc + f2 * younger_daughter->acc;
00430     jerk = f1 * older_daughter->jerk + f2 * younger_daughter->jerk;
00431 
00432     store_old_force();
00433 
00434     pred_pos = f1 * older_daughter->pred_pos + f2 * younger_daughter->pred_pos;
00435     pred_vel = f1 * older_daughter->pred_vel + f2 * younger_daughter->pred_vel;
00436 
00437     // Set up the physical coordinates of the daughters.
00438 
00439     older_daughter->pos -= pos;
00440     older_daughter->vel -= vel;
00441     older_daughter->acc -= acc;
00442     older_daughter->jerk -= jerk;
00443     older_daughter->store_old_force();
00444     older_daughter->pred_pos -= pred_pos;
00445     older_daughter->pred_vel -= pred_vel;
00446 
00447     younger_daughter->pos -= pos;
00448     younger_daughter->vel -= vel;
00449     younger_daughter->acc -= acc;
00450     younger_daughter->jerk -= jerk;
00451     younger_daughter->store_old_force();
00452     younger_daughter->pred_pos -= pred_pos;
00453     younger_daughter->pred_vel -= pred_vel;
00454 }
00455 
00456 void create_binary_from_toplevel_nodes(hdyn * bi, hdyn * bj)
00457 {
00458     if (!bi->is_top_level_node())
00459         err_exit("create_binary_from_toplevel_nodes: bi not top level node");
00460     if (!bj->is_top_level_node())
00461         err_exit("create_binary_from_toplevel_nodes: bj not top level node");
00462 
00463     detach_node_from_general_tree(*bj);
00464     bj->set_younger_sister(NULL);
00465     bj->set_elder_sister(NULL);
00466 
00467     hdyn *new_n = new hdyn();
00468     insert_node_into_binary_tree(*bj, *bi, *new_n);
00469     bj->get_parent()->setup_binary_node();
00470 
00471     label_binary_node(bj->get_parent());
00472 }
00473 
00474 local vector change_of_absolute_something_of_parent(hdyn * node,
00475                                                     vector & dx,
00476                                                     real dm,
00477                                                     hdyn_VMF_ptr something)
00478 {
00479     return ((node->get_mass() + dm) * dx + (node->*something) () * dm)
00480          / (node->get_mass() + node->get_binary_sister()->get_mass() + dm);
00481 }
00482 
00483 local void adjust_relative_something_of_sister(vector
00484                                                  absolute_change_of_parent,
00485                                                real dm,
00486                                                hdyn * node,
00487                                                hdyn_MF_ptr inc_something)
00488 {
00489     hdyn *sister;
00490     sister = (hdyn *) (node->get_binary_sister());
00491 
00492     (sister->*inc_something) (-absolute_change_of_parent);
00493     real dummy = dm;            // to keep the compiler happy
00494 
00495 }
00496 
00497 local void adjust_parent_and_sister(hdyn * node,
00498                                     hdyn * ancestor,
00499                                     vector d_something,
00500                                     real dm,
00501                                     hdyn_VMF_ptr get_something,
00502                                     hdyn_MF_ptr inc_something,
00503                                     int modify_node)
00504 {
00505 //    cerr << "in adjust_parent_and_sister...\n" << flush;
00506 
00507     vector d_parent = 0;
00508 
00509     if (node != ancestor) {
00510         if(node->is_low_level_node()){
00511             d_parent = change_of_absolute_something_of_parent(node,
00512                                                               d_something,
00513                                                               dm,
00514                                                               get_something);
00515             adjust_relative_something_of_sister(d_parent, dm, node,
00516                                                 inc_something);
00517             if (node->get_parent() != ancestor)
00518                 adjust_parent_and_sister(node->get_parent(), ancestor,
00519                                          d_parent, dm,
00520                                          get_something, inc_something, 1);
00521         }
00522     } else {
00523         d_parent = 0;
00524     }
00525 
00526     if (modify_node)
00527         (node->*inc_something) (d_something - d_parent);
00528 }
00529 
00530 local void init_pred_of_parent_and_sister(hdyn * node,
00531                                           hdyn * ancestor)
00532 {
00533     if (node != ancestor) {
00534         if(node->is_low_level_node()){
00535             node->get_parent()->init_pred();
00536             node->get_binary_sister()->init_pred();
00537             if (node->get_parent() != ancestor)
00538             init_pred_of_parent_and_sister(node->get_parent(), ancestor);
00539         }
00540     }
00541     node->init_pred();
00542 }
00543 
00544 local void adjust_parent_mass(hdyn * node,
00545                               hdyn * ancestor,
00546                               real dm)
00547 {
00548     node->set_mass(node->get_mass() + dm);
00549 
00550     if (node->get_parent() != ancestor)
00551         adjust_parent_mass(node->get_parent(), ancestor, dm);
00552 }
00553 
00554 void remove_node_and_correct_upto_ancestor(hdyn * ancestor, hdyn * node)
00555 {
00556     // cerr << "in remove_node_and_correct_upto_ancestor" << endl;
00557 
00558     if (node->is_top_level_node()) {
00559 
00560         detach_node_from_general_tree(*node);
00561 
00562     } else {
00563 
00564         real dm = -node->get_mass();
00565         vector dz;
00566         dz = 0;
00567 
00568         hdyn *parent = node->get_parent();
00569         hdyn *sister = (hdyn *) (node->get_binary_sister());
00570 
00571         adjust_parent_and_sister(node, ancestor, dz, dm,
00572                                  &hdyn::get_pos,
00573                                  &hdyn::inc_pos, 0);
00574         adjust_parent_and_sister(node, ancestor, dz, dm,
00575                                  &hdyn::get_vel,
00576                                  &hdyn::inc_vel, 0);
00577         adjust_parent_and_sister(node, ancestor, dz, dm,
00578                                  &hdyn::get_acc,
00579                                  &hdyn::inc_acc, 0);
00580         adjust_parent_and_sister(node, ancestor, dz, dm,
00581                                  &hdyn::get_old_acc,
00582                                  &hdyn::inc_old_acc, 0);
00583         adjust_parent_and_sister(node, ancestor, dz, dm,
00584                                  &hdyn::get_jerk,
00585                                  &hdyn::inc_jerk, 0);
00586         adjust_parent_and_sister(node, ancestor, dz, dm,
00587                                  &hdyn::get_old_jerk,
00588                                  &hdyn::inc_old_jerk, 0);
00589 
00590         adjust_parent_mass(node, ancestor, dm);
00591         init_pred_of_parent_and_sister(node, ancestor);
00592 
00593         sister->set_pos(parent->get_pos());
00594         sister->set_vel(parent->get_vel());
00595         sister->set_acc(parent->get_acc());
00596         sister->set_jerk(parent->get_jerk());
00597         sister->store_old_force();
00598 
00599         // Sister pred_xxx may be needed before the next step.
00600         // Not currently correct -- fix that here.
00601 
00602         sister->init_pred();
00603 
00604         detach_node_from_binary_tree(*node);
00605 
00606         // parent must be deleted.... (17-Aug-1996)
00607 
00608         delete parent;
00609     }
00610 }
00611 
00612 vector something_relative_to_ancestor(hdyn * bj,
00613                                       hdyn * bi,
00614                                       hdyn_VMF_ptr get_something)
00615 {
00616 #ifdef DEBUG
00617     cerr << "something_relative_to_ancestor\n";
00618 #endif
00619     vector d_something = 0;
00620     for (hdyn * b = bi; b != bj; b = b->get_parent()) {
00621         if (b == NULL) {
00622             cerr << "something_relative_to_ancestor: Error, bj is not the "
00623                  << "ancestor of bi\n";
00624             exit(1);
00625         }
00626         d_something += (b->*get_something)();
00627     }
00628     return d_something;
00629 }
00630 
00631 vector hdyn_something_relative_to_root(hdyn * bi,
00632                                        hdyn_VMF_ptr get_something)
00633 {
00634 #ifdef DEBUG
00635     cerr << "hdyn_something_relative_to_root\n";
00636 #endif
00637 
00638     // Special treatment of top-level nodes:
00639 
00640     if (bi->is_top_level_node())
00641 
00642         return (bi->*get_something) ();
00643 
00644     else {
00645 
00646         vector d_something = 0;
00647         for (hdyn * b = bi; b->get_parent() != NULL; b = b->get_parent()) {
00648             if (b == NULL)
00649                 err_exit(
00650                 "hdyn_something_relative_to_root: bj not the ancestor of bi");
00651             d_something += (b->*get_something)();
00652         }
00653         return d_something;
00654     }
00655 }
00656 
00657 local vector relative_something(hdyn * ancestor,
00658                                 hdyn * bj,
00659                                 hdyn * bi,
00660                                 hdyn_VMF_ptr get_something)
00661 {
00662 #ifdef DEBUG
00663     cerr << "relative_something\n";
00664 #endif
00665     return (something_relative_to_ancestor(ancestor, bj, get_something)
00666             - something_relative_to_ancestor(ancestor, bi, get_something));
00667 }
00668 
00669 local void calculate_new_physical_quantities(hdyn * ancestor,
00670                                              hdyn * node,
00671                                              hdyn * new_sister)
00672 {
00673     node->set_pos(relative_something(ancestor, node, new_sister,
00674                                      &hdyn::get_pos));
00675     node->set_vel(relative_something(ancestor, node, new_sister,
00676                                      &hdyn::get_vel));
00677     node->set_acc(relative_something(ancestor, node, new_sister,
00678                                      &hdyn::get_acc));
00679     node->set_old_acc(relative_something(ancestor, node, new_sister,
00680                                          &hdyn::get_old_acc));
00681     node->set_jerk(relative_something(ancestor, node, new_sister,
00682                                       &hdyn::get_jerk));
00683     node->set_old_jerk(relative_something(ancestor, node, new_sister,
00684                                           &hdyn::get_old_jerk));
00685 }
00686 
00687 local void insert_node_and_correct_upto_ancestor(hdyn * ancestor,
00688                                                  hdyn * node,
00689                                                  hdyn * new_sister)
00690 {
00691     if (new_sister->get_parent() == NULL) {
00692         add_node(*node, *new_sister);   // Convention for new top-level node
00693 
00694     } else {
00695         hdyn *new_n = new hdyn();
00696         insert_node_into_binary_tree(*node, *new_sister, *new_n);
00697 
00698         hdyn *parent = node->get_parent();
00699         parent->set_mass(new_sister->get_mass());
00700         parent->set_pos(new_sister->get_pos());
00701         parent->set_vel(new_sister->get_vel());
00702         parent->set_acc(new_sister->get_acc());
00703         parent->set_jerk(new_sister->get_jerk());
00704         parent->store_old_force();
00705         new_sister->set_pos(0);
00706         new_sister->set_vel(0);
00707         new_sister->set_acc(0);
00708         new_sister->set_jerk(0);
00709         new_sister->set_old_acc(0);
00710         new_sister->set_old_jerk(0);
00711 
00712         // Attach the node with zero mass (so no corrections needed),
00713         // then correct the mass to the proper value.
00714 
00715         real dm = node->get_mass();
00716         node->set_mass(0);
00717         static vector dx, dv, da, dj;
00718         dx = dv = da = dj = 0;
00719 
00720         adjust_parent_and_sister(node, ancestor, dx, dm,
00721                                  &hdyn::get_pos,
00722                                  &hdyn::inc_pos, 1);
00723         adjust_parent_and_sister(node, ancestor, dv, dm,
00724                                  &hdyn::get_vel,
00725                                  &hdyn::inc_vel, 1);
00726         adjust_parent_and_sister(node, ancestor, da, dm,
00727                                  &hdyn::get_acc,
00728                                  &hdyn::inc_acc, 1);
00729         adjust_parent_and_sister(node, ancestor, da, dm,
00730                                  &hdyn::get_old_acc,
00731                                  &hdyn::inc_old_acc, 1);
00732         adjust_parent_and_sister(node, ancestor, dj, dm,
00733                                  &hdyn::get_jerk,
00734                                  &hdyn::inc_jerk, 1);
00735         adjust_parent_and_sister(node, ancestor, dj, dm,
00736                                  &hdyn::get_old_jerk,
00737                                  &hdyn::inc_old_jerk, 1);
00738         adjust_parent_mass(node, ancestor, dm);
00739         init_pred_of_parent_and_sister(node, ancestor);
00740 
00741         label_binary_node(ancestor);
00742     }
00743 }
00744 
00745 void correct_leaf_for_change_of_mass(hdyn * node, real dm)
00746 {
00747     hdyn *parent = node->get_parent();
00748     hdyn *ancestor = node->get_root();
00749     static vector dx, dv, da, dj;
00750     dx = dv = da = dj = 0;
00751     if(node->is_low_level_node()){
00752         adjust_parent_and_sister(node, ancestor, dx, dm,
00753                                  &hdyn::get_pos,
00754                                  &hdyn::inc_pos, 1);
00755         adjust_parent_and_sister(node, ancestor, dv, dm,
00756                                  &hdyn::get_vel,
00757                                  &hdyn::inc_vel, 1);
00758         adjust_parent_and_sister(node, ancestor, da, dm,
00759                                  &hdyn::get_acc,
00760                                  &hdyn::inc_acc, 1);
00761         adjust_parent_and_sister(node, ancestor, da, dm,
00762                                  &hdyn::get_old_acc,
00763                                  &hdyn::inc_old_acc, 1);
00764         adjust_parent_and_sister(node, ancestor, dj, dm,
00765                                  &hdyn::get_jerk,
00766                                  &hdyn::inc_jerk, 1);
00767         adjust_parent_and_sister(node, ancestor, dj, dm,
00768                                  &hdyn::get_old_jerk,
00769                                  &hdyn::inc_old_jerk, 1);
00770         adjust_parent_mass(node, ancestor, dm);
00771         init_pred_of_parent_and_sister(node, ancestor);
00772     }else{
00773         node->set_mass(node->get_mass() + dm);
00774     }
00775 }
00776 
00777 void correct_leaf_for_change_of_vector(hdyn * node,
00778                                        vector d_something,
00779                                        hdyn_VMF_ptr get_something,
00780                                        hdyn_MF_ptr inc_something)
00781 {
00782     hdyn *ancestor = node->get_root();
00783     adjust_parent_and_sister(node, ancestor, d_something, 0,
00784                              get_something, inc_something, 1);
00785     init_pred_of_parent_and_sister(node, ancestor);
00786 }
00787 
00788 
00789 // move_node: reattach a node as a sister of another node.
00790 
00791 void move_node(hdyn * node_to_move,
00792                hdyn * place_to_insert)
00793 {
00794     // cerr << "in move_node" << endl << flush;
00795 
00796     // Make a copy of node_to_move prior to any tree operations.
00797     // Code here is all very murky, mainly to avoid compiler problems...
00798 
00799     static hdyn tmp;            // NECESSARY to avoid an occasional unexplained
00800                                 // Bus Error on return with g++ version 2.2.2!!
00801 
00802     tmp = *node_to_move;        // NECESSARY to ensure proper initialization!
00803 
00804     hdyn *ancestor = (hdyn *) common_ancestor(node_to_move, place_to_insert);
00805 
00806     calculate_new_physical_quantities(ancestor, &tmp, place_to_insert);
00807 
00808     // The parent of node_to_move will be deleted by function
00809     // remove_node_and_correct_upto_ancestor.  Correct any nn
00810     // references, at least within the local clump, before the
00811     // tree structure is modified.
00812 
00813     for_all_nodes(hdyn, node_to_move->get_top_level_node(), bb) {
00814         if (bb->get_nn() == node_to_move->get_parent()) {
00815             bb->set_nn(NULL);
00816             // cerr << "set nn of " << bb->format_label() << " NULL" << endl;
00817         }
00818     }
00819 
00820     hdyn* temp_ancestor = NULL;
00821     if (node_to_move->get_parent() == ancestor) {
00822 
00823         // In this case, the ancestor will be removed, so we must reset it.
00824 
00825         temp_ancestor = (hdyn *) node_to_move->get_binary_sister();
00826     }
00827     remove_node_and_correct_upto_ancestor(ancestor, node_to_move);
00828 
00829     real temp_mass;
00830     if (temp_ancestor) {
00831         temp_mass = ancestor->get_mass();
00832         ancestor = temp_ancestor;
00833     }
00834 
00835     // Restore the original node_to_move.
00836 
00837     *node_to_move = tmp;
00838     insert_node_and_correct_upto_ancestor(ancestor, node_to_move,
00839                                           place_to_insert);
00840 
00841     if (temp_ancestor)
00842         ancestor->set_mass(temp_mass);
00843 
00844     // Bus error sometimes occurs AFTER last executable line here...
00845 
00846     // Update times and labels:
00847 
00848     node_to_move->get_parent()->set_time(place_to_insert->get_time());
00849     node_to_move->get_parent()->set_timestep(place_to_insert->get_timestep());
00850 
00851     // Clear all pointers in tmp.  (Don't attempt to follow *any*
00852     // of them, as they point into *node_to_move.)
00853 
00854     tmp.null_pointers();
00855 
00856     // cerr << "leaving move_node" << endl << flush;
00857 }
00858 
00859 
00860 // move_node_by_index: reattach node #i as a sister of node #j
00861 
00862 void move_node_by_index(int i, int j, hdyn * root)
00863 {
00864     move_node((hdyn *) node_with_index(i, root),
00865               (hdyn *) node_with_index(j, root));
00866 }
00867 
00868 //-------------------------------------------------------------------------
00869 
00870 // These are hdyn versions of functions in dyn/util/dyn_di.C -- they
00871 // perform the same operations and have the same calling sequences,
00872 // but they also set the hdyn::pot member data.
00873 
00874 local void accumulate_energies(hdyn * root, hdyn * b, real eps2,
00875                                real & epot, real & ekin, real & etot,
00876                                bool cm = false)
00877 {
00878     if (b->get_oldest_daughter()) {
00879 
00880         // Start/continue the recursion.  In the cm = true case,
00881         // expand daughters only if root really is root.
00882 
00883         if (!cm || b->is_root())
00884             for_all_daughters(hdyn, b, bb)
00885                 accumulate_energies(root, bb, eps2, epot, ekin, etot, cm);
00886 
00887         etot = epot + ekin;
00888     }
00889 
00890     // Do the work.
00891 
00892     if (!b->is_root()) {
00893         real m = b->get_mass();
00894 
00895         // pot_on_general_node is a dyn function returning the total
00896         // potential per unit mass of a top-level node, or the potential
00897         // per unit mass of a low-level node due to its binary sister
00898 
00899         real pot = m * pot_on_general_node(root, b, eps2, cm);
00900 
00901         epot += 0.5 * m * pot_on_general_node(root, b, eps2, cm);
00902         ekin += 0.5 * m * square(b->get_vel());
00903 
00904         b->set_pot(pot);        // will only be used for top-level nodes
00905     }
00906 }
00907 
00908 void calculate_energies(hdyn * root, real eps2,
00909                         real & epot, real & ekin, real & etot,
00910                         bool cm)        // default = false
00911 {
00912     epot = ekin = etot = 0;
00913     accumulate_energies(root, root, eps2, epot, ekin, etot, cm);
00914 }
00915 
00916 #endif