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