---

encountervelocity.C

Go to the documentation of this file.

//  sys_stats.C: Print out various diagnostic statistics on a system.

#include "dyn.h"
#include "star/stdfunc.h"

// cerr << "    [R]: " << Rsun_pc / r_conv_star_to_dyn << " pc\n";
#define Rsun_pc 2.255e-8                 // R_sun/parsec = 6.960e+10/3.086e+18

#define ALL_i ni = n->get_oldest_daughter(); ni != NULL; \
                                             ni = ni->get_younger_sister()
#define j_ABOVE_i nj = ni->get_younger_sister(); nj != NULL; \
                                             nj = nj->get_younger_sister()

#ifdef TOOLBOX

local real conv_v_dyn_to_star(real v, real rf, real tf) {

//              Internal velocity is km/s
//              Internal size is solar raii
//              Internal time is Myr.
//              km/s to Rsun/Myr
//      real to_Rsun_Myr = cnsts.physics(km_per_s) * cnsts.physics(Myear)
//                     / cnsts.parameters(solar_radius);

      real myear = 3.15e+13;
      real km_s = 1.0e+5;
      real R_sun = 6.960e+10;
      real to_Rsun_Myr = km_s * myear / R_sun;
      real to_dyn      = rf/tf;
      
      return v/(to_Rsun_Myr * to_dyn);
   }


typedef  struct
{
    real  radius;
    real  mass;
} rm_pair, *rm_pair_ptr;

//-----------------------------------------------------------------------------
//  compare_radii  --  compare the radii of two particles
//-----------------------------------------------------------------------------

local int compare_radii(const void * pi, const void * pj)
{
    if (((rm_pair_ptr) pi)->radius > ((rm_pair_ptr) pj)->radius)
        return(1);
    else if (((rm_pair_ptr)pi)->radius < ((rm_pair_ptr)pj)->radius)
        return(-1);
    else
        return(0);
}

//-----------------------------------------------------------------------------
//  compute_general_mass_radii  --  Get the massradii for all particles.
//-----------------------------------------------------------------------------

static real nonlin_masses[9] = {0.005, 0.01, 0.02, 0.05, 0.1,
                                0.25, 0.5, 0.75, 0.9};

local void compute_general_mass_radii(dyn * b, vector dc_pos, int nzones,
                                      real r_lagr[],
                                      real m_cutoff, real M_cutoff,
                                      bool nonlin = false, 
                                      boolfn bf=NULL)
{
    if (nzones < 2) return;
    if (nonlin && nzones != 10) return;         // Special case

    // Note: nzones specifies the number of radial zones to consider.
    //       However, we only store radii corresponding to the nzones-1 
    //       interior zone boundaries.

    real* mass_percent = new real[nzones-1];
    if (mass_percent == NULL) {
        cerr << "compute_general_mass_radii: "
             << "not enough memory left for mass_percent\n";
        return;
    }

    int n = 0;
    if (bf == NULL) {
      for_all_daughters(dyn, b, bb)
        if (bb->get_mass()>=m_cutoff && bb->get_mass()<=M_cutoff)
          n++;
    }
    else {
        for_all_daughters(dyn, b, bb)
          if (bb->get_mass()>=m_cutoff && bb->get_mass()<=M_cutoff)
            if ((*bf)(bb))
              n++;
    }

    rm_pair_ptr rm_table = new rm_pair[n];

    if (rm_table == NULL) {
        cerr << "compute_general_mass_radii: "
             << "not enough memory left for rm_table\n" << flush;
        return;
    }

    // Use the geometric center if the density center is unknown.

    //vector dc_pos = 0;
    //if (find_qmatch(b->get_dyn_story(), "density_center_pos")) 
    //dc_pos = getvq(b->get_dyn_story(), "density_center_pos");

    // Set up an array of (radius, mass) pairs.  Also find the total
    // mass of all nodes under consideration.

    real total_mass = 0;

    int i = 0;
    for_all_daughters(dyn, b, bi) {
      if (bf == NULL || (*bf)(bi)) {
        if (bi->get_mass()>=m_cutoff && bi->get_mass()<=M_cutoff) {
          total_mass += bi->get_mass();
          rm_table[i].radius = abs(bi->get_pos() - dc_pos);
          rm_table[i].mass = bi->get_mass();
          i++;
//          PRL(i);
        }
      }
    }   

    // Sort the array by radius.

    qsort((void *)rm_table, (size_t)n, sizeof(rm_pair), compare_radii);

    // Determine Lagrangian radii.

    //    cerr << "Determine Lagrangian radii" << endl;
    int k;
    for (k = 0; k < nzones-1; k++) {
        if (!nonlin) 
            mass_percent[k] = ((k + 1) / (real)nzones) * total_mass;
        else
            mass_percent[k] = nonlin_masses[k] * total_mass;
    }

    //    real *rlagr = new real[nzones-1];
    real cumulative_mass = 0.0;
    i = 0;


    for (k = 0; k < nzones-1; k++) {

      while (cumulative_mass < mass_percent[k]) {
        cumulative_mass += rm_table[i++].mass;
      }

      r_lagr[k] = rm_table[i-1].radius;
    }


    // Place the data in the root dyn story.

#if 0
    if (bf == NULL)
        putiq(b->get_dyn_story(), "boolfn", 0);
    else
        putiq(b->get_dyn_story(), "boolfn", 1);

    putiq(b->get_dyn_story(), "n_nodes", n);

    putiq(b->get_dyn_story(), "n_lagr", nzones-1);
    putra(b->get_dyn_story(), "r_lagr", rlagr, nzones-1);
#endif

    delete mass_percent;
    delete rm_table;
//    delete rlagr;
    
    for (k = 0; k < nzones-1; k ++) 
      cerr << " " << b->get_starbase()->conv_r_dyn_to_star(r_lagr[k])
                     * Rsun_pc;
    cerr << " [pc]" << endl;
    
}

