Main Page   Class Hierarchy   Data Structures   File List   Data Fields   Globals  

molecules.C

Go to the documentation of this file.
00001 
00005 
00006 //.............................................................................
00007 //    version 1:  July 1989   Piet Hut               email: piet@iassns.bitnet
00008 //                            Institute for Advanced Study, Princeton, NJ, USA
00009 //    version 2:  Dec 1992   Piet Hut  --  adopted to the new C++-based Starlab
00010 //.............................................................................
00011 //     NOTE: this program is an adoptation of a Newton0 diagnostics piece.
00012 //.............................................................................
00013 
00014 #include "dyn.h"
00015 
00016 #ifdef TOOLBOX
00017 
00018 #define  BUFF_LENGTH    128
00019 
00020 /*-----------------------------------------------------------------------------
00021  *  macros for the hierarchical decomposition of a few-body system:
00022  *      the following macro definitions govern the bookkeeping;
00023  *      they should migrate to their own header file in due time.
00024  *             NOTE: only one-digit star numbering implemented, i.e.
00025  *                   only valid for an N-body sub-system with N < 11
00026  *-----------------------------------------------------------------------------
00027  */
00028 #define  NBODY_MAX    10                     /* maximum number of particles  */
00029 #define  NNODE_MAX    (2 * NBODY_MAX - 1)    /* maximum number of nodes      */
00030 
00031 /*-----------------------------------------------------------------------------
00032  *  the following macros are used to access
00033  *                int  admin_table[NNODE_MAX][N_ADMIN];
00034  *  the entrees for the second variable are:
00035  *                                           number of daughter nodes   
00036  *                                           covered or visible         
00037  *                                           stablve, unstable, barely unbound
00038  *                                              or strongly unbound
00039  *                                           index of first daughter  
00040  *                                           index of second daughter 
00041  *                                           ...
00042  *                                           index of last daughter   
00043  *-----------------------------------------------------------------------------
00044  */
00045 #define  N_ADMIN      (3 + NBODY_MAX)
00046 
00047 #define  Ndaughter(ptr, offset)          ((ptr)[(offset)][0])
00048 
00049 #define  Visibility(ptr, offset)         ((ptr)[(offset)][1])
00050 #define  COVERED      1
00051 #define  VISIBLE      2
00052 
00053 #define  Stability(ptr, offset)          ((ptr)[(offset)][2])
00054 #define  STABLE            1                          /* bound, and stable   */
00055 #define  UNSTABLE          2                          /* bound, and unstable */
00056 #define  BARELY_UNBOUND    3                          /* barely unbound      */
00057 #define  STRONGLY_UNBOUND  4                          /* barely unbound      */
00058 #define  UNKNOWN           5
00059 
00060 #define  Daughter(ptr, offset, i)        ((ptr)[(offset)][3 + (i)])
00061 
00062 /*-----------------------------------------------------------------------------
00063  *  the following macros are used to access
00064  *                real  struct_table[NNODE_MAX][N_STRUCT];
00065  *  the entrees of the second variable are:
00066  *                                          radius
00067  *                                          a: semimajor axis
00068  *                                          e: eccentricity
00069  *-----------------------------------------------------------------------------
00070  */
00071 #define  N_STRUCT     3
00072 
00073 #define  Radius(ptr, offset)    ((ptr)[(offset)][0])
00074 #define  Apair(ptr, offset)     ((ptr)[(offset)][1])
00075 #define  Epair(ptr, offset)     ((ptr)[(offset)][2])
00076 
00077 local void  append_object(char hier_string[BUFF_LENGTH],
00078                           char tmp_string[BUFF_LENGTH]);
00079 local int  count_visible_nodes(int admin_table[NNODE_MAX][N_ADMIN],
00080                                    int nnode);
00081 local void  dissolve_node(dyn * b, int * nnodeptr, int i,
00082                           int admin_table[NNODE_MAX][N_ADMIN],
00083                           real struct_table[NNODE_MAX][N_STRUCT],
00084                           real dist_table[NNODE_MAX][NNODE_MAX]);
00085 local bool  dissolve_visible_strongly_unbound_nodes(dyn *b, int *nnodeptr,
00086                                         int admin_table[NNODE_MAX][N_ADMIN],
00087                                         real struct_table[NNODE_MAX][N_STRUCT],
00088                                         real dist_table[NNODE_MAX][NNODE_MAX]);
00089 local void  drop_tailnode(int * nnodeptr,
00090                               int admin_table[NNODE_MAX][N_ADMIN]);
00091 local char *end_of_string(char * the_string);
00092 local void  find_namenumbers(dyn *, char **);
00093 local bool  find_new_node(dyn * b, int * nnodeptr,
00094                           int admin_table[NNODE_MAX][N_ADMIN],
00095                           real struct_table[NNODE_MAX][N_STRUCT],
00096                           real dist_table[NNODE_MAX][NNODE_MAX], int k);
00097 local void get_binary_parameters(dyn *g1, dyn *g2, real * aptr, real * eptr);
00098 local void  init_admin_table(int admin_table[NNODE_MAX][N_ADMIN]);
00099 local void init_dist_table(dyn *b, real dist_table[NNODE_MAX][NNODE_MAX]);
00100 local void  init_struct_table(real  struct_table[NNODE_MAX][N_STRUCT]);
00101 local void  init_tuple(int tuple[NBODY_MAX],
00102                      int admin_table[NNODE_MAX][N_ADMIN], int n);
00103 local void  install_com_dynamics(dyn * b, int nnode,
00104                                int tuple[NBODY_MAX], int k);
00105 local void  install_node(int tuple[NBODY_MAX], dyn * b, int * nnodeptr,
00106                        int admin_table[NNODE_MAX][N_ADMIN],
00107                        real struct_table[NNODE_MAX][N_STRUCT],
00108                        real dist_table[NNODE_MAX][NNODE_MAX],
00109                        real radius, int k);
00110 local bool is_approximately_isolated_tuple(int tuple[NBODY_MAX], int nnode,
00111                                          int admin_table[NNODE_MAX][N_ADMIN],
00112                                          real dist_table[NNODE_MAX][NNODE_MAX],
00113                                          int k);
00114 local bool  is_isolated_tuple(int tuple[NBODY_MAX], dyn * b, int nnode,
00115                               int admin_table[NNODE_MAX][N_ADMIN],
00116                               real struct_table[NNODE_MAX][N_STRUCT],
00117                               real * radiusptr, int k);
00118 local bool  is_new_node(int tuple[NBODY_MAX], dyn * b, int * nnodeptr,
00119                         int admin_table[NNODE_MAX][N_ADMIN],
00120                         real struct_table[NNODE_MAX][N_STRUCT],
00121                         real dist_table[NNODE_MAX][NNODE_MAX], int k);
00122 local void  iterate_dissolve_visible_strongly_unbound_nodes(dyn * b,
00123                                                             int * nnodeptr,
00124                                         int admin_table[NNODE_MAX][N_ADMIN],
00125                                         real struct_table[NNODE_MAX][N_STRUCT],
00126                                         real dist_table[NNODE_MAX][NNODE_MAX]);
00127 local void  iterate_find_new_node(dyn * b, int * nnodeptr,
00128                                   int admin_table[NNODE_MAX][N_ADMIN],
00129                                   real struct_table[NNODE_MAX][N_STRUCT],
00130                                   real dist_table[NNODE_MAX][NNODE_MAX],int k);
00131 local void  map_to_integers(char hier_string[BUFF_LENGTH],
00132                             int unordered[BUFF_LENGTH], int * nptr);
00133 local int  next_element(int element, int tuple[NBODY_MAX], int n);
00134 local bool  next_tuple(int tuple[NBODY_MAX], int k, int n,
00135                      int ordered_tuple[NBODY_MAX]);
00136 local int  old_ranking(int old_array[BUFF_LENGTH],
00137                            int new_array[BUFF_LENGTH], int i, int n);
00138 local void  print_binary_parameters(char node_report[BUFF_LENGTH],
00139                                     real struct_table[NNODE_MAX][N_STRUCT],
00140                                     int member);
00141 local void  print_group_internal_structure(int nnode,
00142                                         int admin_table[NNODE_MAX][N_ADMIN],
00143                                         real struct_table[NNODE_MAX][N_STRUCT],
00144                                         char * name_table[NNODE_MAX]);
00145 local void  print_node(char node_report[BUFF_LENGTH],
00146                        int admin_table[NNODE_MAX][N_ADMIN],
00147                        real struct_table[NNODE_MAX][N_STRUCT],
00148                        char * name_table[NNODE_MAX],
00149                        int member);
00150 local void  print_radius(char node_report[BUFF_LENGTH],
00151                        real struct_table[NNODE_MAX][N_STRUCT], int member);
00152 local void  select_object(char tmp_string[BUFF_LENGTH],
00153                           char hier_string[BUFF_LENGTH], int member);
00154 local void  sort_int_array(int old_array[BUFF_LENGTH],
00155                            int new_array[BUFF_LENGTH], int n);
00156 local void  swap_body_nodes(dyn *, dyn *);
00157 local void  swap_nodes(dyn * b, int nnode, int i1, int i2,
00158                      int admin_table[NNODE_MAX][N_ADMIN],
00159                      real struct_table[NNODE_MAX][N_STRUCT],
00160                      real dist_table[NNODE_MAX][NNODE_MAX]);
00161 local real  tuple_kin_energy(int tuple[NBODY_MAX], dyn * b, int k);
00162 local real  tuple_pot_energy(int tuple[NBODY_MAX], dyn * b,
00163                            real dist_table[NNODE_MAX][NNODE_MAX], int k);
00164 local real  tuple_radius(int tuple[NBODY_MAX], dyn * b,
00165                          real struct_table[NNODE_MAX][N_STRUCT],
00166                          vector & com_pos, int k);   
00167 local void  unwrap_string(char substring[BUFF_LENGTH],
00168                           char * head_ptr, char * tail_ptr);
00169 local void  update_admin_table(int tuple[NBODY_MAX], dyn * b, int nnode,
00170                              int admin_table[NNODE_MAX][N_ADMIN],
00171                              real struct_table[NNODE_MAX][N_STRUCT],
00172                              real dist_table[NNODE_MAX][NNODE_MAX], int k);
00173 local void  update_dist_table(dyn * b, int nnode,
00174                             real dist_table[NNODE_MAX][NNODE_MAX]);
00175 local void  update_struct_table(int nnode,
00176                               real struct_table[NNODE_MAX][N_STRUCT],
00177                               real radius, real a, real e, int k);
00178 local void  wrap_string(char substring[BUFF_LENGTH],
00179                         char head_char, char tail_char);
00180 
00181 /*-----------------------------------------------------------------------------
00182  *  find_group_internal_structure  --  
00183  *-----------------------------------------------------------------------------
00184  */
00185 void  find_group_internal_structure(dyn * b
00186 //                                  , bool s_flag
00187                                     )
00188 //bool  s_flag;                      /* if TRUE: simple, streamlined output  */
00189     {
00190 //    int  pid;                           /* group id */
00191     int  nnode;
00192     int *nnodeptr;
00193     int  admin_table[NNODE_MAX][N_ADMIN];
00194     char *name_table[NNODE_MAX];
00195     real  struct_table[NNODE_MAX][N_STRUCT];
00196     real  dist_table[NNODE_MAX][NNODE_MAX];
00197     dyn *bi;
00198 
00199     for (nnode = 0, bi = b->get_oldest_daughter(); bi != NULL;
00200          bi = bi->get_younger_sister())
00201         nnode++;    
00202 
00203     if (nnode > NBODY_MAX)
00204         {
00205         cerr << "find_group_internal_structure: n = " << nnode
00206              << " > NBODY_MAX = " << NBODY_MAX << "\n";
00207         exit(1);
00208         }
00209 
00210     find_namenumbers(b, name_table);
00211 
00212     init_dist_table(b, dist_table);
00213     init_admin_table(admin_table);
00214     init_struct_table(struct_table);
00215 
00216     nnodeptr = &nnode;
00217 
00218     iterate_find_new_node(b, nnodeptr, admin_table, struct_table, dist_table,
00219                           2);
00220 
00221     iterate_dissolve_visible_strongly_unbound_nodes(b, nnodeptr, admin_table,
00222                                                     struct_table, dist_table);
00223 
00224 //    if (!s_flag && findiq(Pdata(Pdown(pn)), "group", &pid))
00225 //        printf("group #%d\n", pid);
00226 
00227     print_group_internal_structure(nnode, admin_table, struct_table,
00228                                    name_table);
00229     for (int i = 0; i < NNODE_MAX; i++)
00230         if (name_table[i])
00231             delete name_table[i];
00232     }
00233 
00234 /*-----------------------------------------------------------------------------
00235  *  find_namenumbers  --  
00236  *                        NOTE: implemented temporarily through index numbers
00237  *-----------------------------------------------------------------------------
00238  */
00239 local void  find_namenumbers(dyn * b, char * name_table[NNODE_MAX])
00240     {
00241     int  i;
00242     dyn * bi;
00243 
00244     for (i = 0, bi = b->get_oldest_daughter(); bi != NULL;
00245          i++, bi = bi->get_younger_sister())
00246         {
00247         if (bi->get_index() >= 0)
00248             {
00249             name_table[i] = new char[BUFF_LENGTH];
00250             sprintf(name_table[i], "%d", bi->get_index());
00251             }
00252         else
00253             name_table[i] = NULL;
00254         }
00255     }
00256 
00257 /*-----------------------------------------------------------------------------
00258  *  iterate_find_new_node  --  
00259  *-----------------------------------------------------------------------------
00260  */
00261 local void  iterate_find_new_node(dyn * b, int * nnodeptr,
00262                                   int admin_table[NNODE_MAX][N_ADMIN],
00263                                   real struct_table[NNODE_MAX][N_STRUCT],
00264                                   real dist_table[NNODE_MAX][NNODE_MAX], int k)
00265     {
00266     int  n_visible_nodes;
00267 
00268     n_visible_nodes = count_visible_nodes(admin_table, *nnodeptr);
00269     if (k > n_visible_nodes)
00270         return;
00271 
00272     if (find_new_node(b, nnodeptr, admin_table, struct_table, dist_table,
00273                       k) == TRUE)
00274         iterate_find_new_node(b, nnodeptr, admin_table, struct_table,
00275                               dist_table, 2);
00276     else
00277         iterate_find_new_node(b, nnodeptr, admin_table, struct_table,
00278                               dist_table, k+1);
00279     }
00280 
00281 /*-----------------------------------------------------------------------------
00282  *  find_new_node  --  
00283  *-----------------------------------------------------------------------------
00284  */
00285 local bool  find_new_node(dyn * b, int * nnodeptr,
00286                           int admin_table[NNODE_MAX][N_ADMIN],
00287                           real struct_table[NNODE_MAX][N_STRUCT],
00288                           real dist_table[NNODE_MAX][NNODE_MAX], int k)
00289     {
00290     int  i;
00291     int  n;                                       /* number of visible nodes */
00292     int  tuple[NBODY_MAX];
00293     int  ordered_tuple[NBODY_MAX];
00294 
00295     n = count_visible_nodes(admin_table, *nnodeptr);
00296     if (k > n)
00297         return(FALSE);
00298 
00299     init_tuple(tuple, admin_table, n);
00300     for (i = 0; i < n; i++)
00301         ordered_tuple[i] = tuple[i];
00302 
00303     do
00304         if (is_new_node(tuple, b, nnodeptr, admin_table, struct_table,
00305                         dist_table, k) == TRUE)
00306             return(TRUE);
00307     while
00308         (next_tuple(tuple, k, n, ordered_tuple) == TRUE);
00309         
00310     return(FALSE);
00311     }
00312 
00313 /*-----------------------------------------------------------------------------
00314  *  is_new_node  --  checks wether a node is ilsolated
00315  *              note: old version:
00316  *                   if this particular tuple is not isolated, or does not form
00317  *                     a bound subsystem, then the value  FALSE  is returned;
00318  *                   if the tuple is isolated and bound,  TRUE  is returned
00319  *                     after the following actions have been carried out:
00320  *                        1) the members of this tuple are covered,
00321  *                        2) a new visible node is assigned to this tuple,
00322  *                        3) both tables are updated.
00323  *-----------------------------------------------------------------------------
00324  */
00325 local bool  is_new_node(int tuple[NBODY_MAX], dyn * b, int * nnodeptr,
00326                         int admin_table[NNODE_MAX][N_ADMIN],
00327                         real struct_table[NNODE_MAX][N_STRUCT],
00328                         real dist_table[NNODE_MAX][NNODE_MAX], int k)
00329     {
00330     real  radius;
00331 
00332     if (k < 2)
00333         {
00334         cerr << "is_new_node: k = " << k << " < 2\n";
00335         exit(1);
00336         }
00337 
00338     if (is_approximately_isolated_tuple(tuple, *nnodeptr, admin_table,
00339                                         dist_table, k) == FALSE)
00340         return(FALSE);
00341 
00342     if (is_isolated_tuple(tuple, b, *nnodeptr, admin_table, struct_table,
00343                           &radius, k) == FALSE)
00344         return(FALSE);
00345 
00346     install_node(tuple, b, nnodeptr, admin_table, struct_table, dist_table,
00347                  radius, k);
00348 
00349     return(TRUE);
00350     }
00351 
00352 /*-----------------------------------------------------------------------------
00353  *  is_approximately_isolated_tuple  --  a tuple is defined to be approximately
00354  *                                       isolated if the minimum distance
00355  *                                       between a member of the tuple and a
00356  *                                       non-member exceeds the maximum
00357  *                                       distance between a pair of members.
00358  *-----------------------------------------------------------------------------
00359  */
00360 local bool  is_approximately_isolated_tuple(int tuple[NBODY_MAX], int nnode,
00361                                          int admin_table[NNODE_MAX][N_ADMIN],
00362                                          real dist_table[NNODE_MAX][NNODE_MAX],
00363                                          int k)
00364     {
00365     int  i, j, jt;
00366     real  max_internal_pair_distance;
00367     int  visible_non_tuple[NBODY_MAX];
00368     int  n_visible_non_tuple;
00369 
00370     if (k < 2)
00371         {
00372         cerr << "is_approximately_isolated_tuple: k = " << k << " < 2\n";
00373         exit(1);
00374         }
00375 
00376     max_internal_pair_distance = 0.0;
00377 
00378     for (i = 0; i < k - 1; i++)
00379         for (j = i + 1; j < k; j++)
00380             if (dist_table[tuple[i]][tuple[j]] > max_internal_pair_distance)
00381                 max_internal_pair_distance = dist_table[tuple[i]][tuple[j]];
00382 
00383     jt = 0;
00384     j = 0;
00385     for (i = 0; i < nnode; i++)
00386         if (Visibility(admin_table, i) == VISIBLE)
00387             if (jt >= k)                            /* i.e.:  i > tuple[k-1] */
00388                 visible_non_tuple[j++] = i;
00389             else if (i < tuple[jt])
00390                 visible_non_tuple[j++] = i;
00391             else if (i == tuple[jt])
00392                 jt++;
00393             else
00394                 {
00395                 cerr << 
00396                       "is_approximately_isolated_tuple: non-visible tuple?\n";
00397                 exit(1);
00398                 }
00399 
00400     n_visible_non_tuple = j;
00401     
00402     for (i = 0; i < n_visible_non_tuple; i++)
00403         for (j = 0; j < k; j++)
00404             if (dist_table[visible_non_tuple[i]][tuple[j]]
00405                 < max_internal_pair_distance)
00406                 return(FALSE);
00407 
00408     return(TRUE);    
00409     }
00410 
00411 /*-----------------------------------------------------------------------------
00412  *  is_isolated_tuple  -- two objects, with radius Ri and Rj, are considered
00413  *                        isolated, if each point inside a sphere belonging
00414  *                        to an object is closer to any point in its own sphere
00415  *                        than to any point in the other sphere. This implies
00416  *                        that the distance between the sphere surfaces should
00417  *                        be at least as large as the largest of the two 
00418  *                        diameters, i.e. 2.0*max(Ri, Rj), and the distance
00419  *                        between the two centers of mass of the objects should
00420  *                        be at least  
00421  *                                   Ri + Rj + 2.0*max(Ri, Rj)
00422  *-----------------------------------------------------------------------------
00423  */
00424 local bool  is_isolated_tuple(int tuple[NBODY_MAX], dyn * b, int nnode,
00425                               int admin_table[NNODE_MAX][N_ADMIN],
00426                               real struct_table[NNODE_MAX][N_STRUCT],
00427                               real * radiusptr, int k)
00428     {
00429     int  i, j, jt;
00430     real  radius;
00431     real  radius_i;
00432     vector  com_pos;
00433     real  distance_to_com;
00434     int  visible_non_tuple[NBODY_MAX];
00435     int  n_visible_non_tuple;
00436     int  pi_offset;
00437     dyn * bi;
00438 
00439     radius = tuple_radius(tuple, b, struct_table, com_pos, k);
00440     *radiusptr = radius;
00441 
00442     jt = 0;
00443     j = 0;
00444 
00445     for (i = 0; i < nnode; i++)
00446         if (Visibility(admin_table, i) == VISIBLE)
00447             if (jt >= k)                            /* i.e.:  i > tuple[k-1] */
00448                 visible_non_tuple[j++] = i;
00449             else if (i < tuple[jt])
00450                 visible_non_tuple[j++] = i;
00451             else if (i == tuple[jt])
00452                 jt++;
00453             else
00454                 {
00455                 cerr << "is_isolated_tuple: non-visible tuple?\n";
00456                 exit(1);
00457                 }
00458 
00459     n_visible_non_tuple = j;
00460     for (i = 0; i < n_visible_non_tuple; i++)
00461         {
00462         bi = b->get_oldest_daughter();
00463         pi_offset = visible_non_tuple[i];
00464         while (pi_offset--)
00465             bi = bi->get_younger_sister();
00466         distance_to_com = abs(bi->get_pos() - com_pos);
00467         radius_i = Radius(struct_table, visible_non_tuple[i]);
00468         if (distance_to_com < radius + radius_i + 2.0 * max(radius, radius_i))
00469             return(FALSE);
00470         }
00471 
00472     return(TRUE);
00473     }
00474 
00475 /*-----------------------------------------------------------------------------
00476  *  how_bound_tuple  --  
00477  *-----------------------------------------------------------------------------
00478  */
00479 local void  how_bound_tuple(int tuple[NBODY_MAX], dyn * b,
00480                             real dist_table[NNODE_MAX][NNODE_MAX],
00481                             int k, bool * isbound, bool * isbarelyunbound)
00482     {
00483     real  kinetic_energy;
00484     real  abs_potential_energy;
00485 
00486     kinetic_energy = tuple_kin_energy(tuple, b, k);
00487     abs_potential_energy = -tuple_pot_energy(tuple, b, dist_table, k);
00488 
00489     *isbound = *isbarelyunbound = FALSE;
00490 
00491     if (kinetic_energy < abs_potential_energy)
00492         *isbound = TRUE;
00493     else if (kinetic_energy < 1.5 * abs_potential_energy)
00494         *isbarelyunbound = TRUE;
00495     }
00496 
00497 /*-----------------------------------------------------------------------------
00498  *  tuple_pot_energy  --  return the total potential energy of all the nodes in
00499  *                        the tuple, due to their mutual gravitational
00500  *                        attraction.
00501  *-----------------------------------------------------------------------------
00502  */
00503 local real  tuple_pot_energy(int tuple[NBODY_MAX], dyn * b,
00504                              real dist_table[NNODE_MAX][NNODE_MAX], int k)
00505     {
00506     int  i, j;
00507     real  pot_energy;
00508     int  pi_offset;
00509     int  pj_offset;
00510     dyn * bi;
00511     dyn * bj;
00512 
00513     pot_energy = 0.0;
00514 
00515     for (i = 0; i < k - 1; i++)
00516         {
00517         bi = b->get_oldest_daughter();
00518         pi_offset = tuple[i];
00519         while (pi_offset--)
00520             bi = bi->get_younger_sister();
00521 
00522         for (j = i + 1; j < k; j++)
00523             {
00524             bj = b->get_oldest_daughter();
00525             pj_offset = tuple[j];
00526             while (pj_offset--)
00527                 bj = bj->get_younger_sister();
00528 
00529             pot_energy -= ( bi->get_mass() * bj->get_mass() )
00530                           / dist_table[ tuple[i] ][ tuple[j] ];
00531             }
00532         }
00533 
00534     return(pot_energy);
00535     }
00536 
00537 /*-----------------------------------------------------------------------------
00538  *  tuple_kin_energy  --  return the total kinetic energy of all the nodes in
00539  *                        the tuple, in the center of mass frame of the tuple.
00540  *-----------------------------------------------------------------------------
00541  */
00542 local real  tuple_kin_energy(int tuple[NBODY_MAX], dyn * b, int k)
00543     {
00544     int  i;
00545     real  velocity_squared;
00546     real  kin_energy;
00547     real  total_mass;
00548     vector  mass_times_vel;
00549     vector  rel_vel;                  /* velocity with respect to c.o.m. */
00550     vector  com_vel;
00551     int  pi_offset;
00552     dyn * bi;
00553 
00554     total_mass = 0;
00555 
00556     com_vel = 0;
00557     for (i = 0; i < k; i++)
00558         {
00559         bi = b->get_oldest_daughter();
00560         pi_offset = tuple[i];
00561         while (pi_offset--)
00562             bi = bi->get_younger_sister();
00563 
00564         total_mass += bi->get_mass();
00565         mass_times_vel = bi->get_mass() * bi->get_vel();
00566         com_vel += mass_times_vel;
00567         }
00568     com_vel /= total_mass;                       /* the real c.o.m. velocity */
00569 
00570     kin_energy = 0;
00571 
00572     for (i = 0; i < k; i++)
00573         {
00574         bi = b->get_oldest_daughter();
00575         pi_offset = tuple[i];
00576         while (pi_offset--)
00577             bi = bi->get_younger_sister();
00578 
00579         rel_vel = bi->get_vel() - com_vel;
00580         velocity_squared = rel_vel * rel_vel;
00581         kin_energy += bi->get_mass() * velocity_squared;     /* 2 x too high */
00582         }
00583 
00584     kin_energy *= 0.5;                 /* corrects for the factor 2 left out */
00585 
00586     return(kin_energy);
00587     }
00588 
00589 #define  LARGE_INTEGER   100000
00590 /*-----------------------------------------------------------------------------
00591  *  representative_body  --  returns the lowest of the numbers of the atoms
00592  *                           which make up the node  i .
00593  *-----------------------------------------------------------------------------
00594  */
00595 local int  representative_body(int admin_table[NNODE_MAX][N_ADMIN], int i)
00596     {
00597     int  j;
00598     int  number_of_daughters;
00599     int  daughter_j_representative;
00600     int  lowest_number;
00601 
00602     number_of_daughters = Ndaughter(admin_table, i);
00603 
00604     if (number_of_daughters != 0)
00605         {
00606         lowest_number = LARGE_INTEGER;
00607         for (j = 0; j < number_of_daughters; j++)
00608             {
00609             daughter_j_representative = representative_body(admin_table,
00610                                                     Daughter(admin_table,i,j));
00611             if (lowest_number > daughter_j_representative)
00612                 lowest_number = daughter_j_representative;
00613             }
00614         return(lowest_number);
00615         }
00616     else
00617         return(i);
00618     }
00619 
00620 /*-----------------------------------------------------------------------------
00621  *  init_dist_table  --  computes all node-node distances between their
00622  *                       centers of mass, and installs these values in the
00623  *                       distance table.
00624  *                       note: a moderate amount of inefficiency is allowed,
00625  *                             for the sake of simplicity.
00626  *                                 For example, the table is  n x n ,  rather 
00627  *                                 than n x (n-1) / 2  [symmetry is not used];
00628  *                                 and the distances, rather than the cheaper
00629  *                                 squared distances, are computed.
00630  *-----------------------------------------------------------------------------
00631  */
00632 local void  init_dist_table(dyn * b, real dist_table[NNODE_MAX][NNODE_MAX])
00633     {
00634     int  i, j;
00635     dyn *bi, *bj;
00636     real  distance;
00637 
00638     for (i = 0, bi = b->get_oldest_daughter(); bi != NULL;
00639          i++, bi = bi->get_younger_sister())
00640         dist_table[i][i] = 0.0;
00641 
00642     for (i = 0, bi = b->get_oldest_daughter(); bi != NULL;
00643                  i++, bi = bi->get_younger_sister())
00644         for (j = i+1, bj = bi->get_younger_sister(); bj != NULL;
00645                      j++, bj = bj->get_younger_sister())
00646             {
00647             distance = abs(bi->get_pos() - bj->get_pos());
00648             dist_table[i][j] = dist_table[j][i] = distance;
00649             }
00650     }
00651 
00652 /*-----------------------------------------------------------------------------
00653  *  init_admin_table  --  
00654  *-----------------------------------------------------------------------------
00655  */
00656 local void  init_admin_table(int admin_table[NNODE_MAX][N_ADMIN])
00657     {
00658     int  i;
00659 
00660     for (i = 0; i < NNODE_MAX; i++)      /* initially the following holds:   */
00661         {
00662         Ndaughter(admin_table, i) = 0;         /* no daughter nodes present  */
00663         Visibility(admin_table, i) = VISIBLE;  /* and all nodes are visible  */
00664         Stability(admin_table, i) = STABLE;    /* and stable point particles */
00665         }
00666     }
00667 
00668 /*-----------------------------------------------------------------------------
00669  *  init_struct_table  --  
00670  *-----------------------------------------------------------------------------
00671  */
00672 local void  init_struct_table(real  struct_table[NNODE_MAX][N_STRUCT])
00673     {
00674     int  i;
00675 
00676     for (i = 0; i < NNODE_MAX; i++)
00677         Radius(struct_table, i) = 0.0;
00678     }
00679 
00680 /*-----------------------------------------------------------------------------
00681  *  count_visible_nodes  --  
00682  *-----------------------------------------------------------------------------
00683  */
00684 local int  count_visible_nodes(int admin_table[NNODE_MAX][N_ADMIN], int nnode)
00685     {
00686     int  i;
00687     int  n_visible_nodes;
00688 
00689     n_visible_nodes = 0;
00690 
00691     for (i = 0; i < nnode; i++)
00692         if (Visibility(admin_table, i) == VISIBLE)
00693             n_visible_nodes++;
00694 
00695     return( n_visible_nodes );
00696     }
00697 
00698 /*-----------------------------------------------------------------------------
00699  *  init_tuple  -  find the indices of the first  n  visible nodes, and list
00700  *                 these indices at the beginning of the  tuple[]  array.
00701  *-----------------------------------------------------------------------------
00702  */
00703 local void  init_tuple(int tuple[NBODY_MAX],
00704                        int admin_table[NNODE_MAX][N_ADMIN], int n)
00705     {
00706     int  i, j;
00707 
00708     if (n > NBODY_MAX)
00709         {
00710         cerr << "init_tuple: n = " << n << " > NBODY_MAX = "
00711              << NBODY_MAX << "\n";
00712         exit(1);
00713         }
00714 
00715     i = j = 0;
00716 
00717     while(i < NNODE_MAX && j < n)
00718         {
00719         if (Visibility(admin_table, i) == VISIBLE)
00720             tuple[j++] = i;
00721         i++;
00722         }
00723     if (j != n)
00724         {
00725         cerr << "init_tuple: j = " << j << " != n = " << n << "\n";
00726         exit(1);
00727         }
00728     }
00729 
00730 /*-----------------------------------------------------------------------------
00731  *  next_tuple  -  search for the next k-tuple from among the ordered list of
00732  *                 indices of visible nodes (provided in  ordered_tuple[] );
00733  *                 if there is a next k-tuple,
00734  *                     then replace the first k entrees in tuple[] with the
00735  *                     new k-tuple, and return TRUE;
00736  *                 else
00737  *                     return FALSE.
00738  *           note: the contents of tuple[] beyond the first k entrees, i.e.
00739  *                     tuple[k], tuple[k+1], ..., tuple[NBODY_MAX -1]
00740  *                 is undetermined, and should not be accessed.
00741  *-----------------------------------------------------------------------------
00742  */
00743 local bool  next_tuple(int tuple[NBODY_MAX], int k, int n,
00744                        int ordered_tuple[NBODY_MAX])
00745     {
00746     if (k > n)
00747         {
00748         cerr << "next_tuple: k = " << k << " > n = " << n << "\n";
00749         exit(1);
00750         }
00751 
00752     if (k < 1)
00753         return(FALSE);
00754 
00755     if (tuple[k-1] < ordered_tuple[n-1])
00756         {
00757         tuple[k-1] = next_element(tuple[k-1], ordered_tuple, n);
00758         return(TRUE);
00759         }
00760     else if (next_tuple(tuple, k-1, n-1, ordered_tuple) == TRUE)
00761         {
00762         tuple[k-1] = next_element(tuple[k-2], ordered_tuple, n);
00763         return(TRUE);
00764         }
00765     else
00766         return(FALSE);
00767     }
00768 
00769 /*-----------------------------------------------------------------------------
00770  *  next_element  -  return the next element from an array tuple[]
00771  *-----------------------------------------------------------------------------
00772  */
00773 local int  next_element(int element, int tuple[NBODY_MAX], int n)
00774     {
00775     int  i;
00776 
00777     for (i = 0; i < n - 1; i++)
00778         if (element == tuple[i])
00779             return(tuple[i+1]);
00780 
00781     int shhh = 1;           // to make the SG compiler shut up
00782 
00783     if (element == tuple[n-1])
00784         {
00785         cerr << 
00786               "next_element: element provided is the last element of tuple\n";
00787         exit(1);
00788         }
00789     else if (shhh)           // to make the SG compiler shut up
00790         {
00791         cerr << "next_element: element provided is not an element of tuple\n";
00792         exit(1);
00793         }
00794     else                    // to make the SG compiler shut up
00795         return 0;           // to make the g++ compiler happy
00796 
00797     return 0;               // to make HP g++ happy
00798     }
00799 
00800 /*-----------------------------------------------------------------------------
00801  *  tuple_radius  --  defined as the largest value of
00802  *                    (distance from the center of mass of tuple to a member
00803  *                     + internal radius of that member);
00804  *                    as a side effect, the center-of-mass positon is returned
00805  *                    in the argument com_pos[].
00806  *           new note:
00807  *                    for a k_tuple with  k = 2 but unbound, same treatment as
00808  *                    for k > 2.
00809  *           old note:
00810  *                    for a k_tuple with  k > 2  the instantaneous radius
00811  *                    is computed; 
00812  *                    for a pair (k = 2), the maximum radius is computed in the
00813  *                    approximation of an unperturbed ellipse,
00814  *                    i.e. the maximum value of
00815  *                    (distance from the center of mass of the binary to the
00816  *                    apocenter of a member + internal radius of that member).
00817  *-----------------------------------------------------------------------------
00818  */
00819 local real  tuple_radius(int tuple[NBODY_MAX], dyn * b,
00820                          real struct_table[NNODE_MAX][N_STRUCT],
00821                          vector & com_pos, int k)
00822     {
00823     int  i;
00824     real  member_radius;
00825     real  max_member_radius;
00826     real  total_mass;
00827     real  lever_arm_factor;          /* in case of a binary, i.e. k = 2      */
00828     vector  mass_times_pos;
00829     real  a, e;                      /* binary parameters for the case k = 2 */
00830     int  member_offset;
00831     int  g0_offset;
00832     int  g1_offset;
00833     dyn * member;
00834     dyn * g0;
00835     dyn * g1;
00836 
00837     total_mass = 0.0;
00838     com_pos = 0;
00839     for (i = 0; i < k; i++)
00840         {
00841         member = b->get_oldest_daughter();
00842         member_offset = tuple[i];
00843         while (member_offset--)
00844             member = member->get_younger_sister();
00845 
00846         total_mass += member->get_mass();
00847         mass_times_pos = member->get_mass() * member->get_pos();
00848         com_pos += mass_times_pos;
00849         }
00850     com_pos /= total_mass;                       /* the real c.o.m. position */
00851 
00852     if (k == 2)
00853         {
00854         g0 = g1 = b->get_oldest_daughter();
00855         g0_offset = tuple[0];
00856         g1_offset = tuple[1];
00857         while (g0_offset--)
00858             g0 = g0->get_younger_sister();
00859         while (g1_offset--)
00860             g1 = g1->get_younger_sister();
00861 
00862         get_binary_parameters(g0, g1, &a, &e);
00863         }
00864 
00865     max_member_radius = 0.0;
00866 
00867     for (i = 0; i < k; i++)
00868         {
00869         if (k == 2 && a > 0.0)
00870             {
00871             if (i == 0)
00872                 lever_arm_factor = g1->get_mass() / total_mass;
00873             else
00874                 lever_arm_factor = g0->get_mass() / total_mass;
00875             member_radius = a * (1.0 + e) * lever_arm_factor;
00876             }
00877         else
00878             {
00879             member = b->get_oldest_daughter();
00880             member_offset = tuple[i];
00881             while (member_offset--)
00882                 member = member->get_younger_sister();
00883 
00884             member_radius = abs(member->get_pos() - com_pos);
00885             }
00886         member_radius += Radius(struct_table, tuple[i]);
00887         if (member_radius > max_member_radius)
00888             max_member_radius = member_radius;
00889         }
00890 
00891     return(max_member_radius);
00892     }
00893 
00894 /*-----------------------------------------------------------------------------
00895  *  install_node  --  does all the bookkeeping needed for the introduction of
00896  *                    a new node.
00897  *                    note: incrementing the node counter should be postponed
00898  *                          to the last line, since all functions called
00899  *                          presume that the new node is  nodes[*nnodeptr] .
00900  *                     note:
00901  *                          update_struct_table()  has to be invoked before
00902  *                          invoking  update_admin_table() , so that the
00903  *                          physical parameters of the new node are availalbe
00904  *                          to determine stability of the new node.
00905  *-----------------------------------------------------------------------------
00906  */
00907 local void  install_node(int tuple[NBODY_MAX], dyn * b, int * nnodeptr,
00908                          int admin_table[NNODE_MAX][N_ADMIN],
00909                          real struct_table[NNODE_MAX][N_STRUCT],
00910                          real dist_table[NNODE_MAX][NNODE_MAX],
00911                          real radius, int k)
00912     {
00913     real  a;                                   /* semimajor axis of a binary */
00914     real  e;                                   /* eccentricity of a binary   */
00915     int  g0_offset;
00916     int  g1_offset;
00917     dyn * g0;
00918     dyn * g1;
00919 
00920     if (*nnodeptr >= NNODE_MAX)
00921         {
00922         cerr << "install_node: *nnodeptr = " << *nnodeptr
00923              << " >= NNODE_MAX = " << NNODE_MAX << "\n";
00924         exit(1);
00925         }
00926 
00927     if (k == 2)
00928         { 
00929         g0 = g1 = b->get_oldest_daughter();
00930         g0_offset = tuple[0];
00931         g1_offset = tuple[1];
00932         while (g0_offset--)
00933             g0 = g0->get_younger_sister();
00934         while (g1_offset--)
00935             g1 = g1->get_younger_sister();
00936 
00937         get_binary_parameters(g0, g1, &a, &e);
00938         }
00939 
00940     install_com_dynamics(b, *nnodeptr, tuple, k);
00941 
00942     update_dist_table(b, *nnodeptr, dist_table);
00943     update_struct_table(*nnodeptr, struct_table, radius, a, e, k);
00944     update_admin_table(tuple, b, *nnodeptr, admin_table, struct_table,
00945                        dist_table, k);
00946     *nnodeptr += 1;
00947     }
00948 
00949 /*-----------------------------------------------------------------------------
00950  *  update_dist_table  --  
00951  *-----------------------------------------------------------------------------
00952  */
00953 local void  update_dist_table(dyn * b, int nnode,
00954                               real dist_table[NNODE_MAX][NNODE_MAX])
00955     {
00956     int  i;
00957     real  distance;
00958     dyn * bi;
00959     dyn * glast;
00960 
00961     dist_table[nnode][nnode] = 0.0;
00962 
00963     glast = b->get_oldest_daughter();
00964     for (i = 0; i < nnode; i++)
00965         glast = glast->get_younger_sister();
00966 
00967     for (i = 0, bi = b->get_oldest_daughter(); i < nnode;
00968          i++, bi = bi->get_younger_sister())
00969         {
00970         distance = abs(bi->get_pos() - glast->get_pos());
00971         dist_table[i][nnode] = dist_table[nnode][i] = distance;
00972         }
00973     }
00974 
00975 #define  STABLE_SEPARATION_FACTOR     2.0              /* somewhat arbitrary */
00976 
00977 /*-----------------------------------------------------------------------------
00978  *  update_admin_table  --  determines the values for the visibility and
00979  *                          stability of a new node, and enters those values
00980  *                          in the  admin_table[] .
00981  *                          Check the node for being barely or strongly 
00982  *                          unbound; if not, then bound, so check for
00983  *                          stable or unstable.
00984  *                     note:
00985  *                          the use of how_bound_tuple() may well be overkill,
00986  *                          but as long as it works, no problem for now.
00987  *    note: oldversion:
00988  *                     note:
00989  *                          a k-tuple with  k > 2  is considered to be always
00990  *                          unstable; a binary (k = 2) is stable if both
00991  *                          members are internally stable and in addition
00992  *                          the pericenter exceeds the sum of the radii of the
00993  *                          members by a factor  STABLE_SEPARATION_FACTOR .
00994  *-----------------------------------------------------------------------------
00995  */
00996 local void  update_admin_table(int tuple[NBODY_MAX], dyn * b, int nnode,
00997                                int admin_table[NNODE_MAX][N_ADMIN],
00998                                real struct_table[NNODE_MAX][N_STRUCT],
00999                                real dist_table[NNODE_MAX][NNODE_MAX], int k)
01000     {
01001     int  i;
01002     bool  isbound;
01003     bool  isbarelyunbound;
01004     real  pericenter_distance;
01005     real  sum_of_radii;
01006 
01007     Ndaughter(admin_table, nnode) = k;
01008     Visibility(admin_table, nnode) = VISIBLE;
01009 
01010     for (i = 0; i < k; i++)
01011         {
01012         Daughter(admin_table, nnode, i) = tuple[i];
01013         Visibility(admin_table, tuple[i]) = COVERED;
01014         }
01015 
01016     how_bound_tuple(tuple, b, dist_table, k, &isbound, &isbarelyunbound);
01017 
01018     if (isbarelyunbound == TRUE)
01019         Stability(admin_table, nnode) = BARELY_UNBOUND;
01020     else if (isbound == FALSE)
01021         Stability(admin_table, nnode) = STRONGLY_UNBOUND;
01022     else if (k > 2)
01023         Stability(admin_table, nnode) = UNSTABLE;
01024     else
01025         {
01026         Stability(admin_table, nnode) = STABLE;
01027         for (i = 0; i < k; i++)
01028             if (Ndaughter(admin_table, Daughter(admin_table, nnode, i)) > 0)
01029                 if (Stability(admin_table, Daughter(admin_table, nnode, i))
01030                     != STABLE)
01031                     Stability(admin_table, nnode) = UNSTABLE;
01032 
01033         pericenter_distance =
01034             Apair(struct_table, nnode) * (1.0 - Epair(struct_table, nnode));
01035         sum_of_radii = Radius(struct_table, Daughter(admin_table, nnode, 0))
01036                      + Radius(struct_table, Daughter(admin_table, nnode, 1));
01037 
01038         if (pericenter_distance < STABLE_SEPARATION_FACTOR * sum_of_radii)
01039             Stability(admin_table, nnode) = UNSTABLE;
01040         }
01041     }
01042 
01043 /*-----------------------------------------------------------------------------
01044  *  update_struct_table  --  
01045  *-----------------------------------------------------------------------------
01046  */
01047 local void  update_struct_table(int nnode,
01048                                 real struct_table[NNODE_MAX][N_STRUCT],
01049                                 real radius, real a, real e, int k)
01050     {
01051     Radius(struct_table, nnode) = radius;
01052 
01053     if (k == 2)
01054         {
01055         Apair(struct_table, nnode) = a;
01056         Epair(struct_table, nnode) = e;
01057         }
01058     }
01059 
01060 /*-----------------------------------------------------------------------------
01061  *  get_binary_parameters  --  
01062  *-----------------------------------------------------------------------------
01063  */
01064 local void get_binary_parameters(dyn *g1, dyn *g2, real * aptr, real * eptr)
01065 //realptr  aptr;             /* pointer to a, the semimajor axis of a binary */
01066 //realptr  eptr;             /* pointer to e, the eccentricity of a binary   */
01067     {
01068     real  delta_r, delta_v, m_sum;
01069     vector  r_rel, v_rel, r_out_v;
01070     real  r_out_v_squared;
01071 
01072     delta_r = abs(g1->get_pos() - g2->get_pos());
01073     delta_v = abs(g1->get_vel() - g2->get_vel());
01074     m_sum = g1->get_mass() + g2->get_mass();
01075 
01076     *aptr = 1.0 / (2.0/delta_r - delta_v*delta_v / m_sum);
01077 
01078     r_rel = g1->get_pos() - g2->get_pos();
01079     v_rel = g1->get_vel() - g2->get_vel();
01080     r_out_v = r_rel ^ v_rel;
01081     r_out_v_squared = r_out_v * r_out_v;
01082 
01083     *eptr = sqrt(max(0.0, 1.0 - (r_out_v_squared / (m_sum * *aptr))));
01084     }
01085 
01086 /*-----------------------------------------------------------------------------
01087  *  iterate_dissolve_visible_strongly_unbound_nodes  --  
01088  *-----------------------------------------------------------------------------
01089  */
01090 local void  iterate_dissolve_visible_strongly_unbound_nodes(dyn * b,
01091                                                             int * nnodeptr,
01092                                         int admin_table[NNODE_MAX][N_ADMIN],
01093                                         real struct_table[NNODE_MAX][N_STRUCT],
01094                                         real dist_table[NNODE_MAX][NNODE_MAX])
01095     {
01096     while (dissolve_visible_strongly_unbound_nodes(b, nnodeptr,admin_table,
01097                                                    struct_table, dist_table))
01098         ;
01099     }     
01100 
01101 /*-----------------------------------------------------------------------------
01102  *  dissolve_visible_strongly_unbound_nodes  --  
01103  *-----------------------------------------------------------------------------
01104  */
01105 local bool  dissolve_visible_strongly_unbound_nodes(dyn * b, int * nnodeptr,
01106                                         int admin_table[NNODE_MAX][N_ADMIN],
01107                                         real struct_table[NNODE_MAX][N_STRUCT],
01108                                         real dist_table[NNODE_MAX][NNODE_MAX])
01109     {
01110     int  i;
01111     bool  new_activity;
01112 
01113     new_activity = FALSE;
01114     for (i = 0; i < *nnodeptr; i++)
01115         if (Visibility(admin_table, i) == VISIBLE &&
01116             Stability(admin_table, i) == STRONGLY_UNBOUND)
01117             {
01118             dissolve_node(b, nnodeptr, i, admin_table, struct_table,
01119                           dist_table);
01120             new_activity = TRUE;
01121             break;                        /* only one dissolution at a time, */
01122             }                             /* to retain integrity of nnodeptr */
01123 
01124     return(new_activity);
01125     }
01126 
01127 /*-----------------------------------------------------------------------------
01128  *  dissolve_node  --  
01129  *-----------------------------------------------------------------------------
01130  */
01131 local void  dissolve_node(dyn * b, int * nnodeptr, int i,
01132                           int admin_table[NNODE_MAX][N_ADMIN],
01133                           real struct_table[NNODE_MAX][N_STRUCT],
01134                           real dist_table[NNODE_MAX][NNODE_MAX])
01135     {
01136     if (i != *nnodeptr-1)
01137         swap_nodes(b, *nnodeptr, i, *nnodeptr-1, admin_table, struct_table,
01138                    dist_table);
01139 
01140     drop_tailnode(nnodeptr, admin_table);
01141     }
01142 
01143 /*-----------------------------------------------------------------------------
01144  *  swap_nodes  --  
01145  *-----------------------------------------------------------------------------
01146  */
01147 local void  swap_nodes(dyn * b, int nnode, int i1, int i2,
01148                        int admin_table[NNODE_MAX][N_ADMIN],
01149                        real struct_table[NNODE_MAX][N_STRUCT],
01150                        real dist_table[NNODE_MAX][NNODE_MAX])
01151     {
01152     int  tmp_int;
01153     real  tmp_real;
01154     int  j;
01155     int  i;
01156     dyn * bi1;
01157     dyn * bi2;
01158 
01159     if (i1 >= nnode)
01160         {
01161         cerr << "swap_nodes: i1 = " << i1 << " > nnode = " << nnode << "\n";
01162         exit(1);
01163         }
01164     if (i2 >= nnode)
01165         {
01166         cerr << "swap_nodes: i2 = " << i2 << " > nnode = " << nnode << "\n";
01167         exit(1);
01168         }
01169 
01170     bi1 = bi2 = b->get_oldest_daughter();
01171     while (i1--)
01172         bi1 = bi1->get_younger_sister();
01173     while (i2--)
01174         bi2 = bi2->get_younger_sister();
01175 
01176     swap_body_nodes(bi1, bi2);
01177 
01178     tmp_int = Ndaughter(admin_table, i1);
01179     Ndaughter(admin_table, i1) = Ndaughter(admin_table, i2);
01180     Ndaughter(admin_table, i2) = tmp_int;
01181 
01182     tmp_int = Visibility(admin_table, i1);
01183     Visibility(admin_table, i1) = Visibility(admin_table, i2);
01184     Visibility(admin_table, i2) = tmp_int;
01185 
01186     tmp_int = Stability(admin_table, i1);
01187     Stability(admin_table, i1) = Stability(admin_table, i2);
01188     Stability(admin_table, i2) = tmp_int;
01189     
01190     j = max(Ndaughter(admin_table, i1), Ndaughter(admin_table, i2));
01191     while (j-- > 0)           /* this swaps also unused numbers, if the two  */
01192         {                     /* Ndaughter values are unequal, but who cares */
01193         tmp_int = Daughter(admin_table, i1, j);
01194         Daughter(admin_table, i1, j) = Daughter(admin_table, i2, j);
01195         Daughter(admin_table, i2, j) = tmp_int;
01196         }
01197 /*
01198  * Now switch also the pointers to i1 and i2 from the outside,
01199  * including the cases where i = i1 and i = i2 !
01200  */
01201     for (i = 0; i < nnode; i++)
01202         if (Ndaughter(admin_table, i) > 0)
01203             for (j = 0; j < Ndaughter(admin_table, i); j++)
01204                 {
01205                 if (Daughter(admin_table, i, j) == i1)
01206                     Daughter(admin_table, i, j) = i2;
01207                 else if (Daughter(admin_table, i, j) == i2)
01208                     Daughter(admin_table, i, j) = i1;
01209                 }
01210 //
01211 // NOTE: 921219: the body of the above loop used to be as follows
01212 //       (obviously different from what originally was intended).
01213 //       The question is now whether the above repair will introduce
01214 //       new bugs or not! 
01215 //
01216 //              {
01217 //              if (Daughter(admin_table, i, j) == i1)
01218 //                  Daughter(admin_table, i, j) == i2;
01219 //              else if (Daughter(admin_table, i, j) == i2)
01220 //                  Daughter(admin_table, i, j) == i1;
01221 //              }
01222 //
01223     tmp_real = Radius(struct_table, i1);
01224     Radius(struct_table, i1) = Radius(struct_table, i2);
01225     Radius(struct_table, i2) = tmp_real;
01226     
01227     tmp_real = Apair(struct_table, i1);
01228     Apair(struct_table, i1) = Apair(struct_table, i2);
01229     Apair(struct_table, i2) = tmp_real;
01230     
01231     tmp_real = Epair(struct_table, i1);
01232     Epair(struct_table, i1) = Epair(struct_table, i2);
01233     Epair(struct_table, i2) = tmp_real;
01234     
01235     for (i = 0; i < nnode; i++)
01236         {
01237         tmp_real = dist_table[i][i1];
01238         dist_table[i][i1] = dist_table[i][i2];
01239         dist_table[i][i2] = tmp_real;
01240         }
01241 
01242     for (j = 0; j < nnode; j++)
01243         {
01244         tmp_real = dist_table[i1][j];
01245         dist_table[i1][j] = dist_table[i2][j];
01246         dist_table[i2][j] = tmp_real;
01247         }
01248     }
01249 
01250 /*-----------------------------------------------------------------------------
01251  *  swap_body_nodes  --  swap two sister nodes
01252  *-----------------------------------------------------------------------------
01253  */
01254 local void  swap_body_nodes(dyn * bi, dyn * bj)
01255     {
01256     dyn *bir, *bil, *bjr, *bjl, *bb;
01257 
01258     if ((bb = bi->get_parent()) == NULL)
01259         {
01260         cerr <<  "swap_body_nodes: node #" << bi << " has no parent\n";
01261         exit(1);
01262         }
01263     
01264     if (bi->get_parent() != bj->get_parent())
01265         {
01266         cerr << "nodes #" << bi << " and #" << bj
01267              << " have different parents\n";
01268         exit(1);
01269         }
01270 
01271     if (bi == bj)
01272         return;
01273 
01274     bir = bi->get_younger_sister();
01275     bil = bi->get_elder_sister();
01276     bjr = bj->get_younger_sister();
01277     bjl = bj->get_elder_sister();
01278 
01279     if (bil == bj)
01280         {
01281         bj->set_elder_sister(bi);
01282         bi->set_younger_sister(bj);
01283 
01284         bi->set_elder_sister(bjl);
01285         if (bjl)
01286             bjl->set_younger_sister(bi);
01287         else
01288             bb->set_oldest_daughter(bi);
01289 
01290         bj->set_younger_sister(bir);
01291         if (bir)
01292             bir->set_elder_sister(bj);
01293         }
01294     else if (bjl == bi)
01295         {
01296         bi->set_elder_sister(bj);
01297         bj->set_younger_sister(bi);
01298 
01299         bj->set_elder_sister(bil);
01300         if (bil)
01301             bil->set_younger_sister(bj);
01302         else
01303             bb->set_oldest_daughter(bj);
01304 
01305         bi->set_younger_sister(bjr);
01306         if (bjr)
01307             bjr->set_elder_sister(bi);
01308         }
01309     else
01310         {
01311         bi->set_elder_sister(bjl);
01312         bi->set_younger_sister(bjr);
01313         bj->set_elder_sister(bil);
01314         bj->set_younger_sister(bir);
01315 
01316         if (bil)
01317             bil->set_younger_sister(bj);
01318         else
01319             bb->set_oldest_daughter(bj);
01320 
01321         if (bjl)
01322             bjl->set_younger_sister(bi);
01323         else
01324             bb->set_oldest_daughter(bi);
01325  
01326         if (bir)
01327             bir->set_elder_sister(bj);
01328 
01329         if (bjr)
01330             bjr->set_elder_sister(bi);
01331         }
01332     }
01333 
01334 /*-----------------------------------------------------------------------------
01335  *  drop_tailnode  --  
01336  *-----------------------------------------------------------------------------
01337  */
01338 local void  drop_tailnode(int * nnodeptr, int admin_table[NNODE_MAX][N_ADMIN])
01339     {
01340     int  j;
01341 
01342     if (Visibility(admin_table, *nnodeptr-1) != VISIBLE)
01343         {
01344         cerr << "drop_tailnode: can't drop an invisible tail\n";
01345         exit(1);
01346         }
01347 
01348     for (j = 0; j < Ndaughter(admin_table, *nnodeptr-1); j++)
01349         Visibility(admin_table, Daughter(admin_table, *nnodeptr-1, j))
01350             = VISIBLE;
01351 
01352     (*nnodeptr)--;
01353     }
01354 
01355 /*-----------------------------------------------------------------------------
01356  *  tail_insert  --  
01357  *-----------------------------------------------------------------------------
01358  */
01359 local void tail_insert(dyn * b, dyn * new_node)
01360     {
01361     dyn *bi, *b_prev;
01362 
01363     b_prev = b->get_oldest_daughter();
01364 
01365     if (b_prev == NULL)
01366         {
01367         cerr << "tail_insert: b_prev == NULL\n";
01368         exit(1);
01369         }
01370 
01371     while (bi = b_prev->get_younger_sister())
01372         b_prev = bi;
01373         
01374     b_prev->set_younger_sister(new_node);
01375     }
01376 
01377 /*-----------------------------------------------------------------------------
01378  *  install_com_dynamics  --  compute and install the values of the  mass,
01379  *                            position and velocity for a new node, with the
01380  *                            index  nnode , and having  k  daughters,
01381  *                            whose indices are contained in the first  k
01382  *                            entrees of  tuple[] .
01383  *-----------------------------------------------------------------------------
01384  */
01385 local void  install_com_dynamics(dyn * b, int nnode,
01386                                  int tuple[NBODY_MAX], int k)
01387     {
01388     int  i;
01389     real  total_mass;
01390     vector  com_pos;
01391     vector  com_vel;
01392     vector  mass_times_pos;
01393     vector  mass_times_vel;
01394     int  member_offset;
01395     dyn * member;
01396     dyn * new_node;
01397 
01398 // just to keep the compiler happy:
01399     nnode++; nnode--;  
01400 
01401     total_mass = 0.0;
01402     com_pos = com_vel = 0;
01403     for (i = 0; i < k; i++)
01404         {
01405         member = b->get_oldest_daughter();
01406         member_offset = tuple[i];
01407         while (member_offset--)
01408             member = member->get_younger_sister();
01409 
01410         total_mass += member->get_mass();
01411         mass_times_pos = member->get_mass() * member->get_pos();
01412         com_pos += mass_times_pos;
01413         mass_times_vel = member->get_mass() * member->get_vel();
01414         com_vel += mass_times_vel;
01415         }
01416     com_pos /= total_mass;                       /* the real c.o.m. position */
01417     com_vel /= total_mass;                       /* the real c.o.m. position */
01418 
01419     new_node = new dyn;
01420     new_node->set_mass(total_mass);
01421     new_node->set_pos(com_pos);
01422     new_node->set_vel(com_vel);
01423 
01424     tail_insert(b, new_node);
01425     }
01426 
01427 /*-----------------------------------------------------------------------------
01428  *  print_group_internal_structure  --  
01429  *-----------------------------------------------------------------------------
01430  */
01431 local void  print_group_internal_structure(int nnode,
01432                                         int admin_table[NNODE_MAX][N_ADMIN],
01433                                         real struct_table[NNODE_MAX][N_STRUCT],
01434                                         char * name_table[NNODE_MAX])
01435     {
01436     int  i;
01437     char  node_report[BUFF_LENGTH];
01438 
01439     for (i = 0; i < nnode; i++)
01440         if (Ndaughter(admin_table, i) > 0)
01441             {
01442             sprintf(node_report, "  ", i);
01443 
01444             print_node(node_report, admin_table, struct_table, name_table, i);
01445 
01446             if ((Ndaughter(admin_table, i) == 2))
01447                 print_binary_parameters(node_report, struct_table, i);
01448             else
01449                 print_radius(node_report, struct_table, i);
01450 
01451             printf(node_report);
01452             }
01453     }
01454 
01455 /*-----------------------------------------------------------------------------
01456  *  print_node  --  
01457  *-----------------------------------------------------------------------------
01458  */
01459 local void  print_node(char node_report[BUFF_LENGTH],
01460                        int admin_table[NNODE_MAX][N_ADMIN],
01461                        real struct_table[NNODE_MAX][N_STRUCT],
01462                        char * name_table[NNODE_MAX],
01463                        int member)
01464     {
01465     int  i;
01466 
01467     if (Ndaughter(admin_table, member) == 0)
01468         {
01469         if (name_table[member])
01470             {
01471             sprintf(end_of_string(node_report), "%s", name_table[member]);
01472 //            *(end_of_string(node_report) - 1) = NULL;        /* erases '\n\' */
01473             }
01474         }
01475     else
01476         {
01477         if (Stability(admin_table, member) == STABLE)
01478             sprintf(end_of_string(node_report), "[");
01479         else if (Stability(admin_table, member) == UNSTABLE)
01480             sprintf(end_of_string(node_report), "(");
01481         else if (Stability(admin_table, member) == BARELY_UNBOUND)
01482             sprintf(end_of_string(node_report), "<");
01483         else if (Stability(admin_table, member) == STRONGLY_UNBOUND)
01484             sprintf(end_of_string(node_report), "{");
01485         else if (Stability(admin_table, member) == UNKNOWN)
01486             cerr << "print_node: stability of member " << member
01487                 << " is UNKNOWN\n";
01488         else
01489             cerr << "print_node: invalid value "
01490                  << Stability(admin_table, member)
01491                  << " for stability of member " << member << "\n";
01492         
01493         for (i = 0; i < Ndaughter(admin_table, member); i++)
01494             {
01495             if (i > 0)
01496                 sprintf(end_of_string(node_report), ", ");
01497             print_node(node_report, admin_table, struct_table, name_table,
01498                        Daughter(admin_table, member, i));
01499             }
01500 
01501         if (Stability(admin_table, member) == STABLE)
01502             sprintf(end_of_string(node_report), "]");
01503         else if (Stability(admin_table, member) == UNSTABLE)
01504             sprintf(end_of_string(node_report), ")");
01505         else if (Stability(admin_table, member) == BARELY_UNBOUND)
01506             sprintf(end_of_string(node_report), ">");
01507         else if (Stability(admin_table, member) == STRONGLY_UNBOUND)
01508             sprintf(end_of_string(node_report), "}");
01509         else if (Stability(admin_table, member) == UNKNOWN)
01510             cerr << "print_node: stability of member " << member
01511                 << " is UNKNOWN\n";
01512         else
01513             cerr << "print_node: invalid value "
01514                  << Stability(admin_table, member)
01515                  << " for stability of member " << member << "\n";
01516         }
01517     }
01518 
01519 /*-----------------------------------------------------------------------------
01520  *  print_binary_parameters  --  
01521  *-----------------------------------------------------------------------------
01522  */
01523 local void  print_binary_parameters(char node_report[BUFF_LENGTH],
01524                                     real struct_table[NNODE_MAX][N_STRUCT],
01525                                     int member)
01526     {
01527     int  offset_in_node_report;
01528     real  semimajor_axis;
01529 
01530     offset_in_node_report = end_of_string(node_report) - node_report;
01531 /*
01532  * if the radius information starts in the 49th column, the output will be
01533  * nicely lined up with the distances.
01534  */
01535     while (offset_in_node_report < 49)
01536         node_report[offset_in_node_report++] = ' ';
01537 
01538     semimajor_axis = Apair(struct_table, member);
01539 
01540 #if 0
01541 
01542     if (semimajor_axis < 10.0 && semimajor_axis > 0.0)
01543         sprintf(node_report + offset_in_node_report,"a = %.16f  ;  e = %.16f\n",
01544                 semimajor_axis, Epair(struct_table, member));
01545     else if (semimajor_axis < 100.0 && semimajor_axis > -10.0)
01546         sprintf(node_report + offset_in_node_report,"a = %.15f  ;  e = %.16f\n",
01547                 semimajor_axis, Epair(struct_table, member));
01548     else if (semimajor_axis < 1000.0 && semimajor_axis > -100.0)
01549         sprintf(node_report + offset_in_node_report,"a = %.14f  ;  e = %.16f\n",
01550                 semimajor_axis, Epair(struct_table, member));
01551     else if (semimajor_axis < 10000.0 && semimajor_axis > -1000.0)
01552         sprintf(node_report + offset_in_node_report,"a = %.13f  ;  e = %.16f\n",
01553                 semimajor_axis, Epair(struct_table, member));
01554     else if (semimajor_axis < 100000.0 && semimajor_axis > -10000.0)
01555         sprintf(node_report + offset_in_node_report,"a = %.12f  ;  e = %.16f\n",
01556                 semimajor_axis, Epair(struct_table, member));
01557     else if (semimajor_axis < 1000000.0 && semimajor_axis > -100000.0)
01558         sprintf(node_report + offset_in_node_report,"a = %.11f  ;  e = %.16f\n",
01559                 semimajor_axis, Epair(struct_table, member));
01560     else
01561         sprintf(node_report + offset_in_node_report,"a = %.10f  ;  e = %.16f\n",
01562                 semimajor_axis, Epair(struct_table, member));
01563 
01564 #else
01565 
01566 
01567     if (semimajor_axis < 10.0 && semimajor_axis > 0.0)
01568         sprintf(node_report + offset_in_node_report,"a = %.6f  ;  e = %8f\n",
01569                 semimajor_axis, Epair(struct_table, member));
01570     else if (semimajor_axis < 100.0 && semimajor_axis > -10.0)
01571         sprintf(node_report + offset_in_node_report,"a = %.5f  ;  e = %8f\n",
01572                 semimajor_axis, Epair(struct_table, member));
01573     else if (semimajor_axis < 1000.0 && semimajor_axis > -100.0)
01574         sprintf(node_report + offset_in_node_report,"a = %.4f  ;  e = %8f\n",
01575                 semimajor_axis, Epair(struct_table, member));
01576     else if (semimajor_axis < 10000.0 && semimajor_axis > -1000.0)
01577         sprintf(node_report + offset_in_node_report,"a = %.3f  ;  e = %8f\n",
01578                 semimajor_axis, Epair(struct_table, member));
01579     else if (semimajor_axis < 100000.0 && semimajor_axis > -10000.0)
01580         sprintf(node_report + offset_in_node_report,"a = %.2f  ;  e = %8f\n",
01581                 semimajor_axis, Epair(struct_table, member));
01582     else if (semimajor_axis < 1000000.0 && semimajor_axis > -100000.0)
01583         sprintf(node_report + offset_in_node_report,"a = %.1f  ;  e = %8f\n",
01584                 semimajor_axis, Epair(struct_table, member));
01585     else
01586         sprintf(node_report + offset_in_node_report,"a = %.0f  ;  e = %8f\n",
01587                 semimajor_axis, Epair(struct_table, member));
01588 
01589 #endif
01590     }
01591 
01592 /*-----------------------------------------------------------------------------
01593  *  print_radius  --  
01594  * NOTE: CLEANUP THE SILLY OUTPUT SWITCHES WITH A NEW FIXED_WIDTH %f PROCEDURE
01595  *-----------------------------------------------------------------------------
01596  */
01597 local void  print_radius(char node_report[BUFF_LENGTH],
01598                          real struct_table[NNODE_MAX][N_STRUCT], int member)
01599     {
01600     int  offset_in_node_report;
01601     real  r_node;
01602 
01603     offset_in_node_report = end_of_string(node_report) - node_report;
01604 /*
01605  * if the radius information starts in the 49th column, the output will be
01606  * nicely lined up with the distances.
01607  */
01608     while (offset_in_node_report < 49)
01609         node_report[offset_in_node_report++] = ' ';
01610 
01611     r_node = Radius(struct_table, member);
01612     if (r_node < 10.0)
01613         sprintf(node_report + offset_in_node_report, "R = %.6f\n", r_node);
01614     else if (r_node < 100.0)
01615         sprintf(node_report + offset_in_node_report, "R = %.5f\n", r_node);
01616     else if (r_node < 1000.0)
01617         sprintf(node_report + offset_in_node_report, "R = %.4f\n", r_node);
01618     else if (r_node < 10000.0)
01619         sprintf(node_report + offset_in_node_report, "R = %.3f\n", r_node);
01620     else if (r_node < 100000.0)
01621         sprintf(node_report + offset_in_node_report, "R = %.2f\n", r_node);
01622     else if (r_node < 1000000.0)
01623         sprintf(node_report + offset_in_node_report, "R = %.1f\n", r_node);
01624     else
01625         sprintf(node_report + offset_in_node_report, "R = %.0f\n", r_node);
01626     }
01627 
01628 /*-----------------------------------------------------------------------------
01629  *  how_much_offspring  --  computes the number of real particles below the
01630  *                          the node  member .
01631  *-----------------------------------------------------------------------------
01632  */
01633 local int  how_much_offspring(int admin_table[NNODE_MAX][N_ADMIN], int member)
01634     {
01635     int  i, k, n_total;
01636 
01637     k = Ndaughter(admin_table, member);
01638     n_total = 0;
01639     
01640     if (k > 0)
01641         for (i = 0; i < k; i++)
01642             n_total += how_much_offspring(admin_table,
01643                                           Daughter(admin_table, member, i));
01644     else
01645         n_total = 1;
01646 
01647     return(n_total);
01648     }
01649 
01650 #define  BUFF_SAFETY_MARGIN  28
01651 
01652 /*-----------------------------------------------------------------------------
01653  *  end_of_string  --  
01654  *-----------------------------------------------------------------------------
01655  */
01656 local char *end_of_string(char * the_string)
01657     {
01658     int  char_position;
01659 
01660     char_position = 0;
01661     while(the_string[char_position] != '\0')
01662         char_position++;
01663 
01664     if(char_position > BUFF_LENGTH - 1 - BUFF_SAFETY_MARGIN)
01665         {
01666         cerr << "end_of_string: too little room left in buffer\n";
01667         exit(1);
01668         }
01669 
01670     return(the_string + char_position);
01671     }
01672 
01673 /*-----------------------------------------------------------------------------
01674  *  bring_to_normal_form  --  
01675  *                   example: 3,(7,[5,[4,8]],0),6,(2,1);  ==>
01676  *                            (0,[[4,8],5],7),(1,2),3,6;
01677  *                      NOTE: only one-digit star numbering implemented, i.e.
01678  *                            only valid for an N-body sub-system with N < 11
01679  *                  OH, WELL: a much cleaner way, of course, would be to first
01680  *                            map the string to an integer array of symbols,
01681  *                            where, e.g.,  ";" --> -1 , "(" --> -2 , 
01682  *                            ")" --> -3 , "[" --> -4 , "]" --> -5 , etc,
01683  *                            and the delimiting , is simply skipped since the
01684  *                            array elements are automatically separated.
01685  *-----------------------------------------------------------------------------
01686  */
01687 local void  bring_to_normal_form(char old_hier_string[BUFF_LENGTH])
01688     {
01689     int  i, j;
01690     int  number_of_objects;
01691     int  ordered[BUFF_LENGTH];
01692     int  unordered[BUFF_LENGTH];
01693     char  new_hier_string[BUFF_LENGTH];
01694     char  substring[BUFF_LENGTH];
01695     char  head_char, tail_char;
01696 
01697     if (old_hier_string[1] == ';')
01698         return;                     /* a singleton is already in normal form */
01699 
01700     new_hier_string[0] = '\0';         /* initialize with a null_string, for */
01701     substring[0] = '\0';              /* end_of_string()  to work properly  */
01702 
01703     map_to_integers(old_hier_string, unordered, &number_of_objects);
01704 
01705     sort_int_array(unordered, ordered, number_of_objects);
01706 
01707     for (i = 0; i < number_of_objects; i++)
01708         {
01709         j = old_ranking(unordered, ordered, i, number_of_objects);
01710         select_object(substring, old_hier_string, j);
01711         unwrap_string(substring, &head_char, &tail_char);
01712         bring_to_normal_form(substring);
01713         wrap_string(substring, head_char, tail_char);
01714         append_object(new_hier_string, substring);
01715         }
01716 /*
01717  * copy the new string back unto the old string:
01718  */
01719     i = 0;    
01720     while (new_hier_string[i] != '\0' && i < BUFF_LENGTH)
01721         {
01722         old_hier_string[i] = new_hier_string[i];
01723         i++;
01724         }
01725     if (i >= BUFF_LENGTH)
01726         {
01727         cerr << "bring_to_normal_form: buffer length exceeded\n";
01728         exit(1);
01729         }
01730     if (old_hier_string[i] != '\0')
01731         {
01732         cerr << 
01733               "bring_to_normal_form: new and old hier. string diff. length\n";
01734         exit(1);
01735         }
01736     }
01737 
01738 /*-----------------------------------------------------------------------------
01739  *  map_to_integers  --  
01740  *                 NOTE: only one-digit star numbering implemented, i.e.
01741  *                       only valid for an N-body sub-system with N < 11
01742  *-----------------------------------------------------------------------------
01743  */
01744 local void  map_to_integers(char hier_string[BUFF_LENGTH],
01745                             int unordered[BUFF_LENGTH], int * nptr)
01746 //int *nptr;                           /* will count total number of objects */
01747     {
01748     int  i;
01749     int  lowest_number;
01750     int  level;
01751     char  c;
01752 
01753     *nptr = 0;
01754     level = 0;
01755     lowest_number = 10;        /* higher than highest allowed number */
01756     i = 0;
01757     while (i < BUFF_LENGTH)
01758         {
01759         c = hier_string[i++];
01760 
01761         if (c >= '0' && c <= '9')
01762             {
01763             if (c - '0' < lowest_number)
01764                 lowest_number = c - '0';
01765             }
01766         else if (c == '(' || c == '[' || c == '{' || c == '<')
01767             level++;
01768         else if (c == ')' || c == ']' || c == '}' || c == '>')
01769             level--;
01770         else if (level == 0 && (c == ',' || c == ';'))
01771             {
01772             unordered[*nptr] = lowest_number;
01773             (*nptr)++;
01774             lowest_number = 10;
01775             }
01776 
01777         if (c == ';')
01778             break;
01779         }
01780     if (i >= BUFF_LENGTH)
01781         {
01782         cerr << "map_to_integers: buffer length exceeded\n";
01783         exit(1);
01784         }
01785     }
01786 
01787 /*-----------------------------------------------------------------------------
01788  *  sort_int_array  --  
01789  *                NOTE: only one-digit star numbering implemented, i.e.
01790  *                      only valid for an N-body sub-system with N < 11
01791  *-----------------------------------------------------------------------------
01792  */
01793 local void  sort_int_array(int old_array[BUFF_LENGTH],
01794                            int new_array[BUFF_LENGTH], int n)
01795     {
01796     int  i, j;
01797     int  dummy;
01798 
01799     for (i = 0; i < n; i++)
01800         new_array[i] = old_array[i];
01801 
01802     for (i = 0; i < n - 1; i++)         /* lazy programmer's sorting method  */
01803         for (j = i+1; j < n; j++)
01804             if (new_array[j] < new_array[i])
01805                 {
01806                 dummy = new_array[i];
01807                 new_array[i] = new_array[j];
01808                 new_array[j] = dummy;
01809                 }
01810     }
01811 
01812 /*-----------------------------------------------------------------------------
01813  *  old_ranking  --  
01814  *               NOTE: only one-digit star numbering implemented, i.e.
01815  *                     only valid for an N-body sub-system with N < 11
01816  *-----------------------------------------------------------------------------
01817  */
01818 local int  old_ranking(int old_array[BUFF_LENGTH], int new_array[BUFF_LENGTH],
01819                         int i, int n)
01820     {
01821     int  j;
01822 
01823     j = 0;
01824     while (j < n)
01825         if (old_array[j++] == new_array[i])
01826             break;
01827     if (j > n)
01828         {
01829         cerr << "old_ranking: used buffer length exceeded\n";
01830         exit(1);
01831         }
01832 
01833     return(--j);
01834     }
01835 
01836 /*-----------------------------------------------------------------------------
01837  *  select_object  --  
01838  *               NOTE: only one-digit star numbering implemented, i.e.
01839  *                     only valid for an N-body sub-system with N < 11
01840  *-----------------------------------------------------------------------------
01841  */
01842 local void  select_object(char tmp_string[BUFF_LENGTH],
01843                           char hier_string[BUFF_LENGTH], int member)
01844     {
01845     int  i, j;
01846     int  level;
01847     char  c;
01848 
01849     level = 0;
01850     i = 0;
01851     while (member > 0)
01852         {
01853         c = hier_string[i];
01854         if (c == '(' || c == '[' || c == '{' || c == '<')
01855             level++;
01856         else if (c == ')' || c == ']' || c == '}' || c == '>')
01857             level--;
01858         else if (c == ',' && level == 0)
01859             member--;
01860         else if (c == ';')
01861             cerr << "select_object: reached end of string\n";
01862         if (++i >= BUFF_LENGTH)
01863             cerr << "select_object: at position #1: buffer length exceeded\n";
01864         }
01865 
01866     if (level != 0)
01867         {
01868         cerr << "select_object: number of (,[,{,< and of ),],},> unbalanced\n";
01869         exit(1);
01870         }
01871     j = 0;
01872     while (c = hier_string[i])
01873         {
01874         if (c == '(' || c == '[' || c == '{' || c == '<')
01875             level++;
01876         else if (c == ')' || c == ']' || c == '}' || c == '>')
01877             level--;
01878         else if ((c == ',' && level == 0) || c == ';')
01879             {
01880             sprintf(tmp_string + j, ";");      /* end with ';' and  NULL  !! */
01881             if (level != 0)
01882                 {
01883                 cerr << "select_object: number of (,[, etc. unbalanced\n";
01884                 exit(1);
01885                 }
01886             break;
01887             }
01888         tmp_string[j++] = c;
01889 
01890         if (++i >= BUFF_LENGTH)
01891             cerr << "select_object: at position #2: buffer length exceeded\n";
01892         }
01893     }
01894 
01895 /*-----------------------------------------------------------------------------
01896  *  append_object  --  
01897  *               NOTE: only one-digit star numbering implemented, i.e.
01898  *                     only valid for an N-body sub-system with N < 11
01899  *-----------------------------------------------------------------------------
01900  */
01901 local void  append_object(char hier_string[BUFF_LENGTH],
01902                           char tmp_string[BUFF_LENGTH])
01903     {
01904     char *last_endptr;
01905     char *new_startingptr;
01906     char *end_of_string(char *);
01907 
01908     if (hier_string[0] == '\0')
01909         new_startingptr = hier_string;
01910     else
01911         {
01912         last_endptr = end_of_string(hier_string) - 1;
01913         if (*last_endptr == ';')
01914             *last_endptr = ',';                      /* insert new delimiter */
01915         else        
01916             cerr << "append_object: hier_string not properly terminated\n";
01917         new_startingptr = last_endptr + 1;
01918         }
01919 
01920     sprintf(new_startingptr, tmp_string);
01921     }
01922 
01923 /*-----------------------------------------------------------------------------
01924  *  unwrap_string  --  takes off the leading and trailing bracket,
01925  *                     if present, and only if the string 'substring' contains
01926  *                     only one object (the first two 'return;' statements
01927  *                     take care of these two exceptions).
01928  *                     If no enclosing brackets are present, the head and tail
01929  *                     characters '*head_ptr' and '*tail_ptr' are assigned the
01930  *                     value '*' to indicate the absence of brackets.
01931  *               NOTE: only one-digit star numbering implemented, i.e.
01932  *                     only valid for an N-body sub-system with N < 11
01933  *-----------------------------------------------------------------------------
01934  */
01935 local void  unwrap_string(char substring[BUFF_LENGTH],
01936                           char * head_ptr, char * tail_ptr)
01937     {
01938     int  i;
01939     int  level;
01940     char  c;
01941 
01942     c = substring[0];
01943 
01944     if (c == '(' || c == '[' || c == '{' || c == '<')
01945         *head_ptr = c;
01946     else
01947         {
01948         *head_ptr = *tail_ptr = '*';
01949         return;                         /* no enclosing brackets of any kind */
01950         }
01951 
01952     level = 1;
01953     i = 1;
01954     for (;;)
01955         {
01956         c = substring[i];
01957         if (c == '(' || c == '[' || c == '{' || c == '<')
01958             level++;
01959         else if (c == ')' || c == ']' || c == '}' || c == '>')
01960             level--;
01961 
01962         if (level == 0)
01963             {
01964             if (substring[i+1] != ';')    /* more than one object,           */
01965                 {                         /* therefore no enclosing brackets */
01966                 *head_ptr = *tail_ptr = '*';
01967                 return;
01968                 }
01969             else
01970                 break;
01971             }
01972 
01973         if (++i >= BUFF_LENGTH)
01974             cerr << "unwrap_string: buffer length exceeded\n";
01975         }
01976 
01977     i = 1;
01978     while ((c = substring[i]) != ';')
01979         {
01980         substring[i-1] = c;
01981         i++;
01982         }
01983     *tail_ptr = substring[i-2];
01984     substring[i-2] = ';';
01985     substring[i-1] = '\0';                    /* proper string ending */
01986     }
01987 
01988 /*-----------------------------------------------------------------------------
01989  *  wrap_string  --  puts the leading and trailing bracket in place again,
01990  *                   that is, if there are any brackets.
01991  *       convention: the character '*' for the value of 'head_char' and 
01992  *                   'tail_char' indicates the absence of brackets.
01993  *             NOTE: only one-digit star numbering implemented, i.e.
01994  *                   only valid for an N-body sub-system with N < 11
01995  *-----------------------------------------------------------------------------
01996  */
01997 local void  wrap_string(char substring[BUFF_LENGTH],
01998                         char head_char, char tail_char)
01999     {
02000     int  i;
02001     char  last_c;
02002     char  c;
02003 
02004     if (head_char == '*')
02005         return;
02006 
02007     last_c = head_char;
02008     i = 0;
02009     while((c = substring[i]) != ';')
02010         {
02011         substring[i] = last_c;
02012         last_c = c;
02013         i++;
02014         }
02015     substring[i] = last_c;
02016     substring[++i] = tail_char;
02017     substring[++i] = ';';
02018     substring[++i] = '\0';
02019     }
02020 
02021 main(int argc, char ** argv)
02022 {
02023     // silly trick to keep the compiler happy:
02024     char dummy = **argv;
02025 
02026     check_help();
02027 
02028     if (argc != 1) {
02029         cerr << "molecules: no arguments allowed\n";
02030         exit(1);
02031     } 
02032 
02033     dyn *b =  get_dyn(cin);
02034 
02035     printf("\n");
02036 
02037     while (b) {                 // return quietly if no new body is found
02038 
02039         dyn *bi;
02040 
02041         for (bi = b->get_oldest_daughter(); bi != NULL;
02042              bi = bi->get_younger_sister())
02043             if (bi->get_oldest_daughter() != NULL)
02044                 {
02045                 cerr << "molecules: flat tree required, sorry!\n";
02046                 exit(1);
02047                 }
02048 
02049         if (find_qmatch(b->get_dyn_story(), "t"))
02050             printf("Time = %15g   \n", getrq(b->get_dyn_story(), "t"));
02051 
02052         find_group_internal_structure(b);
02053 
02054         bi = b->get_oldest_daughter();
02055         while (bi) {
02056             dyn * tmp = bi->get_younger_sister();
02057             delete bi;
02058             bi = tmp;
02059         }
02060         delete b;
02061 
02062         b = get_dyn(cin);
02063     }
02064 }
02065 
02066 #endif
02067 
02068 // endof: molecules.C
Generated at Sun Feb 24 09:57:10 2002 for STARLAB by doxygen1.2.6 written by Dimitri van Heesch, © 1997-2001