new_stats.C Source File

00001 
00002 // new_stats.C: Determine merger statistics for 3-body scattering.
00003 //              Steve McMillan, Charles Bailyn, Alison Procter (1995).
00004 
00005 //              This version doesn't compute cross sections.  Instead,
00006 //              it prints out raw data in a form that can be summed
00007 //              after the fact.
00008 
00009 // Starlab application:  get_sigma3.
00010 
00011 #include "sigma3.h"
00012 #define DEBUG 0
00013 
00014 #define Grav    6.673e-8
00015 #define Msun    1.989e33
00016 #define AU      1.486e13
00017 #define Rsun    6.960e10
00018 
00019 #define Kms     1e5
00020 #define SECONDS_IN_YEAR 3.16e7
00021 #define CM_IN_PC 3.086e18
00022 
00023 static real r_unit;             // Initial binary semi-major axis in A.U.
00024 static real m_unit;             // Initial total binary mass, in solar masses
00025 static real v_unit;             // Initial v_crit, in km/s
00026 
00027 static int n_merge;
00028 
00029 #define MAX_SAVE 100000
00030 
00031 static intermediate_descriptor3 save_inter[MAX_SAVE];
00032 static final_descriptor3        save_final[MAX_SAVE];
00033 static int                      save_bin[MAX_SAVE];
00034 static real                     save_ecc_0[MAX_SAVE];
00035 static real                     save_sma[MAX_SAVE];
00036 static real                     save_ecc_1[MAX_SAVE];
00037 
00038 // ------------------------------------------------------------------------
00039 
00040 static real m_lower  = 0.1;
00041 static real m_upper  = 0.8;
00042 static real exponent = 1.35;    // Convention: Salpeter = +1.35!
00043 
00044 // Interpolation from real stellar models:
00045 
00046 static int metal = 1;           // Default is solar metallicity
00047 
00048 // Low metallicity, for globular clusters:
00049 
00050 #define N_INT_0 9
00051 static real m_int_0[N_INT_0]
00052                 = {0.1, 0.2, 0.3, 0.4, 0.5,  0.6,  0.7,  0.75, 0.8};
00053 static real r_int_0[N_INT_0]
00054                 = {0.1, 0.2, 0.3, 0.4, 0.51, 0.63, 0.75, 1.09, 1.87};
00055 
00056 // Solar metallicity, for open clusters:
00057 
00058 #define N_INT_1 13
00059 static real m_int_1[N_INT_1]
00060                 = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
00061                    0.8, 0.9, 1.0, 1.1, 1.2, 1.25};
00062 static real r_int_1[N_INT_1]
00063                 = {0.140, 0.216, 0.292, 0.377, 0.465, 0.555, 0.644,
00064                    0.733, 0.836, 0.977, 1.182, 1.466, 1.767};
00065 
00066 local int get_index(real mass, real* m, int n)
00067 {
00068     for (int i = 0; i < n; i++) if (m[i] > mass) return i - 1;
00069     return n - 1;
00070 }
00071 
00072 local real get_quantity(real mass, real* m, real* q, int n)
00073 {
00074     int i = get_index(mass, m, n);
00075 
00076     if (i < 0)
00077         return q[0];
00078     else if (i >= n)
00079         return q[n-1];
00080 
00081     return q[i] + (mass - m[i]) * (q[i+1] - q[i]) / (m[i+1] - m[i]);
00082 }
00083 
00084 local real radius(real mass)    // Mass-radius relation (units: solar)
00085 {
00086     real *r, *m;
00087     int i, n;
00088 
00089     // Choose which interpolation arrays to use:
00090 
00091     if (metal == 0) {
00092         m = m_int_0;
00093         r = r_int_0;
00094         n = N_INT_0;
00095     } else {
00096         m = m_int_1;
00097         r = r_int_1;
00098         n = N_INT_1;
00099     }
00100 
00101     return get_quantity(mass, m, r, n);
00102 }
00103 
00104 // ------------------------------------------------------------------------
00105 
00106 // get_mass: Return a random mass between ml and mu, distributed as a
00107 //           power law with exponent x (i.e. dN/dm \propto m^x).
00108 
00109 local real get_mass(real ml, real mu, real x)
00110 {
00111     real r = randinter(0,1);
00112 
00113     if (x == -1) 
00114         return ml * exp( r*log(mu/ml) );
00115     else
00116         return ml * pow( 1 + r * (pow(mu/ml, 1+x) - 1), 1.0/(1+x));
00117 }
00118 
00119 // ------------------------------------------------------------------------
00120 
00121 // save_stats: Save essential raw data supplied by get_sigma3.
00122 
00123 local void save_stats(scatter_profile& prof,
00124                       initial_state3& init,
00125                       intermediate_state3& inter,
00126                       final_state3& final,
00127                       int rho_bin,
00128                       sigma_out& out)
00129 {
00130     if (   final.descriptor == merger_binary_1
00131         || final.descriptor == merger_escape_1
00132         || final.descriptor == merger_binary_2
00133         || final.descriptor == merger_escape_2
00134         || final.descriptor == merger_binary_3
00135         || final.descriptor == merger_escape_3
00136         || final.descriptor == triple_merger  ) {
00137 
00138         // Save data only in case of merger.
00139 
00140         save_ecc_0[n_merge] = init.ecc;
00141 
00142         save_inter[n_merge] = inter.descriptor;
00143         save_final[n_merge] = final.descriptor;
00144         save_bin[n_merge] = rho_bin;
00145 
00146         if (   final.descriptor == merger_binary_1
00147             || final.descriptor == merger_binary_2
00148             || final.descriptor == merger_binary_3 ) {
00149             save_sma[n_merge] = final.sma;
00150             save_ecc_1[n_merge] = final.ecc;
00151         }
00152 
00153         n_merge++;
00154     }
00155 }
00156 
00157 // dump_stats: print the data saved by save_stats.
00158 
00159 void dump_stats(scatter_profile& prof, sigma_out& out, int verbose)
00160 {
00161     if (n_merge > 0) {
00162 
00163         int i;
00164 
00165         // Header information:
00166 
00167         if (verbose) {
00168 
00169             int p = cerr.precision(4);          // nonstandard precision
00170             cerr << "\n    max_zone = " << out.i_max
00171                  << "  rho_max = " << out.rho_max
00172                  << endl;
00173 
00174             cerr << "    trials per zone = " << out.trials_per_zone;
00175             cerr << ", total hits per zone:";
00176             for (i = 0; i <= out.i_max; i++) cerr << " " << out.n_hit[i];
00177             cerr << endl;
00178 
00179             cerr << "    (v/v_c)^2 * zone weights:";
00180             for (i = 0; i <= out.i_max; i++)
00181                 cerr << " " << prof.v_inf*prof.v_inf * zone_weight(out, i);
00182             cerr << endl;
00183 
00184             cerr << endl <<
00185          "  total merging cross sections [(v/v_c)^2 * sigma / (pi a^2)]:\n\n";
00186             print_sigma3_mergers(out.sigma, out.sigma_err_sq,
00187                              prof.v_inf*prof.v_inf);
00188             cerr.precision(p);
00189 
00190         } else {
00191 
00192             fprintf(stderr, "\n%s%s%s%s%s%s%s%s%s%s%s\n",
00193                     "    m1 ",
00194                     "    m2 ",
00195                     "    m3 ",
00196                     "     a_i ",
00197                     "   e_i ",
00198                     " bin",
00199                     " int",
00200                     " fin",
00201                     "  d(sigma)",
00202                     "    a_f ",
00203                     "   e_f");
00204 
00205         }
00206 
00207         // Details, merger by merger:
00208 
00209         int p = cerr.precision(STD_PRECISION);
00210 
00211         for (i = 0; i < n_merge; i++) {
00212 
00213             if (verbose) {
00214 
00215                 if (i > 0) cerr << endl;
00216                 cerr << "  merger #" << i << endl;
00217 
00218                 cerr << "    initial a = "
00219                      << r_unit << "  e = " << save_ecc_0[i];
00220                 cerr << "   masses: " << m_unit*(1 - prof.m2)
00221                                       << " " <<  m_unit*prof.m2
00222                                       << " " <<  m_unit*prof.m3
00223                      << endl;
00224 
00225                 cerr << "    outcome: " << state_string(save_inter[i])
00226                      << " " << state_string(save_final[i])
00227                      << endl;
00228 
00229                 cerr << "    rho_bin = " << save_bin[i]
00230                      << "  d(sigma) [AU^2] = "
00231                      << zone_weight(out, save_bin[i])*PI*r_unit*r_unit
00232                      << endl;
00233 
00234                 if (   save_final[i] == merger_binary_1
00235                     || save_final[i] == merger_binary_2
00236                     || save_final[i] == merger_binary_3 )
00237                     cerr << "    final binary parameters:  a = "
00238                          << r_unit*save_sma[i]
00239                          << "  e = " << save_ecc_1[i] << endl;
00240 
00241             } else {
00242 
00243                 fprintf(stderr,
00244                         "%7.3f%7.3f%7.3f%9.3f%7.3f %3d %3d %3d%9.3f",
00245                         m_unit*(1 - prof.m2), m_unit*prof.m2, m_unit*prof.m3,
00246                         r_unit, save_ecc_0[i], save_bin[i],
00247                         save_inter[i], save_final[i],
00248                         zone_weight(out, save_bin[i])*PI*r_unit*r_unit);
00249                 if (   save_final[i] == merger_binary_1
00250                     || save_final[i] == merger_binary_2
00251                     || save_final[i] == merger_binary_3 )
00252                     fprintf(stderr, "%9.3f%7.3f",
00253                             r_unit*save_sma[i], save_ecc_1[i]);
00254                 fprintf(stderr, "\n");
00255 
00256             }
00257         }
00258         cerr.precision(p);
00259     }
00260 }
00261 
00262 // ------------------------------------------------------------------------
00263 
00264 main(int argc, char **argv)
00265 {
00266     int  verbose  = 1;
00267 
00268     int  debug  = 0;
00269     int  seed   = 0;
00270     int  n_rand = 0;
00271     real max_trial_density = 1.0;
00272 
00273     real cpu_time_check = 3600; // One check per CPU hour!
00274     real dt_snap = VERY_LARGE_NUMBER;
00275     real snap_cube_size = 10;
00276     int scatter_summary_level = 0;
00277 
00278     real sma = 1.0;     // Initial binary semi-major axis, in A.U.
00279     real v_inf = 10.0;  // Relative velocity at infinity in km/s
00280 
00281     int n_sigma3 = 1;
00282     bool sig_print = FALSE;
00283     bool stat_print = 0;
00284 
00285     scatter_profile prof;
00286     make_standard_profile(prof);
00287 
00288     extern char *poptarg;
00289     int c;
00290     char* param_string = "a:A:c:C:D:d:e:El:L:Mn:N:p:PqQs:u:U:v:Vx:";
00291 
00292     while ((c = pgetopt(argc, argv, param_string)) != -1)
00293         switch(c)
00294             {
00295             case 'a': sma = atof(poptarg);
00296                       break;
00297             case 'A': prof.eta = atof(poptarg);
00298                       break;
00299             case 'c': cpu_time_check = 3600*atof(poptarg);// (Specify in hours)
00300                       break;
00301             case 'C': snap_cube_size = atof(poptarg);
00302                       break;
00303             case 'd': max_trial_density = atof(poptarg);  // Max density for a
00304                       break;                              // single experiment
00305             case 'D': dt_snap = atof(poptarg);
00306                       scatter_summary_level = 2;  // Undo with later "-q/Q"
00307                       break;
00308             case 'e': prof.ecc = atof(poptarg);   // Specify an eccentricity
00309                       prof.ecc_flag = 1;
00310                       break;
00311             case 'E': prof.ecc_flag = 2;          // Gaussian eccentricity
00312                       break;                      // distribution
00313             case 'l':
00314             case 'L': m_lower = atof(poptarg);
00315                       break;
00316             case 'M': metal = 1 - metal;
00317                       break;
00318             case 'n': n_sigma3 = atoi(poptarg);
00319                       break;
00320             case 'N': n_rand = atoi(poptarg);
00321                       break;
00322             case 'p': stat_print = atoi(poptarg);
00323                       if (stat_print <= 0) stat_print = 1;
00324                       break;
00325             case 'P': sig_print = 1 - sig_print;
00326                       stat_print = 1;
00327                       break;
00328             case 'q': if (scatter_summary_level > 0)
00329                           scatter_summary_level = 0;
00330                       else
00331                           scatter_summary_level = 1;
00332                       break;
00333             case 'Q': if (scatter_summary_level > 0)
00334                           scatter_summary_level = 0;
00335                       else
00336                           scatter_summary_level = 2;
00337                       break;
00338             case 's': seed = atoi(poptarg);
00339                       break;
00340             case 'u':
00341             case 'U': m_upper = atof(poptarg);
00342                       break;
00343             case 'v': v_inf = atof(poptarg);
00344                       break;
00345             case 'V': verbose = 1 - verbose;
00346                       break;
00347             case 'x': exponent = atof(poptarg);
00348                       break;
00349             case '?': params_to_usage(cerr, argv[0], param_string);
00350                       exit(1);
00351             }
00352 
00353     cpu_init();
00354     int first_seed = srandinter(seed, n_rand);
00355 
00356     cerr << "\n\t*** "
00357          << "Binary/single-star scattering with stellar collisions."
00358          << " ***\n\n";
00359     cerr << "Binary semi-major-axis a = " << sma
00360          << " AU, velocity at infinity = " << v_inf << " km/s\n";
00361     cerr << "Mass range: " << m_lower << " - " << m_upper
00362          << " solar, power-law exponent = " << exponent
00363          << " (Salpeter = 1.35)\n";
00364     cerr << "Target trial density = " << max_trial_density
00365          << ", total number of mass combinations = " << n_sigma3 << endl;
00366     cerr << "Random seed = " << first_seed << endl;
00367 
00368     sigma_out out;
00369 
00370     // Loop over mass combinations, obtaining cross-sections for each.
00371 
00372     for (int i = 0; i < n_sigma3; i++) {
00373 
00374         // Get masses in solar masses, randomly from the mass function.
00375 
00376         real m1 = get_mass(m_lower, m_upper, -exponent - 1);
00377         real m2 = get_mass(m_lower, m_upper, -exponent - 1);
00378         real m3 = get_mass(m_lower, m_upper, -exponent - 1);
00379 
00380         // Scale masses, v_inf and radii:
00381 
00382         m_unit = m1 + m2;
00383         r_unit = sma;
00384         v_unit = sqrt( (Grav * Msun / AU)
00385                           * m1 * m2 * (m1 + m2 + m3) / ((m1 + m2) * m3)
00386                           / sma ) / Kms;
00387 
00388         prof.m2 = m2 / m_unit;
00389         prof.m3 = m3 / m_unit;
00390         prof.v_inf = v_inf / v_unit;
00391 
00392         prof.r1 = radius(m1) * Rsun / (sma * AU);
00393         prof.r2 = radius(m2) * Rsun / (sma * AU);
00394         prof.r3 = radius(m3) * Rsun / (sma * AU);
00395 
00396         if (verbose) {
00397 
00398             cerr << endl;
00399             for (int j = 0; j < 75; j++) cerr << '-';
00400             cerr << endl;
00401 
00402             int p = cerr.precision(STD_PRECISION);
00403             cerr << "\nMass combination #" << i+1 << ":\n";
00404             cerr << "    m1 = " << m1 << "  m2 = " << m2
00405                 << "  m3 = " << m3 << " solar\n";
00406             cerr << "    r1 = " << radius(m1) << "  r2 = " << radius(m2)
00407                 << "  r3 = " << radius(m3) << " solar\n";
00408             cerr << "    a = " << sma << " A.U. = "
00409                 << sma * AU / Rsun << " solar,";
00410             cerr << "  v_inf = " << v_inf << " km/s = "
00411                 << prof.v_inf << " scaled\n";
00412             cerr.precision(p);
00413         }
00414 
00415         // Note that the only use of get_sigma3 here is to ensure
00416         // that properly sampled data are sent to routine dump_stats.
00417 
00418         n_merge = 0;
00419 
00420         get_sigma3(max_trial_density, prof, out,
00421                    debug, cpu_time_check,
00422                    dt_snap, snap_cube_size,
00423                    scatter_summary_level,
00424                    save_stats);
00425 
00426         dump_stats(prof, out, verbose);
00427     }
00428 }