// Convenient synonyms:

local void  compute_mass_radii_quartiles(dyn * b, vector dc_pos,
                                         real m_cutoff, real M_cutoff)
{
  int nl = 4;
  real *r_lagr = new real[nl-1];
  for (int i=0; i<nl-1; i++)
    r_lagr[i] = 0;

    compute_general_mass_radii(b, dc_pos, nl, r_lagr, m_cutoff, M_cutoff);
}

local void  compute_mass_radii_percentiles(dyn * b, vector dc_pos,
                                           real m_cutoff, real M_cutoff)
{
  int nl = 10;
  real *r_lagr = new real[nl-1];
  for (int i=0; i<nl-1; i++)
    r_lagr[i] = 0;
  
    compute_general_mass_radii(b, dc_pos, nl, r_lagr, m_cutoff, M_cutoff);
}

local bool single_fn(dyn * b)
{

  //  int test = getiq(b->get_log_story(), "mass_doubled");
  //  PRL(test);
  if (getiq(b->get_log_story(), "mass_doubled") == 0 )
    return true;
  else
    return false;
}

local bool double_fn(dyn * b)
{
  if (getiq(b->get_log_story(), "mass_doubled") == 1 )
    return true;
  else
    return false;
}  

local void print_lagrangian_radii(dyn* b, vector dc_pos, int which_lagr,
                                  real m_cutoff, real M_cutoff,
                                  int which = 0)
{
    bool nonlin = false;

    int nl=0, indent=0;
    if (which_lagr == 0) {
        nl = 4;
        indent = 15;
        PRI(4); cerr << "quartiles: ";
    } else if (which_lagr == 1) {
        nl = 10;
        indent = 20;
        PRI(4); cerr << "10-percentiles: ";
    } else {
        nl = 10;
        indent = 26;
        PRI(4); cerr << "selected percentiles: ";
        nonlin = true;
    }
    
    //cerr << "before testnode call" << endl;
    //    if (which != 2 || 3) compute_general_mass_radii(b, nl, nonlin, testnode);

    real* r_lagr = new real[nl-1];
    if (r_lagr==NULL) {
      cerr << "print_lagrangian_radii not enough memory to allocate "
           << "r_lagr[" << nl-1<< "]" << endl;
      return;
    }
    for (int i=0; i<nl-1; i++)
      r_lagr[i] = 0;
    
    if (which == 0) compute_general_mass_radii(b, dc_pos, nl, r_lagr,
                                               m_cutoff, M_cutoff, nonlin); 

    //    cerr << "before single_fn call" << endl;
    //    PRL(which);
    if (which == 2) compute_general_mass_radii(b, dc_pos, nl, r_lagr,
                                               m_cutoff, M_cutoff, nonlin, single_fn);

    //     cerr << "before double_fn call" << endl;
    //    PRL(which);
    if (which == 3) compute_general_mass_radii(b, dc_pos, nl, r_lagr,
                                               m_cutoff, M_cutoff, nonlin, double_fn);

    //    PRC(dc_pos);PRC(m_cutoff);PRC(M_cutoff);PRL(which);

    //for (int k = 0; k < nl; k ++) 
    // cerr << " " << r_lagr[k];
    //cerr <<endl;

#if 0
    if (find_qmatch(b->get_dyn_story(), "n_lagr")) {
        int n_lagr = getiq(b->get_dyn_story(), "n_lagr");
        getra(b->get_dyn_story(), "r_lagr", r_lagr, n_lagr);

        for (int k = 0; k < n_lagr; k += 5) {
          if (k > 0) {
            cerr << endl;
            for (int kk = 0; kk < indent; kk++) cerr << " ";
          }
          for (int i = k; i < k+5 && i < n_lagr; i++)
            cerr << " " << r_lagr[i];
        }
        cerr << endl;
    }
#endif
    
    delete r_lagr;
}

