scatter3.C Source File

00001 
00002        //=======================================================//    _\|/_
00003       //  __  _____           ___                    ___       //      /|\
00004      //  /      |      ^     |   \  |         ^     |   \     //          _\|/_
00005     //   \__    |     / \    |___/  |        / \    |___/    //            /|\
00006    //       \   |    /___\   |  \   |       /___\   |   \   // _\|/_
00007   //     ___/   |   /     \  |   \  |____  /     \  |___/  //   /|\     
00008  //                                                       //            _\|/_
00009 //=======================================================//              /|\
00010 
00058 
00059 // Input to the function is an initial state, output is intermediate
00060 // and final states.
00061 
00062 // No reference is made to scatter profiles in this file.
00063 
00064 // Note the function scatter3 is completely deterministic.
00065 // No randomization is performed at this level.
00066 
00067 // Starlab library function.
00068 
00069 #include "scatter3.h"
00070 
00071 #ifndef TOOLBOX
00072 
00073 // init_to_sdyn3:  Convert an initial state to an sdyn3 for integration.
00074 
00075 local sdyn3 * init_to_sdyn3(initial_state3 & init, final_state3 & final)
00076 {
00077     set_kepler_tolerance(2);
00078 
00079     kepler k1;          // inner binary.
00080 
00081     real peri = 1;      // default value (unimportant unless init.ecc = 1).
00082     if (init.ecc == 1) peri = 0;
00083 
00084     make_standard_kepler(k1, 0, 1, -0.5 / init.sma, init.ecc,
00085                          peri, init.phase.mean_anomaly, 1);  // "1" ==> align
00086                                                              // with x, y axes
00087 
00088 #if 0
00089     cout << endl << "Inner binary:" << endl;
00090     k1.print_all(cout);
00091 #endif
00092 
00093     // Multiply the incoming velocity by vcrit.
00094 
00095     real m2 = init.m2;
00096     real m3 = init.m3;
00097     real mtotal = 1 + m3;
00098     real v_inf = init.v_inf * sqrt( (1 - m2) * m2 * mtotal / m3 );
00099     
00100     real energy3 = .5 * v_inf * v_inf;
00101     real ang_mom3 = init.rho * v_inf;
00102 
00103     real ecc3 = (energy3 == 0 ? 1
00104                               : sqrt( 1 + 2 * energy3 
00105                                             * pow(ang_mom3/mtotal, 2)));
00106 
00107     // Special case: if v_inf = 0, assume rho is periastron.
00108 
00109     real virial_ratio = init.rho;
00110 
00111     kepler k3;          // outer binary.
00112 
00113     // NOTE: passing a mean anomaly of 0 here will cause problems on
00114     // Linux systems if the orbit is linear (i.e. rho = 0), as this
00115     // implies r = 0, v = infinity.  Most systems will accept v = Inf
00116     // (which is actually OK, as we later set the separation to something
00117     // more reasonable).  However, Linux fails with a floating excetion.
00118     // Any nonzero value should be OK, as far as initialization is
00119     // concerned, but any given value of mean_anomaly will present
00120     // difficulties as v_inf --> 0 (because the corresponding separation
00121     // will be outside the desired value, leading to problems with the
00122     // function return_to_radius() below).  Use a negative value of
00123     // mean_anomaly here to keep us on the incoming branch, and take
00124     // special precautions later.
00125 
00126     real mean_anomaly = -0.01;
00127     make_standard_kepler(k3, 0, mtotal, energy3, ecc3, virial_ratio,
00128                          mean_anomaly);
00129 
00130     // Radius for "unperturbed" inner binary:
00131 
00132     real r_unp = (init.sma + init.r3)
00133                      * pow(init.tidal_tol_factor / mtotal, -1/3.0);
00134 
00135     if (r_unp <= k3.get_periastron()) {
00136         final.descriptor = preservation;
00137         final.sma = k1.get_semi_major_axis();
00138         final.ecc = k1.get_eccentricity();
00139         final.outer_separation = k3.get_periastron();
00140         final.escaper = 3;
00141         final.error = 0;
00142         final.time = 0;
00143         final.n_steps = 0;
00144         final.virial_ratio = k3.get_energy() * k3.get_periastron() 
00145                                                 / k3.get_total_mass() - 1;
00146 
00147         // Sufficiently large stellar radii (sum larger than the binary
00148         // periastron) will cause an immediate collision to occur, so we
00149         // must check for it explicitly here.
00150 
00151         if (k1.get_periastron() < init.r1 + init.r2) {
00152             final.descriptor = merger_escape_3;
00153             final.sma = final.ecc = -1;
00154         }
00155 
00156         return NULL;
00157     }
00158 
00159     init.r_init = max(init.r_init_min, min(init.r_init_max, r_unp));
00160 
00161     // NOTE: Can't assume that k3.separation is less than r_init.
00162 
00163     if (k3.get_separation() < init.r_init)
00164         k3.return_to_radius(init.r_init);
00165     else
00166         k3.advance_to_radius(init.r_init);
00167 
00168     k1.transform_to_time(k3.get_time());
00169     set_orientation(k3, init.phase);
00170 
00171 #if 0
00172     cout << endl << "Outer binary:" << endl << flush;
00173     k3.print_all(cout);
00174 #endif
00175 
00176     sdyn3 * b;
00177     b = set_up_dynamics(m2, m3, k1, k3);
00178 
00179 //  Set up radii:
00180 
00181     sdyn3 * bb = b->get_oldest_daughter();
00182     bb->set_radius(init.r1);
00183     bb = bb->get_younger_sister();
00184     bb->set_radius(init.r2);
00185     bb = bb->get_younger_sister();
00186     bb->set_radius(init.r3);
00187 
00188     return  b;
00189 }
00190 
00191 // check_init:  Make sure an initial state is sensible.
00192 
00193 local void check_init(initial_state3 & init)
00194 {
00195     if (init.m2 > 1) err_exit("check_init: m1 < 0");
00196     if (init.m2 < 0) err_exit("check_init: m2 < 0");
00197     if (init.m3 < 0) err_exit("check_init: m3 < 0");
00198     if (init.sma <= 0) err_exit("check_init: semi-major axis <= 0");
00199     if (init.ecc < 0) err_exit("check_init: eccentricity < 0");
00200     if (init.ecc > 1) err_exit("check_init: eccentricity > 1");
00201     if (init.r1 < 0) err_exit("check_init: r1 < 0");
00202     if (init.r2 < 0) err_exit("check_init: r2 < 0");
00203     if (init.r3 < 0) err_exit("check_init: r3 < 0");
00204     if (init.v_inf < 0) err_exit("check_init: v_inf < 0");
00205     if (init.rho < 0) err_exit("check_init: rho < 0");
00206     if (init.eta <= 0) err_exit("check_init: eta <= 0");
00207 }
00208 
00209 // scatter3:  Perform a three-body scattering, initializing from a specified
00210 //            state and returning intermediate- and final-state structures.
00211 
00212 void scatter3(initial_state3 & init,
00213               intermediate_state3 & inter,
00214               final_state3 & final,
00215               real cpu_time_check,
00216               real dt_out,         // diagnostic output interval
00217               real dt_snap,        // snapshot output interval
00218               real snap_cube_size, // limit output to particles within cube
00219               real dt_print,       // print output interval
00220               sdyn3_print_fp p)    // function to print output
00221 {
00222     check_init(init);
00223     inter.id = final.id = init.id;          // For bookkeeping...
00224 
00225     sdyn3 * b = init_to_sdyn3(init, final); // Note that all systems are
00226     initialize_bodies(inter.system);        // initialized as "non-particles"
00227     initialize_bodies(final.system);        // until explicitly set.
00228 
00229     if (!b) {
00230 
00231         // Initial separation was less than pericenter, so just set
00232         // set up the intermediate state and quit.
00233 
00234         inter.n_osc = 1;
00235         inter.n_kepler = 0;
00236         inter.r_min[0] = init.sma * (1 - init.ecc);
00237         inter.r_min[1] = inter.r_min[0];
00238         inter.r_min[2] = final.sma * (final.ecc - 1); // (Approximately)
00239         inter.r_min_min = inter.r_min[0];
00240         inter.descriptor = non_resonance;
00241 
00242         return;
00243     }
00244 
00245     // Save some initial data:
00246 
00247     sdyn3_to_system(b, init.system);
00248 
00249     real e_init = energy(b);
00250     real cpu_init = cpu_time();
00251     real cpu = cpu_init;
00252     
00253     // Evolve the system until the interaction is over or a collision occurs.
00254 
00255     inter.n_kepler = 0;
00256 
00257     do {
00258 
00259         int n_stars = 0;
00260         for_all_daughters(sdyn3, b, bb) n_stars++;
00261 
00262         tree3_evolve(b, CHECK_INTERVAL, dt_out, dt_snap, snap_cube_size,
00263                      init.eta, cpu_time_check, dt_print, p);
00264 
00265         // Check the CPU time.  Note that the printed CPU time is the
00266         // time since this routine was entered.
00267 
00268         if (cpu_time() - cpu > cpu_time_check) {
00269             if (cpu == cpu_init) cerr << endl; // (May be waiting in mid-line!)
00270 
00271             int p = cerr.precision(STD_PRECISION);
00272             cerr << "scatter3:  CPU time = " << cpu_time() - cpu_init
00273                  << "  time = " << b->get_time()
00274                  << "  n_steps = " << b->get_n_steps()
00275                  << endl;
00276             cerr << "           n_osc = " << b->get_n_ssd_osc()
00277                  << "  n_kepler = " << inter.n_kepler
00278                  << endl << flush;
00279             cpu = cpu_time();
00280             cerr.precision(p);
00281         }
00282 
00283         inter.n_osc = b->get_n_ssd_osc();       // For convenience
00284 
00285         sdyn3_to_system(b, inter.system);
00286         merge_collisions(b);
00287 
00288         // See if the number of stars has decreased.
00289 
00290         for_all_daughters(sdyn3, b, bbb) n_stars--;
00291         if (n_stars > 0) {
00292             if (dt_snap < VERY_LARGE_NUMBER
00293                 && system_in_cube(b, snap_cube_size)) {
00294                 put_node(cout, *b);
00295                 cout << flush;
00296             }
00297         }
00298 
00299     } while (!extend_or_end_scatter(b, init, &inter, &final));
00300 
00301     // Set intermediate state:
00302 
00303     inter.r_min_min = VERY_LARGE_NUMBER;
00304     inter.n_stars = 0;
00305 
00306     for_all_daughters(sdyn3, b, bb) {
00307         inter.index[inter.n_stars] = bb->get_index();
00308         inter.r_min[inter.n_stars] = sqrt(bb->get_min_nn_dr2());
00309         inter.r_min_min = min(inter.r_min_min, inter.r_min[inter.n_stars]);
00310         inter.n_stars++;
00311     }
00312 
00313     if (b->get_n_ssd_osc() <= 1)
00314         inter.descriptor = non_resonance;
00315     else {
00316         int n_no_nn_change = 0;
00317         for_all_daughters(sdyn3, b, bbb)
00318             if (bbb->get_nn_change_flag() == 0)
00319                 n_no_nn_change++;
00320         if (n_no_nn_change >= 2)
00321             inter.descriptor = hierarchical_resonance;
00322         else
00323             inter.descriptor = democratic_resonance;
00324     }
00325 
00326     // Set final state:
00327 
00328     final.time = b->get_time();
00329     final.n_steps = (int) b->get_n_steps();
00330     final.error = energy(b) - e_init; // *Absolute* error, note!
00331 
00332     sdyn3_to_system(b, final.system);
00333 
00334     //  Check for integration errors (relax tolerance in the case of mergers).
00335 
00336     if (abs(final.error) > ENERGY_TOLERANCE)
00337         if ((   final.descriptor != merger_binary_1
00338              && final.descriptor != merger_binary_2
00339              && final.descriptor != merger_binary_3
00340              && final.descriptor != merger_escape_1
00341              && final.descriptor != merger_escape_2
00342              && final.descriptor != merger_escape_3
00343              && final.descriptor != triple_merger)
00344             || abs(final.error)
00345                  > MERGER_ENERGY_TOLERANCE * abs(potential_energy(b)))
00346             final.descriptor = error;
00347 
00348     // In all cases, the final system should be a properly configured
00349     // 1-, 2-, or 3-body system.  Print it for the last time, to reflect
00350     // final state, if the "-D" command-line option was specified.
00351 
00352     if (dt_snap < VERY_LARGE_NUMBER && system_in_cube(b, snap_cube_size)) {
00353         put_node(cout, *b);
00354         cout << flush;
00355     }
00356 
00357     // Delete the 3-body sytem.
00358 
00359     sdyn3 * bi = b->get_oldest_daughter();
00360     while (bi) {
00361         sdyn3 * tmp = bi->get_younger_sister();
00362         delete bi;
00363         bi = tmp;
00364     }
00365     delete b;
00366 }
00367 
00368 #else
00369 
00370 main(int argc, char **argv)
00371 {
00372     initial_state3 init;
00373     make_standard_init(init);
00374 
00375     int  seed       = 0;        // seed for random number generator
00376     int n_rand      = 0;        // number of times to invoke the generator
00377                                 // before starting for real
00378     int  n_experiments = 1;     // default: only one run
00379     real dt_out     =           // output time interval
00380           VERY_LARGE_NUMBER;
00381     real dt_snap    =           // output time interval
00382           VERY_LARGE_NUMBER;
00383 
00384     real cpu_time_check = 3600;
00385     real snap_cube_size = 10;
00386 
00387     int planar_flag = 0;
00388     int psi_flag = 0;
00389     real psi = 0;
00390 
00391     bool  b_flag = FALSE;
00392     bool  q_flag = FALSE;
00393     bool  Q_flag = FALSE;
00394 
00395     check_help();
00396 
00397     extern char *poptarg;
00398     int c;
00399     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:";
00400 
00401     while ((c = pgetopt(argc, argv, param_string)) != -1)
00402         switch(c) {
00403 
00404             case 'A': init.eta = atof(poptarg);
00405                       break;
00406             case 'b': b_flag = 1 - b_flag;
00407                       break;
00408             case 'c': cpu_time_check = 3600*atof(poptarg);// (Specify in hours)
00409                       break;
00410             case 'C': snap_cube_size = atof(poptarg);
00411                       break;
00412             case 'd': dt_out = atof(poptarg);
00413                       if (dt_out < 0) dt_out = pow(2.0, dt_out);
00414                       break;
00415             case 'D': dt_snap = atof(poptarg);
00416                       if (dt_snap < 0) dt_snap = pow(2.0, dt_snap);
00417                       break;
00418             case 'e': init.ecc = atof(poptarg);
00419                       break;
00420             case 'g': init.tidal_tol_factor = atof(poptarg);
00421                       break;
00422             case 'L': init.r_init_min = atof(poptarg);
00423                       break;
00424             case 'm': init.m2 = atof(poptarg);
00425                       break;
00426             case 'M': init.m3 = atof(poptarg);
00427                       break;
00428             case 'n': n_experiments = atoi(poptarg);
00429                       break;
00430             case 'N': n_rand = atoi(poptarg);
00431                       break;
00432             case 'o': psi = atof(poptarg);
00433                       psi_flag = 1;
00434                       break;
00435             case 'p': planar_flag = 1;
00436                       break;
00437             case 'P': planar_flag = -1;
00438                       break;
00439             case 'q': q_flag = 1 - q_flag;
00440                       break;
00441             case 'Q': Q_flag = 1 - Q_flag;
00442                       break;
00443             case 'r': init.rho = atof(poptarg);
00444                       break;
00445             case 'R': init.r_stop = atof(poptarg);
00446                       init.r_init_min = init.r_init_max = abs(init.r_stop);
00447                       break;
00448             case 's': seed = atoi(poptarg);
00449                       break;
00450             case 'S': init.r_stop = atof(poptarg);
00451                       break;
00452             case 'U': init.r_init_max = atof(poptarg);
00453                       break;
00454             case 'v': init.v_inf = atof(poptarg);
00455                       break;
00456             case 'x': init.r1 = atof(poptarg);
00457                       break;
00458             case 'y': init.r2 = atof(poptarg);
00459                       break;
00460             case 'z': init.r3 = atof(poptarg);
00461                       break;
00462             case '?': params_to_usage(cerr, argv[0], param_string);
00463                       get_help();
00464         }            
00465 
00466     if (Q_flag) q_flag = TRUE;
00467 
00468     if (init.m2 > 1) {
00469         cerr << "scatter3: init.m2 = " << init.m2 << " > 1" << endl;
00470         exit(1);
00471     }
00472 
00473     cpu_init();
00474     int random_seed = srandinter(seed, n_rand);
00475 
00476     for (int i = 0; i < n_experiments; i++) {
00477 
00478         if (n_experiments > 1) cerr << i+1 << ": ";
00479 
00480         cerr << "Random seed = " << get_initial_seed()
00481              << "  n_rand = " << get_n_rand() << flush;
00482         init.id = get_initial_seed() + get_n_rand();
00483 
00484         // Normally, we just want random angles and phases.
00485         // However, in the planar case, it may be desireable to
00486         // control the orientation of the outer orbit.  The
00487         // angle between the inner and outer orbital axes in
00488         // the planar case is init.phase.psi.  Specify psi
00489         // in *degrees*, with psi = 0 meaning that the outer
00490         // periastron lies along the positive x-axis.
00491 
00492         randomize_angles(init.phase);
00493 
00494         if (planar_flag == 1) {
00495             init.phase.cos_theta = 1;   // Planar prograde
00496             if (psi_flag) init.phase.psi = psi * M_PI / 180.0;
00497         } else if (planar_flag == -1) {
00498             init.phase.cos_theta = -1;  // Planar retrograde
00499             if (psi_flag) init.phase.psi = psi * M_PI / 180.0;
00500         }
00501 
00502         intermediate_state3 inter;
00503         final_state3 final;
00504 
00505         real cpu = cpu_time();  
00506         scatter3(init, inter, final, cpu_time_check,
00507                  dt_out, dt_snap, snap_cube_size);
00508         cpu = cpu_time() - cpu;
00509 
00510         cerr << ":  ";
00511         print_scatter3_outcome(inter, final, cerr);
00512 
00513         if (Q_flag) print_scatter3_summary(inter, final, cpu, cerr);
00514 
00515         if (!q_flag) print_scatter3_report(init, inter, final,
00516                                            cpu, b_flag, cerr);
00517     }
00518 }
00519 
00520 #endif