bound3.C Source File

00001 
00002        //=======================================================//    _\|/_
00003       //  __  _____           ___                    ___       //      /|\
00004      //  /      |      ^     |   \  |         ^     |   \     //          _\|/_
00005     //   \__    |     / \    |___/  |        / \    |___/    //            /|\
00006    //       \   |    /___\   |  \   |       /___\   |   \   // _\|/_
00007   //     ___/   |   /     \  |   \  |____  /     \  |___/  //   /|\     
00008  //                                                       //            _\|/_
00009 //=======================================================//              /|\
00010 
00053 
00054 // Input to the function is an initial state, output is intermediate
00055 // and final states.
00056 
00057 // Note the function bound3 is completely deterministic.
00058 // No randomization is performed at this level.
00059 
00060 // Starlab library function.
00061 
00062 #include "scatter3.h"
00063 
00064 #ifdef TOOLBOX
00065 
00066 // check_init:  Make sure an initial state is sensible.
00067 
00068 local void check_init(initial_state3 & init)
00069 {
00070     if (init.m2 > 1) err_exit("check_init: m2 > 1");
00071     if (init.m2 < 0) err_exit("check_init: m2 < 0");
00072     if (init.m3 > 1) err_exit("check_init: m3 > 0");
00073     if (init.m3 < 0) err_exit("check_init: m3 < 0");
00074     if (init.r1 < 0) err_exit("check_init: r1 < 0");
00075     if (init.r2 < 0) err_exit("check_init: r2 < 0");
00076     if (init.r3 < 0) err_exit("check_init: r3 < 0");
00077     if (init.eta <= 0) err_exit("check_init: eta <= 0");
00078 }
00079 
00080 // bound3:  Perform a three-body experiment, initializing from a specified
00081 //          state and returning intermediate- and final-state structures.
00082 
00083 void bound3(sdyn3 * b,
00084             initial_state3 & init,
00085             intermediate_state3 & inter,
00086             final_state3 & final,
00087             real cpu_time_check,
00088             real dt_out,         // diagnostic output interval
00089             real dt_snap,        // snapshot output interval
00090             real snap_cube_size) // limit output to particles within cube
00091 {
00092     inter.id = final.id = init.id;          // For bookkeeping...
00093 
00094     mksphere(b, 3);                         // Note that all systems are
00095     initialize_bodies(inter.system);        // initialized as "non-particles"
00096     initialize_bodies(final.system);        // until explicitly set.
00097 
00098     // Save some initial data:
00099 
00100     sdyn3_to_system(b, init.system);
00101 
00102 //    put_node(cerr, *b);
00103 //    print_bodies(cerr, init.system, STD_PRECISION);
00104 
00105     real e_init = energy(b);
00106     real cpu_init = cpu_time();
00107     real cpu = cpu_init;
00108     
00109     real kinetic, potential;
00110     get_top_level_energies(b, 0.0, potential, kinetic);
00111 
00112 //    PRC(e_init); PRC(potential); PRL(kinetic);
00113 
00114     // Evolve the system until the interaction is over or a collision occurs.
00115 
00116     inter.n_kepler = 0;
00117 
00118     do {
00119 
00120         int n_stars = 0;
00121         for_all_daughters(sdyn3, b, bb) n_stars++;
00122 
00123         real check_int = CHECK_INTERVAL;
00124 //      check_int = 0.5;
00125 
00126         tree3_evolve(b, check_int, dt_out, dt_snap, snap_cube_size,
00127                      init.eta, cpu_time_check); // , dt_print, p);
00128 
00129         // Check the CPU time.  Note that the printed CPU time is the
00130         // time since this routine was entered.
00131 
00132         if (cpu_time() - cpu > cpu_time_check) {
00133             if (cpu == cpu_init) cerr << endl; // (May be waiting in mid-line!)
00134 
00135             int p = cerr.precision(STD_PRECISION);
00136             cerr << "bound3:  CPU time = " << cpu_time() - cpu_init
00137                  << "  time = " << b->get_time()
00138                  << "  n_steps = " << b->get_n_steps()
00139                  << endl;
00140             cerr << "           n_osc = " << b->get_n_ssd_osc()
00141                  << "  n_kepler = " << inter.n_kepler
00142                  << endl << flush;
00143             cpu = cpu_time();
00144             cerr.precision(p);
00145         }
00146 
00147         inter.n_osc = b->get_n_ssd_osc();       // For convenience
00148 
00149         sdyn3_to_system(b, inter.system);
00150         merge_collisions(b);
00151 
00152         // See if the number of stars has decreased.
00153 
00154         for_all_daughters(sdyn3, b, bbb) n_stars--;
00155         if (n_stars > 0) {
00156             if (dt_snap < VERY_LARGE_NUMBER
00157                 && system_in_cube(b, snap_cube_size)) {
00158                 put_node(cout, *b);
00159                 cout << flush;
00160             }
00161         }
00162 
00163 //      cout << b->get_time() << " " << energy(b)-e_init << endl;
00164 
00165     } while (// b->get_n_steps() < 1000000 &&
00166              !extend_or_end_scatter(b, init, &inter, &final));
00167 
00168     // Set intermediate state:
00169 
00170     inter.r_min_min = VERY_LARGE_NUMBER;
00171     inter.n_stars = 0;
00172 
00173     for_all_daughters(sdyn3, b, bb) {
00174         inter.index[inter.n_stars] = bb->get_index();
00175         inter.r_min[inter.n_stars] = sqrt(bb->get_min_nn_dr2());
00176         inter.r_min_min = min(inter.r_min_min, inter.r_min[inter.n_stars]);
00177         inter.n_stars++;
00178     }
00179 
00180     if (b->get_n_ssd_osc() <= 1)
00181         inter.descriptor = non_resonance;
00182     else {
00183         int n_no_nn_change = 0;
00184         for_all_daughters(sdyn3, b, bbb)
00185             if (bbb->get_nn_change_flag() == 0)
00186                 n_no_nn_change++;
00187         if (n_no_nn_change >= 2)
00188             inter.descriptor = hierarchical_resonance;
00189         else
00190             inter.descriptor = democratic_resonance;
00191     }
00192 
00193     // Set final state:
00194 
00195     final.time = b->get_time();
00196     final.n_steps = (int) b->get_n_steps();
00197     final.error = energy(b) - e_init; // *Absolute* error, note!
00198 
00199     sdyn3_to_system(b, final.system);
00200 
00201     //  Check for integration errors (relax tolerance in the case of mergers).
00202 
00203     if (abs(final.error) > ENERGY_TOLERANCE)
00204         if ((   final.descriptor != merger_binary_1
00205              && final.descriptor != merger_binary_2
00206              && final.descriptor != merger_binary_3
00207              && final.descriptor != merger_escape_1
00208              && final.descriptor != merger_escape_2
00209              && final.descriptor != merger_escape_3
00210              && final.descriptor != triple_merger)
00211             || abs(final.error)
00212                  > MERGER_ENERGY_TOLERANCE * abs(potential_energy(b)))
00213             final.descriptor = error;
00214 
00215     // In all cases, the final system should be a properly configured
00216     // 1-, 2-, or 3-body system.  Print it for the last time, to reflect
00217     // final state, if the "-D" command-line option was specified.
00218 
00219     if (dt_snap < VERY_LARGE_NUMBER && system_in_cube(b, snap_cube_size)) {
00220         put_node(cout, *b);
00221         cout << flush;
00222     }
00223 }
00224 
00225 // Set up a template initial state.  We don't actually use all the
00226 // elements of the initial_state3 structure, but we will clean this
00227 // up later if necessary...
00228 
00229 void make_standard_bound_init(initial_state3 & init)
00230 {
00231     init.m2 = 1.0/3;            // mass of second particle
00232     init.m3 = 1.0/3;            // mass of third particle (m1 + m2 + m3 = 1)
00233     init.r1 = 0;                // radius of star #1
00234     init.r2 = 0;                // radius of star #2
00235     init.r3 = 0;                // radius of star #3
00236 
00237     init.r_stop
00238           = VERY_LARGE_NUMBER;  // final separation
00239     init.tidal_tol_factor = DEFAULT_TIDAL_TOL_FACTOR;
00240     init.eta
00241           = DEFAULT_ETA;        // accuracy parameter
00242 
00243     initialize_bodies(init.system);
00244     init.id = get_initial_seed() + get_n_rand();
00245 }
00246 
00247 main(int argc, char **argv)
00248 {
00249     initial_state3 init;
00250     make_standard_bound_init(init);
00251 
00252     int  seed       = 0;        // seed for random number generator
00253     int n_rand      = 0;        // number of times to invoke the generator
00254                                 // before starting for real
00255     int  n_experiments = 1;     // default: only one run
00256     real dt_out     =           // output time interval
00257           VERY_LARGE_NUMBER;
00258     real dt_snap    =           // output time interval
00259           VERY_LARGE_NUMBER;
00260 
00261     real cpu_time_check = 3600;
00262     real snap_cube_size = 10;
00263 
00264     int planar_flag = 0;
00265     real psi = 0;
00266 
00267     bool  b_flag = FALSE;
00268     bool  q_flag = FALSE;
00269     bool  Q_flag = FALSE;
00270     bool  i_flag = true;
00271     bool  s_flag = false;
00272 
00273     check_help();
00274 
00275     extern char *poptarg;
00276     int c;
00277     char* param_string = "A:bc:C:d:D:e:g:L:m:M:n:N:o:pPqQr:R:s:S:U:v:x:y:z:";
00278 
00279     while ((c = pgetopt(argc, argv, param_string)) != -1)
00280         switch(c) {
00281 
00282             case 'A': init.eta = atof(poptarg);
00283                       break;
00284             case 'b': b_flag = 1 - b_flag;
00285                       break;
00286             case 'c': cpu_time_check = 3600*atof(poptarg);// (Specify in hours)
00287                       break;
00288             case 'C': snap_cube_size = atof(poptarg);
00289                       break;
00290             case 'd': dt_out = atof(poptarg);
00291                       if (dt_out < 0) dt_out = pow(2.0, dt_out);
00292                       break;
00293             case 'D': dt_snap = atof(poptarg);
00294                       if (dt_snap < 0) dt_snap = pow(2.0, dt_snap);
00295                       break;
00296             case 'e': init.ecc = atof(poptarg);
00297                       break;
00298             case 'g': init.tidal_tol_factor = atof(poptarg);
00299                       break;
00300             case 'i': i_flag = !i_flag;
00301                       break;
00302             case 'L': init.r_init_min = atof(poptarg);
00303                       break;
00304             case 'm': init.m2 = atof(poptarg);
00305                       break;
00306             case 'M': init.m3 = atof(poptarg);
00307                       break;
00308             case 'n': n_experiments = atoi(poptarg);
00309                       break;
00310             case 'N': n_rand = atoi(poptarg);
00311                       break;
00312             case 'p': planar_flag = 1;
00313                       break;
00314             case 'P': planar_flag = -1;
00315                       break;
00316             case 'q': q_flag = 1 - q_flag;
00317                       break;
00318             case 'Q': Q_flag = 1 - Q_flag;
00319                       break;
00320             case 'r': init.rho = atof(poptarg);
00321                       break;
00322             case 'R': init.r_stop = atof(poptarg);
00323                       init.r_init_min = init.r_init_max = abs(init.r_stop);
00324                       break;
00325             case 's': seed = atoi(poptarg);
00326                       s_flag = true;
00327                       break;
00328             case 'S': init.r_stop = atof(poptarg);
00329                       break;
00330             case 'U': init.r_init_max = atof(poptarg);
00331                       break;
00332             case 'v': init.v_inf = atof(poptarg);
00333                       break;
00334             case 'x': init.r1 = atof(poptarg);
00335                       break;
00336             case 'y': init.r2 = atof(poptarg);
00337                       break;
00338             case 'z': init.r3 = atof(poptarg);
00339                       break;
00340             case '?': params_to_usage(cerr, argv[0], param_string);
00341                       get_help();
00342         }            
00343 
00344     if (Q_flag) q_flag = TRUE;
00345 
00346     if (init.m2 + init.m3 > 1) {
00347         cerr << "bound3: "; PRC(init.m2); PRL(init.m3);
00348         exit(1);
00349     }
00350 
00351     check_init(init);
00352 
00353     // Create the 3-body tree for all calculations.
00354 
00355     sdyn3 *b, *by, *bo;
00356     b = new sdyn3();
00357     bo = new sdyn3();
00358     if (i_flag) bo->set_label(1);
00359     b->set_oldest_daughter(bo);
00360     bo->set_parent(b);
00361 
00362     // Default from Simo:
00363 
00364     bo->set_mass(1);
00365     bo->set_pos(vector(-0.97000436, 0.24308753, 0));
00366     bo->set_vel(-0.5*vector(0.93240737, 0.86473146, 0));
00367 
00368     for (int i = 1; i < 3; i++) {
00369         by = new sdyn3();
00370         if (i_flag) by->set_label(i+1);
00371         bo->set_younger_sister(by);
00372         by->set_elder_sister(bo);
00373         by->set_parent(b);
00374 
00375         if (i == 1) {
00376             by->set_mass(1);
00377             by->set_pos(-bo->get_pos());
00378             by->set_vel(bo->get_vel());
00379         } else {
00380             by->set_mass(1);
00381             by->set_pos(0);
00382             by->set_vel(-2*bo->get_vel());
00383         }
00384         bo = by;
00385     }
00386 
00387     b->log_history(argc, argv);
00388 
00389     if (s_flag == FALSE) seed = 0;
00390     int random_seed = srandinter(seed, n_rand);
00391 
00392     cpu_init();
00393 
00394     for (int i = 0; i < n_experiments; i++) {
00395 
00396         if (n_experiments > 1) cerr << i+1 << ": ";
00397 
00398         cerr << "Random seed = " << get_initial_seed()
00399              << "  n_rand = " << get_n_rand() << flush;
00400         init.id = get_initial_seed() + get_n_rand();
00401 
00402         intermediate_state3 inter;
00403         final_state3 final;
00404 
00405         real cpu = cpu_time();  
00406         bound3(b, init, inter, final, cpu_time_check,
00407                  dt_out, dt_snap, snap_cube_size);
00408         cpu = cpu_time() - cpu;
00409 
00410         cerr << ":  ";
00411 //      print_bound3_outcome(inter, final, cerr);
00412         print_final(cerr, final);
00413     }
00414 }
00415 
00416 #endif