/*-----------------------------------------------------------------------------
 *  compute_density  --  Get the density for all particles.
 *               note:
 *                    the neighbor_dist_sq[q] array will contain the distances
 *                    to the q-th nearest neighbor (in the form of the squares 
 *                    thereof); therefore  neighbor_dist_sq[0] = 0.0  since
 *                    this indicates the distance of a particle to itself.
 *                    Including this trivial case simplifies the algorithm
 *                    below by making it more homogeneous.
 *-----------------------------------------------------------------------------
 *              Litt.:
 *                    Stefano Casertano and Piet Hut: Astroph.J. 298,80 (1985).
 *                    To use their recipe, k >= 2 is required.
 *-----------------------------------------------------------------------------
 */
local void  compute_density(dyn * b,    /* pointer to an Nbody system */
                            int k,      /* use kth nearest neighbor */
                            int n,
                            real local_density[],
                            real m_cutoff,
                            real M_cutoff)
{
    int  q, r, is = 0;
    real *neighbor_dist_sq;
    real  volume;
    real  delr_sq;
    
    if (k >= n)                           /* no k-th nearest neighbor exists */
        cerr << "compute_density: k = " << k << " >= n = " << n << endl;
    if (k <= 1)
        cerr << "compute_density: k = " << k << " <= 1" << endl;

    neighbor_dist_sq = new real[k+1];

    //  for_all_leaves(dyn, b, d)               // The density is now defined
    for_all_daughters(dyn, b, d)        // ONLY for top-level nodes.
      if (d->get_mass()>=m_cutoff && d->get_mass()<=M_cutoff)
        {
        neighbor_dist_sq[0] = 0.0;
        for (q = 1; q <= k; q++)
            neighbor_dist_sq[q] = VERY_LARGE_NUMBER;

        //      for_all_leaves(dyn, b, dd)
        for_all_daughters(dyn, b, dd)
          if (dd->get_mass()>=m_cutoff && dd->get_mass()<=M_cutoff)
            {
            if (d == dd)
                continue;

            delr_sq = square(something_relative_to_root(d, &dyn::get_pos)
                             - something_relative_to_root(dd, &dyn::get_pos));

            if (delr_sq < neighbor_dist_sq[k])
                for (q = k-1; q >= 0; q--)
                    {
                    if (delr_sq > neighbor_dist_sq[q])
                        {
                        for (r = k; r > q+1; r--)
                            neighbor_dist_sq[r] = neighbor_dist_sq[r-1];
                        neighbor_dist_sq[q+1] = delr_sq;
                        break;
                        }
                    }
            }
            
        volume = (4.0/3.0) * PI * pow(neighbor_dist_sq[k], 1.5);

        real density =  (k - 1) / volume;           /* Ap. J. 298, 80 (1985) */
        story * s = d->get_dyn_story();
        //        putiq(s, "density_k_level", k);
        //        putrq(s, "density", density);
        local_density[is] = density;
        is++;
        }

    delete neighbor_dist_sq;
}



/*-----------------------------------------------------------------------------
 *  compute_mean_cod -- Returns the position and velocity of the density
 *                      center of an N-body system.
 *-----------------------------------------------------------------------------
 */
local void compute_mean_cod(dyn *b,
                            real local_density[],
                            vector& pos, vector& vel, real m_cutoff,
                            real M_cutoff)
{
    real total_weight = 0;
    pos = 0;
    vel = 0;
    int is=0;
    
    for_all_daughters(dyn, b, d) 
      if (d->get_mass()>=m_cutoff && d->get_mass()<=M_cutoff) {
            
        //      real density = getrq(d->get_dyn_story(), "density");
        real density = local_density[is++];
        real dens2 = density * density;        // weight factor

        total_weight += dens2;
        pos += dens2 * d->get_pos();
        vel += dens2 * d->get_vel();
    }   

    pos /= total_weight;
    vel /= total_weight;
}

/*===========================================================================*/

