flyby3.C Source File

00001 
00002 // flyby3.C: Perform a series of three-body scattering experiments,
00003 //           calculating statistics on changes in eccentricity and
00004 //           semi-major axis.
00005 
00006 // Starlab application:  scatter3.
00007 
00008 #include "scatter3.h"
00009 
00010 static real e_init;
00011 
00012 local real energy(sdyn3* b)
00013 {
00014     real k = 0, u = 0, dissipation = 0;
00015 
00016     for_all_daughters(sdyn3, b, bi) {
00017         k += bi->get_mass() * bi->get_vel() * bi->get_vel();
00018         dissipation += bi->get_energy_dissipation();
00019         for (sdyn3 * bj = bi->get_younger_sister(); bj != NULL;
00020              bj = bj->get_younger_sister())
00021           u -= bi->get_mass() * bj->get_mass()
00022                / abs(bi->get_pos() - bj->get_pos());
00023     }
00024 
00025     return 0.5*k + u + dissipation;
00026 }
00027 
00028 // Print hook into low_n3_evolve:
00029 
00030 local void print(sdyn3* b)
00031 {
00032     kepler k;
00033 
00034     // Find the closest pair:
00035 
00036     sdyn3* b1 = NULL;
00037     sdyn3* b2 = NULL;
00038     real rmin2 = VERY_LARGE_NUMBER;
00039 
00040     for_all_daughters(sdyn3, b, bi)
00041         for (sdyn3* bj = bi->get_younger_sister(); bj != NULL;
00042              bj = bj->get_younger_sister()) {
00043             vector rij = bi->get_pos() - bj->get_pos();
00044             real rij2 = rij * rij;
00045             if (rij2 < rmin2) {
00046                 b1 = bi;
00047                 b2 = bj;
00048                 rmin2 = rij2;
00049             }
00050         }
00051     set_kepler_from_sdyn3(k, b1, b2);
00052 
00053     // Identify the third star:
00054 
00055     sdyn3* b3 = NULL;
00056     for_all_daughters(sdyn3, b, bk)
00057         if (b3 == NULL && bk != b1 && bk != b2) b3 = bk;
00058     vector cm = (b1->get_mass()*b1->get_pos() + b2->get_mass()*b2->get_pos())
00059                   / (b1->get_mass() + b2->get_mass());
00060 
00061     if (e_init == VERY_LARGE_NUMBER) e_init = energy(b);
00062     cout << "  t = "  << b->get_time() + b->get_time_offset()
00063          << "  R = "  << abs(b3->get_pos() - cm)
00064          << "  a = "  << k.get_semi_major_axis()
00065          << "  e = "  << k.get_eccentricity()
00066          << "  dE = " << energy(b) - e_init
00067          << endl;
00068 }
00069 
00070 #define N_OUT           10      // feedback
00071 #define N_BINS          24
00072 #define BINS_PER_DECADE  2
00073 #define MIN_STATS        5
00074 
00075 local int realcomp(const void * a, const void * b)      // For use by qsort
00076 {
00077     if (*((real *) a) > *((real *) b))
00078         return 1;
00079     else if (*((real *) a) < *((real *) b))
00080         return -1;
00081     else
00082         return 0;
00083 }
00084 
00085 local void error_print(ostream& s,
00086                        int i,
00087                        intermediate_state3& inter,
00088                        final_state3& final)
00089 
00090 {
00091     s << "  error at trial #" << i
00092       << ":  n_kepler = " << inter.n_kepler
00093       << "  n_osc = " << inter.n_osc
00094       << "  t = " << final.time
00095       << "  n_steps = " << final.n_steps
00096       << ":\n        " << state_string(inter.descriptor)
00097       << "  "  << state_string(final.descriptor)
00098       << endl;
00099 }
00100 
00101 // run_trials: Perform the experiments and get statistics.
00102 
00103 local void run_trials(int n_trials, initial_state3 & init,
00104                       int out_flag, int planar_flag, int dump_flag,
00105                       int debug_flag, real dt_print)
00106 {
00107     intermediate_state3 inter;
00108     final_state3 final;
00109 
00110     int total_ok = 0;
00111     real d_ecc_sum = 0, d_ecc_sq_sum = 0, d_ecc_min = 100, d_ecc_max = -1;
00112     real d_sma_sum = 0, d_sma_sq_sum = 0, d_sma_min = 100, d_sma_max = -1;
00113 
00114     real* d_ecc_list = new real[n_trials];
00115     real* d_sma_list = new real[n_trials];
00116     int d_ecc_histo[N_BINS];
00117     int d_sma_histo[N_BINS];
00118 
00119     real max_error = -VERY_LARGE_NUMBER, min_error = VERY_LARGE_NUMBER;
00120     real error_sum = 0, abs_error_sum = 0, error_sq_sum = 0;
00121 
00122     int i;
00123     for (i = 0; i < N_BINS; i++) d_ecc_histo[i] = d_sma_histo[i] = 0;
00124 
00125     init.r_init_min = init.r_init_max = abs(init.r_stop);
00126 
00127     if (debug_flag) cout << endl;
00128 
00129     for (i = 0; i < n_trials; i++) {
00130 
00131         randomize_angles(init.phase);
00132         if (planar_flag == 1)
00133             init.phase.cos_theta = 1;   // Planar prograde
00134         else if (planar_flag == -1)
00135             init.phase.cos_theta = -1;  // Planar retrograde
00136 
00137         e_init = VERY_LARGE_NUMBER;
00138         scatter3(init, inter, final,
00139                  VERY_LARGE_NUMBER,
00140                  VERY_LARGE_NUMBER,
00141                  VERY_LARGE_NUMBER,
00142                  0,
00143                  dt_print, print);
00144 
00145         // On return, the final state should be "stopped" and the
00146         // intermediate state should be "non_resonance".  Since kepler
00147         // extension is turned off when r_stop < 0 is specified, the 
00148         // only possibilities are (for outer separation sep and relative
00149         // velocity vr):
00150         //
00151         //      1. sep > r_stop with vr > 0 during direct integration
00152         //              ==> n_osc = 1, n_kepler = 0,
00153         //                  final.time is when the criterion was met.
00154         //      2. apocenter occurred during direct integration
00155         //              ==> n_osc = 1, n_kepler = 0, sep < r_stop, vr < 0,
00156         //                  final.time is when vr < 0 was first detected,
00157         //                  (within ~20 time units of apocenter).
00158         //
00159         // (Note that the second case guarantees d_sma > 0.)
00160 
00161 /*
00162         cout << "inter:  n_osc = " << inter.n_osc
00163              << "  n_kepler = " << inter.n_kepler
00164              << "  " << state_string(inter.descriptor) << endl;
00165         cout << "final:  sep = " << final.outer_separation
00166              << "  t = " << final.time << "  n_steps = " << final.n_steps
00167              << "  " << state_string(final.descriptor) << endl;
00168 */
00169 
00170         // Check for errors:
00171 
00172         if (inter.n_kepler > 0
00173             || inter.n_osc != 1
00174             || final.time <= 0
00175             || inter.descriptor != non_resonance
00176             || final.descriptor != stopped) {
00177             
00178             if (debug_flag) 
00179                 error_print(cout, i, inter, final);
00180             else
00181                 error_print(cerr, i, inter, final);
00182 
00183         } else {
00184 
00185             // Gather statistics:
00186 
00187             real d_ecc = final.ecc - init.ecc;
00188             d_ecc_sum += d_ecc;
00189             d_ecc_sq_sum += d_ecc*d_ecc;
00190             d_ecc_min = min(d_ecc_min, d_ecc);
00191             d_ecc_max = max(d_ecc_max, d_ecc);
00192 
00193             // Eccentricity details:
00194 
00195             *(d_ecc_list + total_ok) = d_ecc;
00196 
00197             int ie = (int) (BINS_PER_DECADE * log10(max(1.e-30,abs(d_ecc))))
00198                         + N_BINS/2 - 1;
00199             if (ie < 0)
00200                 ie = 0;
00201             else if (ie >= N_BINS/2)
00202                 ie = N_BINS/2 - 1;
00203 
00204             if (d_ecc > 0) 
00205                 ie += N_BINS/2;
00206             else
00207                 ie = N_BINS/2 -1 - ie;
00208 
00209             d_ecc_histo[ie]++;      
00210 
00211             real d_sma = final.sma - init.sma;
00212             d_sma_sum += d_sma;
00213             d_sma_sq_sum += d_sma*d_sma;
00214             d_sma_min = min(d_sma_min, d_sma);
00215             d_sma_max = max(d_sma_max, d_sma);
00216 
00217             // Semi-major axis details:
00218 
00219             *(d_sma_list + total_ok) = d_sma;
00220 
00221             int ia = (int) (BINS_PER_DECADE * log10(max(1.e-30,abs(d_sma))))
00222                         + N_BINS/2 - 1;
00223             if (ia < 0)
00224                 ia = 0;
00225             else if (ia >= N_BINS/2)
00226                 ia = N_BINS/2 - 1;
00227 
00228             if (d_sma > 0) 
00229                 ia += N_BINS/2;
00230             else
00231                 ia = N_BINS/2 -1 - ia;
00232 
00233             d_sma_histo[ia]++;      
00234 
00235             if (debug_flag) {
00236                 cout << "  trial #" << i
00237                      << ":  d_ecc, d_sma = " << d_ecc << " " << d_sma
00238                      << "  bins " << ie << " " << ia << endl;
00239                 cout << "  dr_minima: " << inter.r_min[0]
00240                      << "  " << inter.r_min[1]
00241                      << "  " << inter.r_min[2] << endl;
00242                 cout << "  energy";
00243                 if (e_init < VERY_LARGE_NUMBER) cout << " = " << e_init << " ";
00244                 cout << " error = " << final.error << endl;
00245             }
00246 
00247             if (final.error > max_error) max_error = final.error;
00248             if (final.error < min_error) min_error = final.error;
00249 
00250             error_sum += final.error;
00251             abs_error_sum += abs(final.error);
00252             error_sq_sum += final.error * final.error;
00253 
00254             total_ok++;
00255             if (out_flag && total_ok % N_OUT == 0) cout << total_ok << endl;
00256         }
00257     }
00258 
00259     if (total_ok >= MIN_STATS) {
00260 
00261         real mean_error = error_sum / total_ok;
00262         real mean_abs_error = abs_error_sum / total_ok;
00263         cout << "\n  total hits = " << total_ok
00264              << "  error range = " << min_error << " to " << max_error
00265              << "\n                    mean  error  = " << mean_error
00266              << "\n                    mean |error| = " << mean_abs_error
00267              << "  sigma = "
00268              << sqrt(error_sq_sum / total_ok - mean_error * mean_error)
00269              << "\n\n";
00270 
00271         real mean_ecc = d_ecc_sum / total_ok;
00272         cout << "  <d_ecc> = " << mean_ecc << "  sigma = "
00273              << sqrt(d_ecc_sq_sum / total_ok - mean_ecc * mean_ecc) << endl;
00274         cout << "  min d_ecc = " << d_ecc_min
00275              << "  max = " << d_ecc_max << endl;
00276 
00277         real mean_sma = d_sma_sum / total_ok;
00278         cout << "\n  <d_sma> = " << mean_sma << "  sigma = "
00279              <<  sqrt(d_sma_sq_sum / total_ok - mean_sma * mean_sma) << endl;
00280         cout << "  min d_sma = " << d_sma_min
00281              << "  max = " << d_sma_max << endl;
00282 
00283         cout << "\nEccentricity details (d_ecc):\n";
00284         cout.precision(4);
00285 
00286         cout << "  10-percentiles:" << endl << " ";
00287         qsort((void*)d_ecc_list, (size_t)total_ok, sizeof(real), realcomp);
00288         real factor = 0.01*total_ok;
00289         for (i = 10; i <= 90; i += 10)
00290             cout << " " << *(d_ecc_list + (int)(factor*i));
00291         cout << endl;
00292 
00293         cout << "  histogram (" << BINS_PER_DECADE << " bins per decade):\n";
00294         int j = 0;
00295         while (j < N_BINS) {    // cout and printf don't mesh on Solaris 2.3!
00296             printf("\t");
00297             for (i = 0; i < BINS_PER_DECADE && j < N_BINS; i++, j++)
00298                 printf("  %5d", d_ecc_histo[j]);        
00299             printf("\n");
00300         }
00301 
00302         cout << "\nSemi-major axis details (d_sma):\n";
00303         cout.precision(4);
00304 
00305         cout << "  10-percentiles:" << endl << " ";
00306         qsort((void*)d_sma_list, (size_t)total_ok, sizeof(real), realcomp);
00307         for (i = 10; i <= 90; i += 10)
00308             cout << " " << *(d_sma_list + (int)(factor*i));
00309         cout << endl;
00310 
00311         cout << "  histogram (" << BINS_PER_DECADE << " bins per decade):\n";
00312         j = 0;
00313         while (j < N_BINS) {
00314             printf("\t");
00315             for (i = 0; i < BINS_PER_DECADE && j < N_BINS; i++, j++)
00316                 printf("  %5d", d_sma_histo[j]);        
00317             printf("\n");
00318         }
00319 
00320         if (dump_flag) {
00321             cout << "\n  Raw data (total = " << total_ok << ", sorted):\n";
00322 
00323             j = 0;
00324             cout << "\n  Eccentricity:\t";
00325             while (j < total_ok) {
00326                 for (i = 0; i < 5 && j < total_ok; i++, j++)
00327                     printf("  %12.8f", *(d_ecc_list + j));
00328                 cout << endl;
00329             }
00330 
00331             j = 0;
00332             cout << "\n  Semi-major axis:\t";
00333             while (j < total_ok) {
00334                 for (i = 0; i < 5 && j < total_ok; i++, j++)
00335                     printf("  %12.8f", *(d_sma_list + j));
00336                 cout << endl;
00337             }
00338         }
00339     }
00340 }
00341 
00342 main(int argc, char **argv)
00343 {
00344     int  n_trials = 100;
00345     int  seed = 0;
00346     int  out_flag = 0;
00347     int  dump_flag = 0;
00348     int  debug_flag = 0;
00349     int  planar_flag = 0;
00350     real dt_print = VERY_LARGE_NUMBER;
00351 
00352     extern char *poptarg;
00353     int c;
00354     char* param_string = "A:dDe:m:M:n:opPr:R:s:t:v:";
00355 
00356     initial_state3 init;
00357     make_standard_init(init);
00358 
00359     init.rho = 5;
00360     init.r_stop = 100;
00361     init.v_inf = 0;
00362 
00363     // Defaults: n_trials = 100
00364     //           ecc   = 0
00365     //           rho   = 5
00366     //           phase = random
00367     //           v_inf = 0 (parabola, rho = periastron)
00368     //           m1 = m2 = m3 = 0.5
00369     //           r_stop = 100
00370     //           no extra output
00371 
00372     while ((c = pgetopt(argc, argv, param_string)) != -1)
00373         switch(c)
00374             {
00375             case 'A': init.eta = atof(poptarg);
00376                       break;
00377             case 'd': dump_flag = 1 - dump_flag;
00378                       break;
00379             case 'D': debug_flag = 1 - debug_flag;
00380                       break;
00381             case 'e': init.ecc = atof(poptarg);
00382                       break;
00383             case 'm': init.m2 = atof(poptarg);
00384                       break;
00385             case 'M': init.m3 = atof(poptarg);
00386                       break;
00387             case 'n': n_trials = atoi(poptarg);
00388                       break;
00389             case 'o': out_flag = 1 - out_flag;
00390                       break;
00391             case 'p': planar_flag = 1;
00392                       break;
00393             case 'P': planar_flag = -1;
00394                       break;
00395             case 'r': init.rho = atof(poptarg);
00396                       break;
00397             case 'R': init.r_stop = -atof(poptarg); // Kludge: stop at apocenter
00398                       break;
00399             case 's': seed = atoi(poptarg);
00400                       break;
00401             case 't': dt_print = atof(poptarg);
00402                       debug_flag = 1;
00403                       break;
00404             case 'v': init.v_inf = atof(poptarg);
00405                       break;
00406             case '?': params_to_usage(cerr, argv[0], param_string);
00407                       exit(1);
00408             }
00409 
00410     srandinter(seed);
00411     real peri = init.rho;
00412 
00413     // Force input rho to mean periastron in ALL cases...
00414     // Remember to scale the velocity by v_crit.
00415 
00416     if (init.v_inf > 0)
00417         init.rho = sqrt(init.rho * (init.rho + 2 * init.m3
00418                                                  / ((1 - init.m2) * init.m2
00419                                                        * pow(init.v_inf, 2))));
00420 
00421     // Basic output on the experiment:
00422 
00423     cout << "  Masses:"
00424          << "  m1 = " << 1 - init.m2
00425          << "  m2 = " << init.m2
00426          << "  m3 = " << init.m3
00427          << endl;
00428     cout << "  v_inf = " << init.v_inf
00429          << "  periastron = " << peri
00430          << "  rho = " << init.rho
00431          << "  ecc = " << init.ecc;
00432     if (planar_flag == 1)
00433         cout << "  planar prograde orbit";
00434     else if (planar_flag == -1)
00435         cout << "  planar retrograde orbit";
00436     cout << endl;
00437 
00438     cout << "  Start/stop at R > " << abs(init.r_stop)
00439          << "  accuracy parameter = " << init.eta
00440          << "  random seed = " << get_initial_seed()
00441          << endl;
00442 
00443     cpu_init();
00444     run_trials(n_trials, init, out_flag, planar_flag, 
00445                dump_flag, debug_flag, dt_print);
00446 }