dyn_util.C Source File

00001 
00002 #include "hdyn.h"
00003 #include "xstarplot.h"
00004 
00005 void convert_relative_to_absolute(hdyn* b)
00006 {
00007     if (b->get_parent()) b->inc_pos(b->get_parent()->get_pos());
00008     for_all_daughters(hdyn, b, bb) convert_relative_to_absolute(bb);
00009 }
00010 
00011 void accumulate_potential_energy(hdyn* bj, hdyn*bi,
00012                                  real& epot, real& rmin, hdynptr& bmin)
00013 
00014 // Determine the potential energy of bi with respect to bj, along
00015 // with bi's nearest neighbor and minimum neighbor distance.
00016 
00017 {
00018     if (bj->get_oldest_daughter() != (hdyn*)NULL)
00019         for (hdyn * bb = bj->get_oldest_daughter(); bb != (hdyn*)NULL;
00020             bb = bb->get_younger_sister()) {
00021             accumulate_potential_energy(bb, bi, epot, rmin, bmin);
00022         }
00023     else
00024         if (bi != bj) {
00025             vector d_pos = bi->get_pos() - bj->get_pos();
00026             real mi = bi->get_mass();
00027             real mj = bj->get_mass();
00028             real r = sqrt(d_pos * d_pos);
00029             epot += -mi*mj/r;
00030             if (r < rmin) {rmin = r; bmin = bj;}
00031         }
00032 }
00033 
00034 // compute_energies: calculate the appropriate color code for particle bi
00035 //                   relative to "particle" bj (usually the root node).
00036 
00037 void compute_energies(hdyn* bj, hdyn* bi, char& c)
00038 {
00039     c = default_color;
00040 
00041     real   epot = 0, rmin = VERY_LARGE_NUMBER;
00042     hdynptr bmin = bi;
00043 
00044     accumulate_potential_energy(bj, bi, epot, rmin, bmin);
00045 
00046     real   mi = bi->get_mass();
00047     real   mj = bmin->get_mass();
00048     vector d_vel = bi->get_vel() - bmin->get_vel();
00049     vector d_pos = bi->get_pos() - bmin->get_pos();
00050     real   r = sqrt(d_pos * d_pos);
00051     real   e0 = (0.5 * d_vel * d_vel - (mi + mj)/r);
00052 
00053     if (e0 < 0) {
00054         real   epot1 = 0, rmin1 = VERY_LARGE_NUMBER;
00055         hdynptr bmin1 = bi;
00056 
00057         accumulate_potential_energy(bj, bmin, epot1, rmin1, bmin1);
00058 
00059         if (bi == bmin1) {
00060             real e  = - mi*mj / r;
00061             vector R_vel = (mi*bi->get_vel()+mj*bmin->get_vel())/(mi+mj);
00062             real ekin = 0.5*(mi+mj)*R_vel*R_vel;
00063 
00064             if (epot + epot1 - 2*e + ekin < 0) c = bound_binary;
00065             else c = unbound_binary;
00066 
00067         } else c = bound_single;
00068 
00069     } else {
00070         vector vel = bi->get_vel();
00071         real ekin = 0.5*bi->get_mass()*vel*vel;
00072         
00073         if (ekin + epot > 0.0) c = unbound_single;
00074     }
00075 }
00076 
00077 int nearest_index(hdyn* b, float r, float s, int kx, int ky)
00078 {
00079     real r2_min = VERY_LARGE_NUMBER;
00080     int index = -1;
00081     for_all_leaves(hdyn, b, bb) {
00082         real dx = bb->get_pos()[kx] - r;
00083         real dy = bb->get_pos()[ky] - s;
00084         real r2 = dx*dx + dy*dy;
00085         if (r2 < r2_min) {
00086             r2_min = r2;
00087             index = bb->get_index();
00088         }
00089     }
00090     return index;
00091 }
00092 
00093 int count_nodes(hdyn* b)
00094 // Return the total number of nodes BELOW the node b.
00095 {
00096     if(b->get_oldest_daughter() == NULL)
00097         return 0;
00098     else {
00099         int  n = 0;
00100         for (hdyn* daughter = b->get_oldest_daughter();
00101              daughter != NULL; daughter = daughter->get_younger_sister())
00102             n += 1 + count_nodes(daughter);
00103         return n;
00104     }
00105 }