local void compute_core_parameters(dyn* b, int k, vector& center,
                                   real& rcore, int& ncore, real& mcore,
                                   real & vcore, int n,
                                   real local_density[], 
                                   real m_cutoff, real M_cutoff)
{
    vector vel;
    int is=0;

    real rcore2 = 0;
    if (rcore<=0) {
      compute_density(b, k, n, local_density, m_cutoff, M_cutoff);
      compute_mean_cod(b, local_density, center, vel, m_cutoff, M_cutoff);

      real total_weight = 0;
      real sum = 0;
      for_all_daughters(dyn, b, bi) 
        if (bi->get_mass()>=m_cutoff && bi->get_mass()<=M_cutoff) {
          real density = local_density[is++];
          real dens2 = density * density;        // weight factor
          //real density = getrq(bi->get_dyn_story(), "density");
          //real dens2 = density * density;              // weight factor

          total_weight += dens2;
          sum += dens2 * square(bi->get_pos() - center);
        }

      rcore2 = sum/total_weight;
      rcore = sqrt(rcore2);
    }
    else
      rcore2 = square(rcore);

    ncore = 0;
    mcore = 0;
    vcore = 0;

    for_all_daughters(dyn, b, bj)
      if (bj->get_mass()>=m_cutoff && bj->get_mass()<=M_cutoff) 
        if (square(bj->get_pos() - center) <= rcore2) {
            ncore++;
            mcore += bj->get_mass();
            vcore += bj->get_vel()*bj->get_vel();
        }
    mcore = mcore/(1.0*ncore);
    vcore = sqrt(vcore)/(1.0*ncore);
}

local void print_core_parameters(dyn* b, vector& density_center, real& rcore,
                                 int & ncore, real & vcore,
                                 int n, real local_density[],
                                 real m_cutoff, real M_cutoff,
                                 bool verbose)
  {
    real mcore;

    compute_core_parameters(b, 6, density_center, rcore, ncore, mcore,
                            vcore, n, local_density,
                            m_cutoff, M_cutoff);

    real r_density_center = sqrt(density_center*density_center);

    vcore = conv_v_dyn_to_star(vcore,
                               b->get_starbase()->conv_r_star_to_dyn(1),
                               b->get_starbase()->conv_t_star_to_dyn(1));
    if (verbose) {

      real M_lower = b->get_starbase()->conv_m_dyn_to_star(m_cutoff);
      real M_upper = b->get_starbase()->conv_m_dyn_to_star(M_cutoff);

      cerr << "n = " << n << ", M (lower/upper) = ( "
                          << M_lower << ",  " << M_upper << " ) [Msun]\n";
      cerr << "     Core: (N, R, m, v) = ( " << ncore << ", "
           << Rsun_pc * b->get_starbase()->conv_r_dyn_to_star(rcore)
           << ", "  
           << b->get_starbase()->conv_m_dyn_to_star(mcore) <<", "
           << vcore
           << " ) [solar/km/s]\n";
      cerr << "     R_to_density_center = "
           << Rsun_pc *
              b->get_starbase()->conv_r_dyn_to_star(r_density_center)
           << " [pc]\n";
      

      /*
      //place the data in the root dyn-story.
      putiq(b->get_dyn_story(), "n_cutoff", n);
      putrq(b->get_dyn_story(), "m_cutoff", m_cutoff);
      putrq(b->get_dyn_story(), "M_cutoff", M_cutoff);
      putvq(b->get_dyn_story(), "density_center_pos", density_center);
      putrq(b->get_dyn_story(), "core_radius", rcore);
      putiq(b->get_dyn_story(), "n_core", ncore);
      putrq(b->get_dyn_story(), "m_core", mcore);
      */
    }
    else {
              PRC(n); PRC(m_cutoff); PRL(M_cutoff);
      PRI(4); PRC(rcore); PRC(ncore); PRL(mcore);
      PRI(4); PRL(r_density_center);
    }
        
}

local real system_energy(dyn* b, real m_cutoff, real M_cutoff)
{
    real kin = 0;
    for_all_leaves(dyn, b, bj)
      if (bj->get_mass()>=m_cutoff && bj->get_mass()<=M_cutoff)
        kin += bj->get_mass()
             * square(something_relative_to_root(bj, &dyn::get_vel));
    kin *= 0.5;

    real pot = 0.0;
    for_all_leaves(dyn, b, bi) 
      if (bi->get_mass()>=m_cutoff && bi->get_mass()<=M_cutoff) {
        real dpot = 0.0;
        for_all_leaves(dyn, b, bj) 
          if (bj->get_mass()>=m_cutoff && bj->get_mass()<=M_cutoff) {
            if (bj == bi) break;
            vector dx = something_relative_to_root(bi, &dyn::get_pos)
                          - something_relative_to_root(bj, &dyn::get_pos);
            dpot += bj->get_mass() / abs(dx);
        }
        pot -= bi->get_mass() * dpot;
    }

    return kin + pot;
}

local void get_energies(dyn * n, real eps2, 
                  real& kinetic_energy, real& potential_energy,
                  real m_cutoff, real M_cutoff)
{
    dyn  * ni, *nj;

    kinetic_energy = 0;
    for (ALL_i) 
      if (ni->get_mass()>=m_cutoff && ni->get_mass()<=M_cutoff) {
      
        vector v = ni->get_vel();
        kinetic_energy += ni->get_mass() * v * v;
    }
    kinetic_energy *= 0.5;

    potential_energy = 0;
    for (ALL_i) 
      if (ni->get_mass()>=m_cutoff && ni->get_mass()<=M_cutoff) {
        real dphi = 0;
        vector xi = ni->get_pos();
        for (j_ABOVE_i) 
          if (nj->get_mass()>=m_cutoff && nj->get_mass()<=M_cutoff) {
            vector xij = nj->get_pos() - xi;
            dphi += nj->get_mass()/sqrt(xij*xij + eps2);
        }
        potential_energy -= ni->get_mass() * dphi;
    }
}


