d_ecc_test.C Source File

00001 
00002 // d_ecc.C: Perform a series of three-body scattering experiments,
00003 //        calculating statistics on the eccentricity change.
00004 
00005 #include "scatter3.h"
00006 
00007 static real e_init;
00008 
00009 local real energy(sdyn3* b)
00010 {
00011     real k = 0, u = 0, dissipation = 0;
00012 
00013     for_all_daughters(sdyn3, b, bi) {
00014         k += bi->get_mass() * bi->get_vel() * bi->get_vel();
00015         dissipation += bi->get_energy_dissipation();
00016         for (sdyn3 * bj = bi->get_younger_sister(); bj != NULL;
00017              bj = bj->get_younger_sister())
00018           u -= bi->get_mass() * bj->get_mass()
00019                / abs(bi->get_pos() - bj->get_pos());
00020     }
00021 
00022     return 0.5*k + u + dissipation;
00023 }
00024 
00025 // Print hook into low_n3_evolve:
00026 
00027 local void print(sdyn3* b)
00028 {
00029     kepler k;
00030 
00031     // Find the closest pair:
00032 
00033     sdyn3* b1 = NULL;
00034     sdyn3* b2 = NULL;
00035     real rmin2 = VERY_LARGE_NUMBER;
00036 
00037     for_all_daughters(sdyn3, b, bi)
00038         for (sdyn3* bj = bi->get_younger_sister(); bj != NULL;
00039              bj = bj->get_younger_sister()) {
00040             vector rij = bi->get_pos() - bj->get_pos();
00041             real rij2 = rij * rij;
00042             if (rij2 < rmin2) {
00043                 b1 = bi;
00044                 b2 = bj;
00045                 rmin2 = rij2;
00046             }
00047         }
00048     set_kepler_from_sdyn3(k, b1, b2);
00049 
00050     // Identify the third star:
00051 
00052     sdyn3* b3 = NULL;
00053     for_all_daughters(sdyn3, b, bk)
00054         if (b3 == NULL && bk != b1 && bk != b2) b3 = bk;
00055     vector cm = (b1->get_mass()*b1->get_pos() + b2->get_mass()*b2->get_pos())
00056                   / (b1->get_mass() + b2->get_mass());
00057 
00058     if (e_init == VERY_LARGE_NUMBER) e_init = energy(b);
00059 
00060     vector x1 = b1->get_pos() - cm;
00061     vector x2 = b2->get_pos() - cm;
00062 
00063     cout << b->get_time() + b->get_time_offset() << " "
00064          << abs(b3->get_pos() - cm) << " "
00065          << k.get_eccentricity() << " "
00066          << energy(b) - e_init << " "
00067          << x1 << " " << x2 << " "
00068          << endl;
00069 }
00070 
00071 #define N_OUT           10      // feedback
00072 #define N_BINS          24
00073 #define BINS_PER_DECADE  4
00074 #define MIN_STATS        5
00075 
00076 local int realcomp(const void * a, const void * b)      // For use by qsort
00077 {
00078     if (*((real *) a) > *((real *) b))
00079         return 1;
00080     else if (*((real *) a) < *((real *) b))
00081         return -1;
00082     else
00083         return 0;
00084 }
00085 
00086 // run_trials: Perform the experiments and get statistics.
00087 
00088 local void run_trials(int n_trials, initial_state3 & init,
00089                       int out_flag, int planar_flag, int dump_flag,
00090                       int debug_flag, real dt_print)
00091 {
00092     intermediate_state3 inter;
00093     final_state3 final;
00094 
00095     int total_ok = 0;
00096     real d_ecc_sum = 0, d_ecc_sq_sum = 0, d_ecc_min = 100, d_ecc_max = -1;
00097 
00098     real* de = new real[n_trials];
00099     int de_histo[N_BINS];
00100 
00101     int i;
00102     for (i = 0; i < N_BINS; i++) de_histo[i] = 0;
00103 
00104     init.r_init_min = init.r_init_max = init.r_stop;
00105 
00106     for (i = 0; i < n_trials; i++) {
00107 
00108         randomize_angles(init.phase);
00109         if (planar_flag == 1)
00110             init.phase.cos_theta = 1;   // Planar prograde
00111         else if (planar_flag == -1)
00112             init.phase.cos_theta = -1;  // Planar retrograde
00113 
00114         e_init = VERY_LARGE_NUMBER;
00115         scatter3(init, inter, final,
00116                  VERY_LARGE_NUMBER,
00117                  VERY_LARGE_NUMBER,
00118                  VERY_LARGE_NUMBER,
00119                  0,
00120                  dt_print, print);
00121 
00122         // Check for return passages:
00123 
00124         if (inter.n_kepler > 0 || inter.n_osc != 1
00125              || inter.descriptor != non_resonance)
00126             cerr << "  error at trial #" << i
00127                  << ":  n_kep = " << inter.n_kepler
00128                  << "  n_osc = " << inter.n_osc
00129                  << ":  " << state_string(inter.descriptor) << endl;
00130         else {
00131 
00132             // Gather statistics:
00133 
00134             real d_ecc = final.ecc - init.ecc;
00135             d_ecc_sum += d_ecc;
00136             d_ecc_sq_sum += d_ecc*d_ecc;
00137             d_ecc_min = min(d_ecc_min, d_ecc);
00138             d_ecc_max = max(d_ecc_max, d_ecc);
00139 
00140             if (debug_flag) {
00141                 cout << "  trial #" << i << ":  d_ecc = " << d_ecc << endl;
00142                 cout << "  minima: " << inter.r_min[0]
00143                      << "  " << inter.r_min[1]
00144                      << "  " << inter.r_min[2] << endl;
00145                 cout << "  energy = " << e_init
00146                      << "  error = " << final.error << endl;
00147             }
00148 
00149             *(de + total_ok) = d_ecc;
00150 
00151             int ie = (int) (BINS_PER_DECADE * log10(d_ecc)) + N_BINS - 1;
00152             if (ie < 0)
00153                 ie = 0;
00154             else if (ie >= N_BINS)
00155                 ie = N_BINS - 1;
00156 
00157             de_histo[ie]++;         
00158 
00159             total_ok++;
00160             if (out_flag && total_ok % N_OUT == 0) cout << total_ok << endl;
00161         }
00162     }
00163 
00164     if (total_ok >= MIN_STATS) {
00165         real mean = d_ecc_sum / total_ok;
00166         cout << "  total hits = " << total_ok << "  <d_ecc> = " << mean
00167              << "  sigma = " << sqrt(d_ecc_sq_sum / total_ok - mean * mean)
00168              << endl;
00169         cout << "  min d_ecc = " << d_ecc_min
00170              << "  max = " << d_ecc_max << endl;
00171 
00172         qsort((void*)de, (size_t)total_ok, sizeof(real), realcomp);
00173 
00174         cout.precision(4);
00175         real factor = 0.01*total_ok;
00176         cout << "  10-percentiles:" << endl << " ";
00177         for (i = 10; i <= 90; i += 10) cout << " " << *(de + (int)(factor*i));
00178         cout << endl;
00179 
00180         cout << "  histogram (" << BINS_PER_DECADE << " bins per decade):\n";
00181         int j = 0;
00182         while (j < N_BINS) {    // cout and printf don't mesh on Solaris 2.3!
00183             printf("\t");
00184             for (i = 0; i < BINS_PER_DECADE && j < N_BINS; i++, j++)
00185                 printf("  %5d", de_histo[j]);   
00186             printf("\n");
00187         }
00188 
00189         if (dump_flag) {
00190             cout << "\n  Raw data (total = " << total_ok << ", sorted):\n";
00191             j = 0;
00192             cout << "\t";
00193             while (j < total_ok) {
00194                 for (i = 0; i < 5 && j < total_ok; i++, j++)
00195                     printf("  %12.8f", *(de + j));
00196                 cout << endl;
00197             }
00198         }
00199     }
00200 }
00201 
00202 main(int argc, char **argv)
00203 {
00204     int  n_trials = 100;
00205     int  seed = 0;
00206     int  out_flag = 0;
00207     int  dump_flag = 0;
00208     int  debug_flag = 0;
00209     int  planar_flag = 0;
00210     real dt_print = VERY_LARGE_NUMBER;
00211 
00212     extern char *poptarg;
00213     int  pgetopt(int, char **, char *), c;
00214 
00215     initial_state3 init;
00216     make_standard_init(init);
00217 
00218     init.rho = 5;
00219     init.r_stop = 100;
00220     init.v_inf = 0;
00221 
00222     // Defaults: n_trials = 100
00223     //           ecc   = 0
00224     //           rho   = 5
00225     //           phase = random
00226     //           v_inf = 0 (parabola, rho = periastron)
00227     //           m1 = m2 = m3 = 0.5
00228     //           r_stop = 100
00229     //           no extra output
00230 
00231     while ((c = pgetopt(argc, argv, "A:dDe:m:M:n:opPr:R:s:t:v:")) != -1)
00232         switch(c)
00233             {
00234             case 'A': init.eta = atof(poptarg);
00235                       break;
00236             case 'd': dump_flag = 1 - dump_flag;
00237                       break;
00238             case 'D': debug_flag = 1 - debug_flag;
00239                       break;
00240             case 'e': init.ecc = atof(poptarg);
00241                       break;
00242             case 'm': init.m2 = atof(poptarg);
00243                       break;
00244             case 'M': init.m3 = atof(poptarg);
00245                       break;
00246             case 'n': n_trials = atoi(poptarg);
00247                       break;
00248             case 'o': out_flag = 1 - out_flag;
00249                       break;
00250             case 'p': planar_flag = 1;
00251                       break;
00252             case 'P': planar_flag = -1;
00253                       break;
00254             case 'r': init.rho = atof(poptarg);
00255                       break;
00256             case 'R': init.r_stop = atof(poptarg);
00257                       break;
00258             case 's': seed = atoi(poptarg);
00259                       break;
00260             case 't': dt_print = atof(poptarg);
00261                       debug_flag = 1;
00262                       break;
00263             case 'v': init.v_inf = atof(poptarg);
00264                       break;
00265             case '?': cerr << "usage: ecc [-A #] [-e #] [-m #] [-M #] [-n #] "
00266                            << "[-o] [-p] [-P] [-r #] [-R #] [-s #] [-t #] "
00267                            << "[-v #]"
00268                            << endl;
00269                       exit(1);
00270             }
00271 
00272     srandinter(seed);
00273     real peri = init.rho;
00274 
00275     // Force input rho to mean periastron in ALL cases...
00276     // Remember to scale the velocity by v_crit.
00277 
00278     if (init.v_inf > 0)
00279         init.rho = sqrt(init.rho * (init.rho + 2 * init.m3
00280                                                  / ((1 - init.m2) * init.m2
00281                                                        * pow(init.v_inf, 2))));
00282 
00283     // Basic output on the experiment:
00284 
00285     cout << "  Masses:"
00286          << "  m1 = " << 1 - init.m2
00287          << "  m2 = " << init.m2
00288          << "  m3 = " << init.m3
00289          << endl;
00290     cout << "  v_inf = " << init.v_inf
00291          << "  periastron = " << peri
00292          << "  rho = " << init.rho
00293          << "  ecc = " << init.ecc;
00294     if (planar_flag == 1)
00295         cout << "  planar prograde orbit";
00296     else if (planar_flag == -1)
00297         cout << "  planar retrograde orbit";
00298     cout << endl;
00299 
00300     cout << "  Start/stop at R > " << init.r_stop
00301          << "  accuracy parameter = " << init.eta
00302          << "  random seed = " << get_initial_seed()
00303          << endl;
00304 
00305     init.r_init_min = init.r_init_max = init.r_stop;
00306 
00307     cpu_init();
00308     run_trials(n_trials, init, out_flag, planar_flag, 
00309                dump_flag, debug_flag, dt_print);
00310 }