rate3.C Source File

00001 
00036 
00037 // Starlab library function, with driver application.
00038 
00039 #include "sigma3.h"
00040 
00041 #ifdef TOOLBOX
00042 
00043 #define Grav    6.673e-8
00044 #define Msun    1.989e33
00045 #define Rsun    6.960e10
00046 #define Kms     1e5
00047 #define SECONDS_IN_YEAR 3.16e7
00048 #define CM_IN_PC 3.086e18
00049 
00050 // summarize_sigma: Output enough information to restart the current get_sigma.
00051 
00052 local void summarize_sigma(real max_dens, scatter_profile & prof,
00053                            real cpu_time_check,
00054                            real dt_snap, real snap_cube_size,
00055                            int scatter_summary_level)
00056 {
00057     int p = cerr.precision(STD_PRECISION);
00058     cerr << "\nsigma3"
00059          << " -s " << get_initial_seed()
00060          << " -N " << get_n_rand()
00061          << " -d " << max_dens;
00062 
00063     cerr.precision(16);                         // Uglify the output
00064     cerr << " -m " << prof.m2
00065          << " -M " << prof.m3
00066          << " -v " << prof.v_inf;
00067     if (prof.ecc_flag) cerr << " -e " << prof.ecc;
00068     cerr << " -x " << prof.r1
00069          << " -y " << prof.r2
00070          << " -z " << prof.r3
00071          << " -A " << prof.eta
00072          << " -g " << prof.tidal_tol_factor;
00073 
00074     cerr.precision(STD_PRECISION);
00075     cerr << " -c " << cpu_time_check/3600;      // (specified in hours!)
00076 
00077     if (dt_snap < VERY_LARGE_NUMBER) {
00078         cerr << " -D " << dt_snap
00079              << " -C " << snap_cube_size;
00080     }
00081 
00082     if (scatter_summary_level == 1)
00083         cerr << " -q ";
00084     else if (scatter_summary_level == 2)
00085         cerr << " -Q ";
00086 
00087     cerr << endl;
00088     cerr.precision(p);
00089 }
00090 
00091 local real stat_weight(real v, real v_th_3d) // $v^3 f(v)$
00092 {
00093     // Return the normalized Maxwellian distribution function at velocity v.
00094 
00095     real x = v/v_th_3d;
00096     return 3 * sqrt(6/PI) * x * x * x * exp(-1.5*x*x);
00097 
00098     // Normalization is such that $\int v^2 f(v) dv = 1$
00099     //                        and $\int v^4 f(v) dv = v_{th,3d}^2$
00100 }
00101 
00102 local int get_rate3(scatter_profile& prof, real v_rel_th,
00103                     int nv, real max_dens,
00104                     real total_sigma[][N_FINAL], real total_err_sq[][N_FINAL],
00105                     real cpu_time_check, real cpu_init, real& cpu_save,
00106                     real dt_snap, real snap_cube_size,
00107                     int scatter_summary_level, int sigma_summary_flag)
00108 {
00109     real v_max = 2 * v_rel_th;
00110     real dv = v_max / nv;
00111 
00112     int ii, jj;
00113 
00114     // Initialize counters:
00115 
00116     int total_trials = 0;
00117 
00118     for (ii = 0; ii < N_INTER; ii++)
00119         for (jj = 0; jj < N_FINAL; jj++) {
00120             total_sigma[ii][jj] = 0;
00121             total_err_sq[ii][jj] = 0;
00122         }
00123 
00124     // Trapezoid rule in v_inf (skip end-points):
00125 
00126     for (int i = 1; i < nv; i++) {
00127 
00128         prof.v_inf = i * dv;
00129 
00130         real cpu_temp = cpu_time();
00131 
00132         if (sigma_summary_flag)
00133             summarize_sigma(max_dens, prof,
00134                             cpu_time_check,
00135                             dt_snap, snap_cube_size,
00136                             scatter_summary_level);
00137 
00138         sigma_out out;
00139         get_sigma3(max_dens, prof, out,
00140                    0, cpu_time_check, dt_snap, snap_cube_size,
00141                    scatter_summary_level);
00142 
00143         int p = cerr.precision(STD_PRECISION);
00144         if (sigma_summary_flag) {
00145             cerr << "\nv_inf = " << prof.v_inf << ":\n";
00146             print_sigma3(out, prof.v_inf * prof.v_inf);
00147         } else
00148             cerr << "\nv_inf = " << prof.v_inf
00149                  << ",  total_trials = " << out.total_trials
00150                  << ",  delta CPU time = " << cpu_time() - cpu_temp << " s"
00151                  << endl;
00152         cerr.precision(p);
00153 
00154         total_trials += out.total_trials;
00155         real weight = stat_weight(prof.v_inf, v_rel_th) * dv / v_rel_th;
00156 
00157         for (ii = 0; ii < N_INTER; ii++)
00158             for (jj = 0; jj < N_FINAL; jj++) {
00159 
00160                 total_sigma[ii][jj] += out.sigma[ii][jj] * weight;
00161 
00162                 // Just add errors in quadrature (sqrt taken at end):
00163 
00164                 total_err_sq[ii][jj] += out.sigma_err_sq[ii][jj]
00165                                         * weight * weight;
00166                 }
00167 
00168         // (Unit of total here is pi a^2 v_rel_th.)
00169 
00170         // Check the CPU time.  Note that the printed CPU time is the
00171         // time since this routine was entered.
00172 
00173         if (cpu_time() - cpu_save > cpu_time_check) {
00174             cpu_save = cpu_time();
00175             cerr << "\nget_rate3:  CPU time = " << cpu_save - cpu_init
00176                  << "  nv = " << nv << "  i = " << i
00177                  << "  v_inf = " << prof.v_inf
00178                  << "  total_trials = " << total_trials
00179                  << endl << flush;
00180         }
00181     }
00182 
00183     return total_trials;
00184 }
00185 
00186 local int rescale_physical(real sma, real v_rel_th, real v_unit,
00187                            real total_sigma[][N_FINAL],
00188                            real total_err_sq[][N_FINAL])
00189 {
00190     real scale_factor = PI * sma * sma * v_rel_th               // to code units
00191                            * Rsun * Rsun * v_unit * Kms         // to cm^3 / s
00192                            * SECONDS_IN_YEAR * 1e9              // to cm^3 / Gyr
00193                            / pow(CM_IN_PC, 3);                  // to pc^3 / Gyr
00194 
00195     real total_max = -1;
00196     real total_max_err = -1;
00197     for (int ii = 0; ii < N_INTER; ii++)
00198         for (int jj = 0; jj < N_FINAL; jj++) {
00199             total_sigma[ii][jj] *= scale_factor;
00200             total_err_sq[ii][jj] *= scale_factor*scale_factor;
00201             if (ii + jj > 0) {
00202                 // (Don't include the top-left element in the array...)
00203                 total_max = max(total_max, total_sigma[ii][jj]);
00204                 total_max_err = max(total_max_err, sqrt(total_err_sq[ii][jj]));
00205             }
00206         }
00207 
00208     // Format the output to maximize significant digits:
00209 
00210     int i_total;
00211     if (total_max < 1)
00212         i_total = -1 - (int) (-log10(total_max));
00213     else
00214         i_total = (int) (log10(total_max));
00215 
00216     return i_total;
00217 }
00218 
00219 local real total_coll(real total_sigma[][N_FINAL])
00220 {
00221     real sigma_coll = 0;        
00222 
00223     for (int i_int = 0; i_int < N_INTER; i_int++)
00224         sigma_coll += total_sigma[i_int][merger_binary_1]
00225                           + total_sigma[i_int][merger_binary_2]
00226                           + total_sigma[i_int][merger_binary_3]
00227                           + total_sigma[i_int][merger_escape_1]
00228                           + total_sigma[i_int][merger_escape_2]
00229                           + total_sigma[i_int][merger_escape_3]
00230                           + total_sigma[i_int][triple_merger];
00231     return sigma_coll;
00232 }
00233 
00234 #define USAGE "usage: rate3 \
00235 [-a #] [-A #] [-c #] [-C #] [-d] [-D] [-e #] [-g #] \
00236 [-I] [-m[123] #] [-m #] [-mv #] [-M #] [-n #] \
00237 [-q] [-Q] [-r[123] # ] [-s #] [-S] [-v #] \n"
00238 
00239 main(int argc, char **argv)
00240 {
00241     int  seed   = 0;
00242     int  n_iter = 1;
00243 
00244     bool intermediate_sigma = TRUE;
00245 
00246     int  scatter_summary_level = 0;     // No low-level output.
00247     int  sigma_summary_flag = 0;
00248     bool physical_units_flag = 1;       // Specify input in physical units.
00249 
00250     real cpu_time_check = 3600;         // One check per CPU hour.
00251     real dt_snap = VERY_LARGE_NUMBER;
00252     real snap_cube_size = 10;
00253 
00254     char* indent = "12345678901234567890";
00255 
00256     // Default dimensionless parameters:
00257 
00258     real m = 0.5, M = 0.5;
00259     bool mM_flag = FALSE;               // To check for incorrect input.
00260 
00261     // Default physical parameters:
00262 
00263     real m1 = 1, m2 = 1, m3 = 1;
00264     bool m123_flag = FALSE;             // To check for incorrect input.
00265     real sma = 1000;                    // 5 A. U. (approx)
00266     real v_rms, mv = 1;
00267 
00268     real r1 = 0, r2 = 0, r3 = 0;
00269     bool v_flag = FALSE;                // Defaults depend on context.
00270 
00271     // Units:   mass = 1 solar mass
00272     //          radius = 1 solar radius
00273     //          velocity = 1 km/s
00274     //          binary separation = 1 solar radius
00275 
00276     // Find which version we are running:
00277 
00278     bool pvm = false;
00279     if (strstr(argv[0], ".pvm")) {
00280 
00281 #ifndef HAS_PVM                                 // Compile-time check.
00282         err_exit("PVM not available");
00283 #endif
00284         if (getenv("PVM_ROOT") == NULL)
00285             err_exit("PVM not available");      // Run time check...
00286 
00287         pvm = true;
00288     }
00289 
00290     check_help();
00291 
00292     scatter_profile prof;
00293     make_standard_profile(prof);
00294 
00295     // Parse the command line (the hard way):
00296 
00297     int i = 0;
00298     while (++i < argc) if (argv[i][0] == '-')
00299         switch (argv[i][1]) {
00300             case 'a': sma = atof(argv[++i]);
00301                       break;
00302             case 'A': prof.eta = atof(argv[++i]);
00303                       break;
00304             case 'c': cpu_time_check = 3600*atof(argv[++i]);  // (In hours)
00305                       break;
00306             case 'C': if (!pvm) 
00307                           snap_cube_size = atof(argv[++i]);
00308                       else
00309                           cerr << "\"-C\" option disallowed in PVM mode\n";
00310                       break;
00311             case 'd': physical_units_flag = 1 - physical_units_flag;
00312                       break;
00313             case 'D': if (!pvm) {
00314                           dt_snap = atof(argv[++i]);
00315                           scatter_summary_level = 2; // Undo with later "-q/Q/S"
00316                           sigma_summary_flag = 1;
00317                       } else
00318                           cerr << "\"-D\" option disallowed in PVM mode\n";
00319                       break;
00320             case 'e': prof.ecc = atof(argv[++i]);
00321                       prof.ecc_flag = 1;
00322                       break;
00323             case 'g': prof.tidal_tol_factor = atof(argv[++i]);
00324                       break;
00325             case 'I': intermediate_sigma = 1 - intermediate_sigma;
00326                       break;
00327             case 'm': switch(argv[i][2]) {
00328                           case '\0':    m = atof(argv[++i]);
00329                                         mM_flag = TRUE;
00330                                         break;
00331                           case '1':     m1 = atof(argv[++i]);
00332                                         m123_flag = TRUE;
00333                                         break;
00334                           case '2':     m2 = atof(argv[++i]);
00335                                         m123_flag = TRUE;
00336                                         break;
00337                           case '3':     m3 = atof(argv[++i]);
00338                                         m123_flag = TRUE;
00339                                         break;
00340                           case 'v':     mv = atof(argv[++i]);
00341                                         break;
00342                           default:      cerr << USAGE;
00343                                         get_help();
00344                       }
00345                       break;
00346             case 'M': M = atof(argv[++i]);
00347                       mM_flag = TRUE;
00348                       break;
00349             case 'n': n_iter = atoi(argv[++i]);
00350                       break;
00351             case 'q': if (scatter_summary_level > 0)
00352                           scatter_summary_level = 0;
00353                       else
00354                           scatter_summary_level = 1;
00355                       break;
00356             case 'Q': if (scatter_summary_level > 0)
00357                           scatter_summary_level = 0;
00358                       else
00359                           scatter_summary_level = 2;
00360                       break;
00361             case 'r': switch(argv[i][2]) {
00362                           case '1':     r1 = atof(argv[++i]);
00363                                         break;
00364                           case '2':     r2 = atof(argv[++i]);
00365                                         break;
00366                           case '3':     r3 = atof(argv[++i]);
00367                                         break;
00368                       }
00369                       if (r3 < 0) err_exit("r3 < 0"); // To avoid bus error in
00370                                                       // optimized g++ 2.3.3!
00371                       break;
00372             case 's': seed = atoi(argv[++i]);
00373                       break;
00374             case 'S': sigma_summary_flag = 1 - sigma_summary_flag;
00375                       break;
00376             case 'v': v_rms = atof(argv[++i]);
00377                       v_flag = TRUE;
00378                       break;
00379             default:  cerr << USAGE;
00380                       get_help();
00381         }
00382 
00383     int first_seed = srandinter(seed);
00384     cerr << "Thermally averaged reaction rates, random seed = "
00385          << first_seed << endl;
00386 
00387     real m_unit, r_unit, v_unit, v_rel_th;
00388 
00389     if (physical_units_flag) {
00390 
00391         // Convert to internal units:
00392 
00393         if (mM_flag) cerr << "Ignoring m/M settings." << endl;
00394         if (!v_flag) v_rms = 10; // Default
00395 
00396         m_unit = m1 + m2;
00397         r_unit = sma;
00398         v_unit = sqrt( (Grav * Msun / Rsun)
00399                            * m1 * m2 * (m1 + m2 + m3) / ((m1 + m2) * m3)
00400                            / sma ) / Kms;
00401         prof.m2 = m2 / m_unit;
00402         prof.m3 = m3 / m_unit;
00403         v_rel_th = sqrt( mv/(m1 + m2) + mv/m3 ) * v_rms / v_unit;
00404 
00405         indent = "                 ";
00406         cerr << "physical units:  m1 = " << m1 << "  m2 = " << m2
00407              << "  m3 = " << m3 << "  mv = " << mv << "  (Msolar)\n";
00408         cerr << indent << "binary semi-major axis = " << sma << "  (Rsolar)\n";
00409         cerr << indent << "v_rms = " << v_rms
00410                        << "  v_crit = " << v_unit << "  (km/s, 3D)\n";
00411 
00412     } else {
00413 
00414         if (m123_flag) cerr << "Ignoring m[123] settings." << endl;
00415         if (!v_flag) v_rms = 1; // Default
00416 
00417         m_unit = 1;
00418         r_unit = 1;
00419         v_unit = 1;
00420         prof.m2 = m;
00421         prof.m3 = M;
00422         v_rel_th = v_rms;
00423 
00424         indent = "                      ";
00425         cerr << "dimensionless units:  m = " << m << "  M = " << M << endl;
00426     }
00427 
00428     prof.r1 = r1 / r_unit;
00429     prof.r2 = r2 / r_unit;
00430     prof.r3 = r3 / r_unit;
00431 
00432     cerr << indent << "r1 = " << r1 << "  r2 = " << r2
00433                    << "  r3 = " << r3;
00434     if (physical_units_flag) cerr << "  (Rsolar)";
00435     cerr << endl;
00436 
00437     cerr << indent << "binary eccentricity = ";
00438     if (prof.ecc_flag)
00439         cerr << prof.ecc << endl;
00440     else
00441         cerr << "thermal [f(e) = 2e]\n";
00442 
00443     cerr << indent << "v_rel_th / v_crit = " << v_rel_th << endl;
00444 
00445     prof.v_inf = -1;
00446     print_profile(cerr, prof);
00447 
00448     cpu_init();
00449     real cpu_init = cpu_time();
00450     real cpu_save = cpu_init;
00451 
00452     real max_dens, init_dens = 4;
00453     int nv, init_nv = 2;
00454 
00455     if (!intermediate_sigma) {
00456         for (; n_iter > 1; n_iter--, init_dens *= 4, init_nv *= 2);
00457         n_iter = 1;
00458     }
00459 
00460     // Loop over pairs of (max_dens, nv):
00461 
00462     for (max_dens = init_dens, nv = init_nv; n_iter > 0;
00463          n_iter--, max_dens *= 4, nv *= 2) {    // Best joint choice of params?
00464 
00465         cerr.precision(STD_PRECISION);
00466 
00467         real total_sigma[N_INTER][N_FINAL];
00468         real total_err_sq[N_INTER][N_FINAL];
00469 
00470         cerr << endl << "max_dens = " << max_dens
00471              << ",  nv = " << nv << endl;
00472 
00473         int total_trials = get_rate3(prof, v_rel_th, nv, max_dens,
00474                                      total_sigma, total_err_sq,
00475                                      cpu_time_check, cpu_init, cpu_save,
00476                                      dt_snap, snap_cube_size,
00477                                      scatter_summary_level,
00478                                      sigma_summary_flag);
00479 
00480         cerr << "\ntotal_trials = " << total_trials << "\n\n";
00481 
00482         cerr << "<v sigma> (unit = pi a^2 v_rel_th):\n\n";
00483         print_sigma3_array(total_sigma, 1);
00484 
00485         cerr << "<v sigma_err> (unit = pi a^2 v_rel_th):\n\n";
00486         print_sigma3_err_array(total_err_sq, 1);
00487 
00488         // More detail:
00489 
00490         real sigma_coll = total_coll(total_sigma);
00491 
00492         cerr << "<v sigma> and <v sigma_err> (details):\n\n";
00493         print_sigma3_nonmergers(total_sigma, total_err_sq, 1);
00494 
00495         if (sigma_coll > 0) 
00496             print_sigma3_mergers(total_sigma, total_err_sq, 1);
00497 
00498         if (physical_units_flag) {
00499 
00500             int i_total = rescale_physical(sma, v_rel_th, v_unit,
00501                                            total_sigma, total_err_sq);
00502             real print_factor = pow(10, -i_total); 
00503         
00504             cerr << "<v sigma> (unit = 10^"
00505                  << i_total << " pc^3/Gyr):\n\n";
00506             print_sigma3_array(total_sigma, print_factor);
00507 
00508             // Use the *same* print_factor for the errors as for sigma:
00509 
00510             cerr << "<v sigma_err> (unit = 10^"
00511                  << i_total << " pc^3/Gyr):\n\n";
00512             print_sigma3_err_array(total_err_sq, print_factor);
00513 
00514             Details:
00515 
00516             cerr << "<v sigma> and <v sigma_err> (details):\n\n";
00517             print_sigma3_nonmergers(total_sigma, total_err_sq, print_factor);
00518 
00519             if (sigma_coll > 0) 
00520                 print_sigma3_mergers(total_sigma, total_err_sq, print_factor);
00521         }
00522     }
00523 }
00524 
00525 #endif