local void print_energies(dyn* b, real& kT, real m_cutoff, real M_cutoff, bool verbose)
{
    // Energies (top-level nodes):

    real kinetic_energy, potential_energy;
    get_energies(b, 0.0, kinetic_energy, potential_energy, m_cutoff, M_cutoff);

    real total_energy = kinetic_energy + potential_energy;
    real virial_ratio = -kinetic_energy / potential_energy;
    kT = -total_energy / (1.5*b->n_daughters());

    if (verbose) {
      PRI(4); cerr << "top-level nodes:\n";
      PRI(8); PRC(kinetic_energy); PRL(potential_energy);
      PRI(8); PRC(total_energy); PRC(kT); PRL(virial_ratio);

      PRI(4); cerr << "total energy (full system) = "
                   << system_energy(b, m_cutoff, M_cutoff) << endl;
    }
    else
      cerr <<"\t"<< kinetic_energy <<" "<< potential_energy
           <<" "<< total_energy <<" "<< kT
           <<" "<< virial_ratio <<" "<< system_energy(b, m_cutoff, M_cutoff) << endl;
}

local void print_relaxation_time(dyn* b, real m_cutoff, real M_cutoff, bool verbose)
{
    // Print the relaxation time

    real potential_energy, kinetic_energy;
    real r_virial, t_relax;
    real total_mass = 0;

    int n=0;
    for_all_leaves(dyn, b, bj) 
      if (bj->get_mass()>=m_cutoff && bj->get_mass()<=M_cutoff) {
        n++;
        total_mass += bj->get_mass();
    }

    get_energies(b, 0.0, kinetic_energy, potential_energy, m_cutoff, M_cutoff);

    r_virial = -0.5 * total_mass * total_mass / potential_energy;

    t_relax = 9.62e-2 * sqrt(r_virial * r_virial * r_virial / total_mass)
              * n / log10(0.4 * n);

    if (verbose) {
      PRI(4); cerr << "r_virial = "
                   << Rsun_pc *
                      b->get_starbase()->conv_r_dyn_to_star(r_virial) << endl;
      PRI(4); cerr << "t_relax = "
                   <<  b->get_starbase()->conv_t_dyn_to_star(t_relax)  << endl;
    }
    else
      cerr <<"\t"<< r_virial <<" "<< t_relax << endl;
}

local void get_composition(dyn* b,
                           real m_cutoff, real M_cutoff,
                           int& n, real& mmean,
                           real& rmean,
                           real& vdisp) {
  
  n=0;
  real T_eff, L_eff, R_eff;
  real m_total=0, v2_disp=0, r_total=0;
  
  for_all_daughters(dyn, b, bi)
    if (bi->get_mass()>=m_cutoff && bi->get_mass()<=M_cutoff) {
      n++;
      m_total += bi->get_mass();
      v2_disp += bi->get_vel()*bi->get_vel();
      //R_eff = getrq(bi->get_star_story(), "R_eff");
      //if (R_eff<=0) {

      T_eff = getrq(bi->get_star_story(), "T_eff");
      L_eff = getrq(bi->get_star_story(), "L_eff");
      R_eff = sqrt(max(0.1, (1130.*L_eff)/pow(T_eff, 4)));
            
      r_total += R_eff;
    }

  if (n>0) {
    mmean = b->get_starbase()->conv_m_dyn_to_star(m_total)/(1.0*n);
    rmean = r_total/(1.0*n);

    vdisp = sqrt(v2_disp)/(1.0*n);
    vdisp = conv_v_dyn_to_star(vdisp,
                               b->get_starbase()->conv_r_star_to_dyn(1),
                               b->get_starbase()->conv_t_star_to_dyn(1));
  }
  
}
//#else

