template.C Source File

00001 
00002 // template.C: Generic program to determine differential cross-sections
00003 //             or perform other detailed calculations in the "sigma3"
00004 //             environment.  Additional statistics are gathered and
00005 //             reduced locally, using hooks to the library routines.
00006 
00007 // Starlab application:  get_sigma3.
00008 
00009 #include "sigma3.h"
00010 
00011 // Global statistics data:
00012 
00013 #define NA 10   // Number of bins in semi-major-axis a
00014 #define NE 10   // Number of bins in eccentricity e
00015 
00016 static int n = 0;
00017 static int counter[NA][NE][N_RHO_ZONE_MAX];     // N_RHO_ZONE_MAX is defined
00018                                                 // in header file sigma3.h
00019 
00020 static real sigma[NA][NE];                      // Cross-sections and errors
00021 static real sigma_err_sq[NA][NE];               // for specific process(es)
00022 
00023 local void initialize_stats()
00024 {
00025     for (int ia = 0; ia < NA; ia++)
00026         for (int je = 0; je < NE; je++)
00027             for (int kr = 0; kr < N_RHO_ZONE_MAX; kr++)
00028                 counter[ia][je][kr] = 0;
00029 }
00030 
00031 // accumulate_stats: Gather statistics supplied by get_sigma3.
00032 //                   Control will be transferred here after each scattering;
00033 //                   all relevent state information and the rho zone used
00034 //                   are provided by get_sigma3.
00035 //                   Customize here to the specific application at hand.
00036 
00037 local void accumulate_stats(scatter_profile& prof,
00038                             initial_state3& init,
00039                             intermediate_state3& inter,
00040                             final_state3& final,
00041                             int rho_zone,
00042                             sigma_out& out)
00043 {
00044     if (final.descriptor == exchange_1 || final.descriptor == exchange_2) {
00045 
00046 //      cerr << "accumulate_stats: " << state_string(inter.descriptor)
00047 //           << " " << state_string(final.descriptor) << endl;
00048 
00049         if (prof.r3 > 0) {
00050 
00051             real ratio = final.sma / prof.r3;
00052 
00053             int ia = (int) (log(ratio) / log(2.0));
00054             if (ia >= NA) ia = NA - 1;
00055 
00056             int je = (int) (10 * final.ecc);
00057             if (je >= NE) je = NE - 1;
00058 
00059             n++;
00060             counter[ia][je][rho_zone]++;
00061 
00062         }
00063     }
00064 }
00065 
00066 // normalize_counts:  Convert raw counts into cross-sections and errors.
00067 //                    Weights are provided by a library routine to minimise
00068 //                    duplication in functionality and to allow the internal
00069 //                    workings of the package to change without affecting
00070 //                    user-written code.
00071 
00072 local void normalize_counts(sigma_out& out)
00073 {
00074     for (int ia = 0; ia < NA; ia++)
00075         for (int je = 0; je < NE; je++) {
00076 
00077             sigma[ia][je] = 0;
00078             sigma_err_sq[ia][je] = 0;
00079 
00080             for (int kr = 0; kr <= out.i_max; kr++) {
00081 
00082                 real weight = zone_weight(out, kr);     // = zone_area
00083                                                         //    / trials_per_zone
00084 
00085                 sigma[ia][je] += weight * counter[ia][je][kr];
00086                 sigma_err_sq[ia][je] += weight * weight * counter[ia][je][kr];
00087             }
00088         }
00089 } 
00090 
00091 local void print_stats(scatter_profile& prof, sigma_out& out)
00092 {
00093     normalize_counts(out);
00094     real v2 = prof.v_inf * prof.v_inf;
00095 
00096     int ia, je;
00097 
00098     cerr << "\nDifferential cross sections:\n";
00099 
00100     cerr << "\nRaw counts (column = log2(sma/r3), row = 10*ecc, total = "
00101          << n << "):\n\n";
00102     for (ia = 0; ia < NA; ia++) {
00103         for (je = 0; je < NE; je++) {
00104             int total = 0;
00105             for (int kr = 0; kr <= out.i_max; kr++)
00106                 total += counter[ia][je][kr];
00107             fprintf(stderr, " %7d", total);
00108         }
00109         cerr << endl;
00110     }
00111 
00112     cerr << "\nv^2 sigma / (pi a^2)\n\n";
00113     for (ia = 0; ia < NA; ia++) {
00114         for (je = 0; je < NE; je++)
00115             fprintf(stderr, " %7.3f", v2 * sigma[ia][je]);
00116         cerr << endl;
00117     }
00118 
00119     cerr << "\nv^2 (sigma error) / (pi a^2)\n\n";
00120     for (ia = 0; ia < NA; ia++) {
00121         for (je = 0; je < NE; je++)
00122             fprintf(stderr, " %7.3f", v2 * sqrt(sigma_err_sq[ia][je]));
00123         cerr << endl;
00124     }
00125 }
00126 
00127 //----------------------------------------------------------------------
00128 //----------------------------------------------------------------------
00129 
00130 // No particular need to modify anything below this line.
00131 
00132 //------------------------------------------------------------------------
00133 //------------------------------------------------------------------------
00134 
00135 main(int argc, char **argv)
00136     {
00137     int  debug  = 0;
00138     int  seed   = 0;
00139     int  n_rand = 0;
00140     real max_trial_density = 1.0;
00141 
00142     real cpu_time_check = 3600; // One check per CPU hour!
00143     real dt_snap = VERY_LARGE_NUMBER;
00144     real snap_cube_size = 10;
00145     int scatter_summary_level = 0;
00146 
00147     bool print_counts = FALSE;
00148 
00149     scatter_profile prof;
00150     make_standard_profile(prof);
00151 
00152     extern char *poptarg;
00153     int c;
00154     char* param_string = "A:c:C:D:d:e:m:M:N:pqQs:v:x:y:z:";
00155 
00156     while ((c = pgetopt(argc, argv, param_string)) != -1)
00157         switch(c)
00158             {
00159             case 'A': prof.eta = atof(poptarg);
00160                       break;
00161             case 'c': cpu_time_check = 3600*atof(poptarg);// (Specify in hours)
00162                       break;
00163             case 'C': snap_cube_size = atof(poptarg);
00164                       break;
00165             case 'd': max_trial_density = atof(poptarg);
00166                       break;
00167             case 'D': dt_snap = atof(poptarg);
00168                       scatter_summary_level = 2;  // Undo with later "-q/Q"
00169                       break;
00170             case 'e': prof.ecc = atof(poptarg);
00171                       prof.ecc_flag = 1;
00172                       break;
00173             case 'm': prof.m2 = atof(poptarg);
00174                       break;
00175             case 'M': prof.m3 = atof(poptarg);
00176                       break;
00177             case 'N': n_rand = atoi(poptarg);
00178                       break;
00179             case 'p': print_counts = 1 - print_counts;
00180                       break;
00181             case 'q': if (scatter_summary_level > 0)
00182                           scatter_summary_level = 0;
00183                       else
00184                           scatter_summary_level = 1;
00185                       break;
00186             case 'Q': if (scatter_summary_level > 0)
00187                           scatter_summary_level = 0;
00188                       else
00189                           scatter_summary_level = 2;
00190                       break;
00191             case 's': seed = atoi(poptarg);
00192                       break;
00193             case 'v': prof.v_inf = atof(poptarg);
00194                       break;
00195             case 'V': debug = atoi(poptarg);
00196                       break;
00197             case 'x': prof.r1 = atof(poptarg);
00198                       break;
00199             case 'y': prof.r2 = atof(poptarg);
00200                       break;
00201             case 'z': prof.r3 = atof(poptarg);
00202                       break;
00203             case '?': params_to_usage(cerr, argv[0], param_string);
00204                       exit(1);
00205             }
00206 
00207     // Debugging:    |debug| = 0 ==> no debugging output
00208     //               |debug| = 1 ==> output only at end (default)
00209     //               |debug| = 2 ==> output at end of each top-level iteration
00210     //               |debug| = 3 ==> output after each sub-trial
00211     //
00212     //                debug  < 0 ==> simple statistics also (default: off)
00213 
00214     cpu_init();
00215 
00216     sigma_out out;
00217     int first_seed = srandinter(seed, n_rand);
00218 
00219     cerr << "Random seed = " << first_seed << endl;
00220     print_profile(cerr, prof);
00221 
00222     initialize_stats();
00223 
00224     // Generic scattering experiment call, but with additional hook to
00225     // transfer control to accumulate_stats (above) after each scattering.
00226 
00227     get_sigma3(max_trial_density, prof, out,
00228                debug, cpu_time_check,
00229                dt_snap, snap_cube_size,
00230                scatter_summary_level,
00231                accumulate_stats);
00232 
00233     // Out is the standard final descriptor of the scattering experiment.
00234     // It contains useful statistical and diagnostic information about
00235     // the details of what went on, but it may not be needed in all cases.
00236 
00237     // Standard "sigma3" output:
00238 
00239     print_sigma3(out, prof.v_inf * prof.v_inf);
00240     if (print_counts) print_all_sigma3_counts(out, cerr);
00241 
00242     // User-specific output:
00243 
00244     print_stats(prof, out);
00245 }