---

compute_max_cod.C

Go to the documentation of this file.

//-----------------------------------------------------------------------------
//   version 1:  May 1989   Piet Hut
//   version 2:  Nov 1994   Piet Hut
//   version 3:  Jul 1996   Steve McMillan & Jun Makino
//.............................................................................
//   non-local function: 
//      compute_max_cod
//.............................................................................
//   see also: density.C
//-----------------------------------------------------------------------------

#include "dyn.h"

#ifndef TOOLBOX

//-----------------------------------------------------------------------------
//  compute_max_cod -- Returns the position and velocity of the density
//                     center of an N-body system.
//-----------------------------------------------------------------------------

#define MAX_COUNT 5

void compute_max_cod(dyn *b, vector& pos, vector& vel)
{
    real max_density = 0;
    bool print_message = true;
    int count = 0;

    pos = 0;
    vel = 0;

//  for_all_leaves(dyn, b, d)
    for_all_daughters(dyn, b, d) {

        real dens_time = getrq(d->get_dyn_story(), "density_time");

        if (dens_time != b->get_system_time() && print_message) {
            warning("compute_max_cod: using out-of-date densities.");
            PRC(d->format_label()), PRL(dens_time);
            if (++count > MAX_COUNT) print_message = false;
        }

        real this_density = getrq(d->get_dyn_story(), "density");

        if (this_density > 0) {
            if (max_density < this_density) {
                max_density = this_density;
                pos = d->get_pos();
                vel = d->get_vel();
            }
        } else if (this_density <= -VERY_LARGE_NUMBER) {
            warning("compute_max_cod: density not set.");
            PRL(d->format_label());
            if (++count > MAX_COUNT) print_message = false;
        }
    }

    putsq(b->get_dyn_story(), "density_center_type", "max");
    putrq(b->get_dyn_story(), "density_center_time", b->get_system_time());
    putvq(b->get_dyn_story(), "density_center_pos", pos);
    putvq(b->get_dyn_story(), "density_center_vel", vel);
}

void compute_max_cod(dyn *b)
{
    vector pos, vel;
    compute_max_cod(b, pos, vel);
}

#else

/*-----------------------------------------------------------------------------
 *  main  --  driver to use compute_max_cod() as a tool
 *-----------------------------------------------------------------------------
 */
main(int argc, char ** argv)
{
    char  *comment;
    dyn * b;
    bool  c_flag = FALSE;       // if TRUE: a comment given on command line

    check_help();

    extern char *poptarg;
    int c;
    char* param_string = "c:";

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

            case 'c': c_flag = TRUE;
                      comment = poptarg;
                      break;
            case '?': params_to_usage(cerr, argv[0], param_string);
                      get_help();
                      exit(1);
        }            

    if ((b = get_dyn(cin)) == NULL)
       err_exit("compute_max_cod: No N-body system on standard input");

    while (b) {

        if (c_flag == TRUE)
            b->log_comment(comment);
        b->log_history(argc, argv);

        compute_max_cod(b);

        put_dyn(cout, *b);
        rmtree(b);
        b = get_dyn(cin);
    }
}

#endif

// endof: compute_max_cod.C

---