main(int argc, char **argv)
{
    bool binaries = true, verbose = true, out = false,
         calc_e = true, Mcut=true, mcut=true;
    int nzones = 10;
    int nvzones = 1000;
    real magnify_r = 1;
    bool auto_detect = false;
    real m_min = 0.1;
    real m_max = 100;

    extern char *poptarg;
    int c;
    char* param_string = "n:Mmeor:v";   // Note: "v" removed because only the
                                        // "verbose" option currently works.

    while ((c = pgetopt(argc, argv, param_string)) != -1)
        switch(c) {

            case 'e': calc_e = !calc_e;
                      break;
            case 'n': nzones = atoi(poptarg);
                      break;
            case 'M': Mcut = !Mcut;
                      break;
            case 'm': mcut = !mcut;
                      break;
            case 'a': auto_detect = !auto_detect;
                      break;
            case 'o': out = !out;
                      break;
            case 'r': magnify_r = atof(poptarg);
                      break;
            case 'v': verbose = !verbose;
                      break;
            case '?': params_to_usage(cerr, argv[0], param_string);
                      exit(1);
        }

    // Loop over input until no more data remain.

    dyn *b;

    real to_volume = 4*PI/3.;
    real n_rho1, n_rhoi, n_rhoj;
    real r_tot, m_tot, v_rel;      
    int n, ncore;
    real mass_int;
    vector density_center;
    real rcore=0, vcore=0;
    real kT;
    real M_cutoff, m_cutoff;
    real m_total;

    real mmean, rmean, vdisp;
    //     real *core_radius = (real *) calloc(nzones+1, 8);
    //     int  *nstars = (int *) calloc(nzones+1, 4);
    //     real *m_mean = (real *) calloc(nzones+1, 8);
    //     real *v_disp = (real *) calloc(nzones+1, 8);
    //     real *r_mean = (real *) calloc(nzones+1, 8);
    
    real *core_radius = new real[nzones+1];
    int  *ncstars = new int[nzones+1];
    int  *nstars = new int[nzones+1];
    real *m_mean = new real[nzones+1];
    real *v_disp = new real[nzones+1];
    real *vc_disp = new real[nzones+1];
    real *r_mean = new real[nzones+1];


    real *pvenc = new real[nvzones];
    for (int i=0; i<nvzones; i++)
      pvenc[i] = 0;
    
    real geometric, grav_foc, encounter_rate, total_rate;
    
    while (b = get_dyn(cin)) {

      total_rate = rcore = vcore=0;
      
      for (int i=0; i<=nzones; i++) {
        core_radius[i] = m_mean[i] = v_disp[i] = vc_disp[i] = r_mean[i] = 0;
        ncstars[i] = nstars[i] = 0;
      }

      b->get_starbase()->print_stellar_evolution_scaling(cerr);
      cerr << "Time = " << b->get_system_time() << endl;
      
      // Count stars and find mass extremes.
      n=0;
      if (auto_detect) {
        m_min = VERY_LARGE_NUMBER;
        m_max = -m_min;
        for_all_daughters(dyn, b, bi) {
          n++;
          m_min = min(m_min, bi->get_mass());
          m_max = max(m_max, bi->get_mass());
        }
      }
      else {
        for_all_daughters(dyn, b, bi) 
          n++;
        m_min = b->get_starbase()->conv_m_star_to_dyn(m_min);
        m_max = b->get_starbase()->conv_m_star_to_dyn(m_max);
      }
      mass_int = (log10(m_max)-log10(m_min))/(1.0*nzones-1);

      //rcore = getrq(b->get_dyn_story(), "core_radius");
      //density_center = getrq(b->get_dyn_story(), "density_center_pos");
      
      real* local_density = new real[n+1];
      if(local_density==NULL) {
        cerr << "Not enough memory to allocate local_density[" << n << "].\n";
        exit(1);
      }
      for(int i=0; i<n; i++)
        local_density[i] = 0;

      print_core_parameters(b, density_center, rcore, ncore, vcore,
                            n, local_density, m_min, m_max, verbose);
      print_energies(b, kT, m_min, m_max, verbose);
      print_relaxation_time(b, m_min, m_max, verbose);
      print_lagrangian_radii(b, density_center, 0, m_min, m_max, 0);

      core_radius[0] = Rsun_pc * b->get_starbase()->conv_r_dyn_to_star(rcore);
      nstars[0] = n;
      ncstars[0] = ncore;
      vc_disp[0] = vcore;

      get_composition(b, m_min, m_max, n, mmean, rmean, vdisp);
      m_mean[0] = mmean;
      v_disp[0] = vdisp;
      r_mean[0] = rmean;
      
      //PRC(nstars[0]);PRC(m_min);PRC(m_max);PRL(mass_int);
      //PRC(rcore);PRC(m_mean[0]);PRC(r_mean[0]);PRC(v_disp[0]);

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

        m_cutoff = pow(10., log10(m_min)+i*mass_int);
        M_cutoff = pow(10., log10(m_min)+(i+1)*mass_int);

        if (!Mcut) M_cutoff = m_max;
        if (!mcut) m_cutoff = m_min;

        for(int is=0; is<nstars[0]; is++)
          local_density[is] = 0;

        mmean = vdisp = rmean = 0;
        get_composition(b, m_cutoff, M_cutoff, n, mmean, rmean, vdisp); 
        nstars[i+1] = n;
        m_mean[i+1] = mmean;
        v_disp[i+1] = vdisp;
        r_mean[i+1] = rmean;
        
        if (n > 6 &&
            !(n==nstars[i] && (!Mcut || !mcut))) {

          rcore = 0;
          print_core_parameters(b, density_center, rcore, ncore, vcore,
                                n, local_density, m_cutoff, M_cutoff, verbose);
          core_radius[i+1] = Rsun_pc * b->get_starbase()
                                        ->conv_r_dyn_to_star(rcore);
          ncstars[i+1] = ncore;
          vc_disp[i+1] = vcore;

          print_energies(b, kT, m_cutoff, M_cutoff, verbose);
          print_relaxation_time(b, m_cutoff, M_cutoff, verbose);
          print_lagrangian_radii(b, density_center, 0, m_cutoff, M_cutoff, 0);

          //PRC(n);PRC(core_radius[0]);PRL(core_radius[i+1]);

          n_rho1 = 1./(to_volume * pow(core_radius[0], 3));
          n_rhoi = ncstars[i+1]/(to_volume * pow(core_radius[i+1], 3));
          r_tot  = r_mean[0] + r_mean[i+1];
          m_tot  = m_mean[0] + m_mean[i+1];
          v_rel  = sqrt(pow(vc_disp[0], 2) + pow(vc_disp[i+1], 2));

          //PRC(vc_disp[0]);PRC(vc_disp[i+1]);PRL(v_rel);
          
          geometric = 6.31e-9 * v_rel * pow(r_tot, 2);
          grav_foc = 3.61e-3 * m_tot * r_tot / v_rel;
          encounter_rate = (n_rho1/1000.) * (n_rhoi/1000)
                         * pow(core_radius[0], 3)
                         * (grav_foc + geometric);

          if (verbose) {
            cerr << "     Stars: (N, r, m, v) = ( "
                 << nstars[i+1] <<", "
                 << r_mean[i+1] <<", "
                 << m_mean[i+1] <<", "
                 << v_disp[i+1] <<" ) [solar/km/s]\n";
            cerr << "     encounter_rate = "
                 << encounter_rate << " [per Myr]\n"
                 << "          (geo, gf) = ( "
                 << 100 * geometric/(geometric+grav_foc)
                 << ", " << 100 * grav_foc/(geometric+grav_foc)
                 << " ) [\%]\n\n";
          }
          else {
            PRI(4); PRC(m_mean[i+1]);PRC(v_disp[i+1]);PRL(r_mean[i+1]);
            PRI(4);PRC(geometric); PRL(grav_foc);
            PRI(4);PRL(encounter_rate);
            cerr << endl;

          }

        }
      }
      
      // Now print the encounter probability in a mass-mass plane.
      cerr << log10(b->get_starbase()
                     ->conv_m_dyn_to_star(m_min)) <<" "
           << log10(b->get_starbase()
                     ->conv_m_dyn_to_star(m_max)) << " "
           << nzones << endl;
      cerr << log10(b->get_starbase()
                     ->conv_m_dyn_to_star(m_min)) <<" "
           << log10(b->get_starbase()
                     ->conv_m_dyn_to_star(m_max)) << " "
           << nzones << endl;
      
      for (int im=1; im<nzones+1; im++) {
        for (int jm=im; jm<nzones+1; jm++) {

          //PRC(im);PRC(jm);PRC(nstars[im]);PRL(nstars[jm]);
          //PRC(core_radius[0]);PRC(core_radius[im]);PRL(core_radius[jm]);
            
          if (core_radius[im]>0 && core_radius[jm]>0) {


            n_rhoi = 1./(to_volume * pow(core_radius[im], 3));
            n_rhoj = 1./(to_volume * pow(core_radius[jm], 3));
            r_tot  = r_mean[im] + r_mean[jm];
            m_tot  = m_mean[im] + m_mean[jm];
            v_rel  = sqrt(pow(vc_disp[im], 2) + pow(vc_disp[jm], 2));

            rcore = min(core_radius[im], core_radius[jm]);

            //PRC(vc_disp[im]);PRC(vc_disp[jm]);PRL(v_rel);

            geometric = 6.31e-15 * v_rel * pow(r_tot, 2);
            grav_foc = 3.61e-9 * m_tot * r_tot / v_rel;

            if (im==jm) 
              encounter_rate = 0.5 * ((ncstars[im]-1) * n_rhoi)
                                   *  (ncstars[jm] * n_rhoj)
                             * pow(rcore, 3)
                             * (grav_foc + geometric);
            else
              encounter_rate = (ncstars[im] * n_rhoi)
                             * (ncstars[jm] * n_rhoj)
                             * pow(rcore, 3)
                             * (grav_foc + geometric);

            cerr << jm <<" "<< im <<" "
                 << m_mean[jm] <<" "<< m_mean[im] <<" "
                 << ncstars[jm] <<" "<< ncstars[im] <<" "
                 << encounter_rate << endl;

          // Compute the probability distribution for the encounter velocity
            real dv = 40/dynamic_cast(real, nvzones);
            for(int iv=0; iv<nvzones; iv++) // weighted on the encounter rate
              pvenc[iv] += encounter_rate   
                         * gravitational_focussed_velocity(m_mean[jm],
                                                           m_mean[im], 
                                                           magnify_r * r_tot,
                                                           iv*dv, v_rel);
          }
          else if (core_radius[im]>0 || core_radius[jm]>0) {

            real ncore_i = ncstars[im];
            real ncore_j = ncstars[jm];
            
            if (ncstars[im]==0) ncore_i = core_radius[0] * nstars[im];
            if (ncstars[jm]==0) ncore_j = core_radius[0] * nstars[jm];
            
            //n_rhoi = 1./(to_volume * pow(core_radius[0], 3));
            //n_rhoj = 1./(to_volume * pow(core_radius[0], 3));
            r_tot  = r_mean[im] + r_mean[jm];
            m_tot  = m_mean[im] + m_mean[jm];
            //v_rel  = sqrt(pow(vc_disp[0], 2) + pow(vc_disp[0], 2));

            if (core_radius[im]>0) {
              n_rhoi = 1./(to_volume * pow(core_radius[im], 3));
              
              if (core_radius[jm]>0) {
                n_rhoj = 1./(to_volume * pow(core_radius[jm], 3));
                rcore = min(core_radius[im], core_radius[jm]);
                v_rel  = sqrt(pow(vc_disp[im], 2) + pow(vc_disp[jm], 2));
              }
              else {
                n_rhoj = 1./(to_volume * pow(core_radius[0], 3));
                rcore = min(core_radius[im], core_radius[0]);
                v_rel  = sqrt(pow(vc_disp[im], 2) + pow(v_disp[jm], 2));
              }
            }
            else {
              n_rhoi = 1./(to_volume * pow(core_radius[0], 3));
              
              if (core_radius[jm]>0) {
                n_rhoj = 1./(to_volume * pow(core_radius[jm], 3));
                rcore = min(core_radius[0], core_radius[jm]);
                v_rel  = sqrt(pow(v_disp[jm], 2) + pow(vc_disp[jm], 2));
              }
              else {
                rcore = core_radius[0];
                n_rhoj = 1./(to_volume * pow(rcore, 3));
                v_rel  = sqrt(pow(v_disp[im], 2) + pow(v_disp[jm], 2));
              }
            }

            //PRC(v_rel);PRC(n_rhoi);PRC(n_rhoj);PRL(rcore);
            //PRC(ncore_i);PRL(ncore_j);
            
            geometric = 6.31e-15 * v_rel * pow(r_tot, 2);
            grav_foc = 3.61e-9 * m_tot * r_tot / v_rel;
            if (im==jm) {
              if (ncore_i>1) 
                encounter_rate = 0.5 * ((ncore_i-1) * n_rhoi)
                                     *  (ncore_j * n_rhoj)
                               * pow(rcore, 3)
                               * (grav_foc + geometric);
              else
                encounter_rate = 0;
            }
            else
              encounter_rate = (ncore_i * n_rhoi)
                             * (ncore_j * n_rhoj)
                             * pow(rcore, 3)
                             * (grav_foc + geometric);

            cerr << jm <<" "<< im <<" "
                 << m_mean[jm] <<" "<< m_mean[im] <<" "
                 << ncore_j <<" "<< ncore_i <<" "
                 << encounter_rate << endl;

          // Compute the probability distribution for the encounter velocity
            real dv = 40/dynamic_cast(real, nvzones);
            for(int iv=0; iv<nvzones; iv++) // weighted on the encounter rate
              pvenc[iv] += encounter_rate   
                         * gravitational_focussed_velocity(m_mean[jm],
                                                           m_mean[im], 
                                                           magnify_r * r_tot,
                                                           iv*dv, v_rel);
          }
          
          total_rate += encounter_rate;
          cerr << flush;
          }
        }
      cerr << "     Total encounter rate = "
           << total_rate << " [per Myr]" << endl;
      
      cerr << 0 << " " << 0 << " "
           << 0 << " " << 0 << " "
           << 0 << " " << 0 << " "
           << 0 << endl;

      real max_pvenc = 0;
      for(int iv=0; iv<nvzones; iv++)
        max_pvenc = max(max_pvenc, pvenc[iv]);
      cerr << "Encounter velocity poribability distribution"<<endl;
      real dv = 40./dynamic_cast(real, nvzones);
      for(int iv=0; iv<nvzones; iv++) {
        pvenc[iv] = pvenc[iv]/max_pvenc;
        if (pvenc[iv]<=1e-7)
          pvenc[iv] = 0;
        cerr << iv*dv << " " << pvenc[iv] << endl;
      }
      
      rmtree(b);
      }
}

#endif

---