scatt_stats_1.C Source File

00001 
00002 // scatt_stats_1.C: Determine specific cross-sections for 3-body scattering.
00003 //                  Steve McMillan, Piet Hut, Fred Rasio (April 1994)
00004 
00005 // Starlab application:  get_sigma3.
00006 
00007 #include "sigma3.h"
00008 
00009 // Global statistics data:
00010 
00011 #define NA 10
00012 #define NE 10
00013 
00014 static int n = 0;
00015 static int counter[NA][NE][N_RHO_ZONE_MAX];
00016 
00017 static real sigma[NA][NE];
00018 static real sigma_err_sq[NA][NE];
00019 
00020 local void initialize_stats()
00021 {
00022     for (int ia = 0; ia < NA; ia++)
00023         for (int je = 0; je < NE; je++)
00024             for (int kr = 0; kr < N_RHO_ZONE_MAX; kr++)
00025                 counter[ia][je][kr] = 0;
00026 }
00027 
00028 // accumulate_stats: Gather statistics supplied by get_sigma3. Customize here
00029 //                   to the specific application at hand.
00030 
00031 local void accumulate_stats(scatter_profile& prof,
00032                             initial_state3& init,
00033                             intermediate_state3& inter,
00034                             final_state3& final,
00035                             int rho_bin,
00036                             sigma_out& out)
00037 {
00038     if (final.descriptor == exchange_1 || final.descriptor == exchange_2) {
00039 
00040 //      cerr << "acc_stats: " << state_string(inter.descriptor)
00041 //           << " " << state_string(final.descriptor) << endl;
00042 
00043         if (prof.r3 > 0) {
00044 
00045             real ratio = final.sma / prof.r3;
00046 
00047             int ia = (int) (log(ratio) / log(2.0));
00048             if (ia >= NA) ia = NA - 1;
00049 
00050             int je = (int) (10 * final.ecc);
00051             if (je >= NE) je = NE - 1;
00052 
00053             n++;
00054             counter[ia][je][rho_bin]++;
00055 
00056 //          cerr << "           a/R = " << ratio
00057 //               << "  e = " << final.ecc << endl;
00058         }
00059     }
00060 }
00061 
00062 // normalize_counts: convert raw counts into cross-sections and errors
00063 
00064 local void normalize_counts(sigma_out& out)
00065 {
00066     for (int ia = 0; ia < NA; ia++)
00067         for (int je = 0; je < NE; je++) {
00068 
00069             sigma[ia][je] = 0;
00070             sigma_err_sq[ia][je] = 0;
00071 
00072             real rho_sq_min, rho_sq_max;
00073             int kr;
00074             for (kr = 0, rho_sq_min = 0, rho_sq_max = out.rho_sq_init;
00075                  kr <= out.i_max;
00076                  kr++, rho_sq_min = rho_sq_max, rho_sq_max *= RHO_SQ_FACTOR) {
00077 
00078                 real factor = (rho_sq_max - rho_sq_min) / out.trials_per_zone;
00079                 sigma[ia][je] += factor * counter[ia][je][kr];
00080                 sigma_err_sq[ia][je] += factor * factor * counter[ia][je][kr];
00081             }
00082         }
00083 } 
00084 
00085 local void print_stats(scatter_profile& prof, sigma_out& out)
00086 {
00087     normalize_counts(out);
00088     real v2 = prof.v_inf * prof.v_inf;
00089 
00090     int ia, je;
00091 
00092     cerr << "\nDifferential cross sections:\n";
00093 
00094     cerr << "\nRaw counts (column = log2(sma/r3), row = 10*ecc, total = "
00095          << n << "):\n\n";
00096     for (ia = 0; ia < NA; ia++) {
00097         for (je = 0; je < NE; je++) {
00098             int total = 0;
00099             for (int kr = 0; kr <= out.i_max; kr++)
00100                 total += counter[ia][je][kr];
00101             fprintf(stderr, " %7d", total);
00102         }
00103         cerr << endl;
00104     }
00105 
00106     cerr << "\nv^2 sigma / (pi a^2)\n\n";
00107     for (ia = 0; ia < NA; ia++) {
00108         for (je = 0; je < NE; je++)
00109             fprintf(stderr, " %7.3f", v2 * sigma[ia][je]);
00110         cerr << endl;
00111     }
00112 
00113     cerr << "\nv^2 (sigma error) / (pi a^2)\n\n";
00114     for (ia = 0; ia < NA; ia++) {
00115         for (je = 0; je < NE; je++)
00116             fprintf(stderr, " %7.3f", v2 * sqrt(sigma_err_sq[ia][je]));
00117         cerr << endl;
00118     }
00119 }
00120 
00121 main(int argc, char **argv)
00122     {
00123     int  debug  = 0;
00124     int  seed   = 0;
00125     int  n_rand = 0;
00126     real max_trial_density = 1.0;
00127 
00128     real cpu_time_check = 3600; // One check per CPU hour!
00129     real dt_snap = VERY_LARGE_NUMBER;
00130     real snap_cube_size = 10;
00131     int scatter_summary_level = 0;
00132 
00133     bool print_counts = FALSE;
00134 
00135     scatter_profile prof;
00136     make_standard_profile(prof);
00137 
00138     extern char *poptarg;
00139     int c;
00140     char* param_string = "A:c:C:d:D:e:m:M:N:pqQs:v:V:x:y:z:";
00141 
00142     while ((c = pgetopt(argc, argv, param_string)) != -1)
00143         switch(c)
00144             {
00145             case 'A': prof.eta = atof(poptarg);
00146                       break;
00147             case 'c': cpu_time_check = 3600*atof(poptarg);// (Specify in hours)
00148                       break;
00149             case 'C': snap_cube_size = atof(poptarg);
00150                       break;
00151             case 'd': max_trial_density = atof(poptarg);
00152                       break;
00153             case 'D': dt_snap = atof(poptarg);
00154                       scatter_summary_level = 2;  // Undo with later "-q/Q"
00155                       break;
00156             case 'e': prof.ecc = atof(poptarg);
00157                       prof.ecc_flag = 1;
00158                       break;
00159             case 'm': prof.m2 = atof(poptarg);
00160                       break;
00161             case 'M': prof.m3 = atof(poptarg);
00162                       break;
00163             case 'N': n_rand = atoi(poptarg);
00164                       break;
00165             case 'p': print_counts = 1 - print_counts;
00166                       break;
00167             case 'q': if (scatter_summary_level > 0)
00168                           scatter_summary_level = 0;
00169                       else
00170                           scatter_summary_level = 1;
00171                       break;
00172             case 'Q': if (scatter_summary_level > 0)
00173                           scatter_summary_level = 0;
00174                       else
00175                           scatter_summary_level = 2;
00176                       break;
00177             case 's': seed = atoi(poptarg);
00178                       break;
00179             case 'v': prof.v_inf = atof(poptarg);
00180                       break;
00181             case 'V': debug = atoi(poptarg);
00182                       break;
00183             case 'x': prof.r1 = atof(poptarg);
00184                       break;
00185             case 'y': prof.r2 = atof(poptarg);
00186                       break;
00187             case 'z': prof.r3 = atof(poptarg);
00188                       break;
00189             case '?': params_to_usage(cerr, argv[0], param_string);
00190                       exit(1);
00191             }
00192 
00193     // Debugging:    |debug| = 0 ==> no debugging output
00194     //               |debug| = 1 ==> output only at end (default)
00195     //               |debug| = 2 ==> output at end of each top-level iteration
00196     //               |debug| = 3 ==> output after each sub-trial
00197     //
00198     //                debug  < 0 ==> simple statistics also (default: off)
00199 
00200     cpu_init();
00201 
00202     sigma_out out;
00203     int first_seed = srandinter(seed, n_rand);
00204 
00205     cerr << "Random seed = " << first_seed << endl;
00206     print_profile(cerr, prof);
00207 
00208     initialize_stats();
00209 
00210     get_sigma3(max_trial_density, prof, out,
00211                debug, cpu_time_check,
00212                dt_snap, snap_cube_size,
00213                scatter_summary_level,
00214                accumulate_stats);
00215 
00216     print_sigma3(out, prof.v_inf * prof.v_inf);
00217     if (print_counts) print_all_sigma3_counts(out, cerr);
00218 
00219     print_stats(prof, out);
00220 }