delta_e.C Source File

00001 
00002 // delta_e.C:   Generate statistics on changes in binary energy
00003 //              following encounters.
00004 
00005 // Starlab application:  get_sigma3.
00006 
00007 #include "sigma3.h"
00008 
00009 #define DEBUG 0
00010 
00011 #define MIN_OSC 5           // Defines a "strong" resonance
00012 
00013 // Global statistics data:
00014 
00015 #define NBIN 5              // defined bins:    0: direct exchange
00016                             //                  1: hier_res
00017                             //                  2: strong dem_res
00018                             //                  3: weak dem_res pres
00019                             //                  4: weak dem_res exch
00020 
00021 #define NR N_RHO_ZONE_MAX   // for brevity (N_RHO_ZONE_MAX defined in sigma3.h)
00022 
00023 // Subdivide hits by bin and radial zone:
00024 
00025 static int  n_total[NBIN][NR];
00026 static real de_total[NBIN][NR];
00027 
00028 // Keep track of individual scatterings in each bin & zone:
00029 
00030 #define NMAX 25000          // should be adequate...
00031 
00032 static int  nscat[NBIN][NR];
00033 static real descat[NBIN][NR][NMAX];
00034 
00035 local void initialize_stats()
00036 {
00037     for (int ibin = 0; ibin < NBIN; ibin++)
00038         for (int kr = 0; kr < NR; kr++) {
00039             n_total[ibin][kr] = 0;
00040             de_total[ibin][kr] = 0;
00041             nscat[ibin][kr] = 0;
00042         }
00043 }
00044 
00045 // accumulate_stats: Gather statistics supplied by get_sigma3.
00046 //                   Control will be transferred here after each scattering;
00047 //                   all relevent state information and the rho zone used
00048 //                   are provided by get_sigma3.
00049 //                   Customize here to the specific application at hand.
00050 
00051 local void accumulate_stats(scatter_profile& prof,
00052                             initial_state3& init,
00053                             intermediate_state3& inter,
00054                             final_state3& final,
00055                             int rho_zone,
00056                             sigma_out& out)
00057 {
00058     if (inter.n_stars == 3 &&
00059         (inter.descriptor == hierarchical_resonance
00060          || inter.descriptor == democratic_resonance
00061          || final.descriptor == exchange_1
00062          || final.descriptor == exchange_2)) {
00063 
00064         real m1 = 1 - init.m2;
00065         real m2 = init.m2;
00066 
00067         real e_init = 0.5*m1*m2 / init.sma;
00068 
00069         if (final.escaper == 1)
00070             m1 = init.m3;
00071         else if (final.escaper == 2)
00072             m2 = init.m3;
00073         real e_final = 0.5*m1*m2 / final.sma;
00074 
00075         real de = e_final/e_init - 1;
00076 
00077         if (DEBUG) {
00078             int p = cerr.precision(STD_PRECISION);
00079             cerr << "accumulate_stats: " << state_string(inter.descriptor)
00080                  << " " << state_string(final.descriptor)
00081                  << ",  n_osc = " << inter.n_osc << endl;
00082             cerr << "  rho_zone = " << rho_zone;
00083             cerr << "  escaper = " << final.escaper;
00084             cerr << "  mass = " << final.system[final.escaper-1].mass;
00085             cerr << "  e_init = " << e_init;
00086             cerr << "  de/e = " << de << endl;
00087             cerr.precision(p);
00088         }
00089 
00090         int ibin = 0;
00091         if (inter.descriptor == hierarchical_resonance)
00092             ibin = 1;
00093         else if (inter.descriptor == democratic_resonance) {
00094             if (inter.n_osc >= MIN_OSC)
00095                 ibin = 2;
00096             else {
00097                 if (final.descriptor == preservation)
00098                     ibin = 3;
00099                 else
00100                     ibin = 4;
00101             }
00102         }
00103 
00104         // Store statistical information.
00105 
00106         n_total[ibin][rho_zone]++;
00107         de_total[ibin][rho_zone] += de;
00108         descat[ibin][rho_zone][nscat[ibin][rho_zone]++] = de;
00109     }
00110 }
00111 
00112 local int realcomp(const void * a, const void * b)      // For use by qsort
00113 {
00114     if (*((real *) a) > *((real *) b))
00115         return 1;
00116     else if (*((real *) a) < *((real *) b))
00117         return -1;
00118     else
00119         return 0;
00120 }
00121 
00122 // print_stats:  Convert raw counts into final data and print them out.
00123 //               Weights are provided by a library routine to minimise
00124 //               duplication in functionality and to allow the internal
00125 //               workings of the package to change without affecting
00126 //               user-written code.
00127 
00128 local void print_stats(scatter_profile& prof, sigma_out& out)
00129 {
00130     char* s[NBIN] = {"non-resonant exchanges       ",
00131                      "hierarchical resonances      ",
00132                      "strong democratic resonances ",
00133                      "weak democratic preservations",
00134                      "weak democratic exchanges    "};
00135 
00136     for (int ibin = 0; ibin < NBIN; ibin++) {
00137 
00138         real sigma = 0;
00139         real de = 0;
00140         int nt = 0;
00141 
00142         cerr << endl << s[ibin] << endl;
00143 
00144         real w_min = VERY_LARGE_NUMBER, w_max = 0;
00145 
00146         for (int kr = 0; kr < out.i_max; kr++) {
00147 
00148             real w = zone_weight(out, kr);  // = zone_area / trials_per_zone
00149 
00150             w_min = min(w_min, w);
00151             w_max = max(w_max, w);
00152 
00153             // May as well determine the cross-section too, as a check.
00154 
00155             sigma += w * n_total[ibin][kr];
00156             de += w * de_total[ibin][kr];       // Note: de_total includes
00157                                                 //       a factor of n_total
00158             nt +=nscat[ibin][kr] ;
00159 
00160             cerr << "    radial zone " << kr << "  (weight = " << w
00161                  << "):  n_total = " << n_total[ibin][kr]
00162                  << ",  <dE/E> = " << de_total[ibin][kr]
00163                                         / max(1, n_total[ibin][kr])
00164                  << endl;
00165 
00166             if (nscat[ibin][kr] > 3) {
00167 
00168                 cerr << "        dE/E quartiles:";
00169 
00170                 qsort((void*)&descat[ibin][kr][0], (size_t)nscat[ibin][kr],
00171                       sizeof(real), realcomp);
00172                 real factor = 0.25*nscat[ibin][kr];
00173                 for (int i = 1; i <= 3; i++)
00174                     cerr << "  " << descat[ibin][kr][(int)(factor*i)];
00175 
00176                 cerr << endl;
00177 
00178                 if (DEBUG)
00179                     for (int i = 0; i < nscat[ibin][kr]; i++)
00180                         cerr << "            " << descat[ibin][kr][i] << endl;
00181 
00182             }
00183         }
00184 
00185         if (sigma > 0) de /= sigma;
00186 
00187         cerr << "<dE/E> for " << s[ibin] << " = " << de
00188              << ",  v^2 sigma = " << prof.v_inf * prof.v_inf * sigma << endl;
00189 
00190         if (nt > 3) {
00191 
00192             // Combine all radial zones.
00193 
00194             cerr << "dE/E quartiles:";
00195 
00196             real* tempde = new real[out.i_max*nt*(int)(w_max/w_min + 0.1)];
00197             nt = 0;
00198 
00199             // Note quick and dirty treatment of zone weights...
00200 
00201             for (int kr = 0; kr <= out.i_max; kr++)
00202                 for (int i = 0; i < nscat[ibin][kr]; i++)
00203                     for (int ww = 0;
00204                          ww < (int)(zone_weight(out, kr)/w_min + 0.1); ww++)
00205                         tempde[nt++] = descat[ibin][kr][i];
00206 
00207             qsort((void*)tempde, (size_t)nt, sizeof(real), realcomp);
00208             real factor = 0.25*nt;
00209             for (int i = 1; i <= 3; i++)
00210                 cerr << "  " << tempde[(int)(factor*i)];
00211 
00212             cerr << endl;
00213 
00214             if (DEBUG)
00215                 for (int i = 0; i < nt; i++)
00216                     cerr << "        " << tempde[i] << endl;
00217 
00218         }
00219     }
00220 }
00221 
00222 //----------------------------------------------------------------------
00223 //----------------------------------------------------------------------
00224 
00225 // Probably no particular need to modify anything below this line.
00226 
00227 //----------------------------------------------------------------------
00228 //----------------------------------------------------------------------
00229 
00230 main(int argc, char **argv)
00231 {
00232     int  debug  = 0;
00233     int  seed   = 0;
00234     int  n_rand = 0;
00235     real max_trial_density = 1.0;
00236 
00237     real cpu_time_check = 3600; // One check per CPU hour!
00238     real dt_snap = VERY_LARGE_NUMBER;
00239     real snap_cube_size = 10;
00240     int scatter_summary_level = 0;
00241 
00242     bool print_counts = FALSE;
00243 
00244     // Find which version we are running:
00245 
00246     bool pvm = false;
00247     if (strstr(argv[0], ".pvm")) {
00248 
00249 #ifndef HAS_PVM                                 // Compile-time check.
00250         err_exit("PVM not available");
00251 #endif
00252         if (getenv("PVM_ROOT") == NULL)
00253             err_exit("PVM not available");      // Run time check...
00254 
00255         pvm = true;
00256     }
00257 
00258     scatter_profile prof;
00259     make_standard_profile(prof);
00260 
00261     extern char *poptarg;
00262     int c;
00263     char* param_string = "A:c:C:D:d:e:m:M:N:pqQs:V:v:x:y:z:";
00264 
00265     while ((c = pgetopt(argc, argv, param_string)) != -1)
00266         switch(c) {
00267 
00268             case 'A': prof.eta = atof(poptarg);
00269                       break;
00270             case 'c': cpu_time_check = 3600*atof(poptarg);// (Specify in hours)
00271                       break;
00272             case 'C': if (!pvm) 
00273                           snap_cube_size = atof(poptarg);
00274                       else
00275                           cerr << "\"-C\" option disallowed in PVM mode\n";
00276                       break;
00277             case 'd': max_trial_density = atof(poptarg);
00278                       break;
00279             case 'D': if (!pvm) {
00280                           dt_snap = atof(poptarg);
00281                           scatter_summary_level = 2;  // Undo with later "-q/Q"
00282                       } else
00283                           cerr << "\"-D\" option disallowed in PVM mode\n";
00284                       break;
00285             case 'e': prof.ecc = atof(poptarg);
00286                       prof.ecc_flag = 1;
00287                       break;
00288             case 'm': prof.m2 = atof(poptarg);
00289                       break;
00290             case 'M': prof.m3 = atof(poptarg);
00291                       break;
00292             case 'N': n_rand = atoi(poptarg);
00293                       break;
00294             case 'p': print_counts = 1 - print_counts;
00295                       break;
00296             case 'q': if (scatter_summary_level > 0)
00297                           scatter_summary_level = 0;
00298                       else
00299                           scatter_summary_level = 1;
00300                       break;
00301             case 'Q': if (scatter_summary_level > 0)
00302                           scatter_summary_level = 0;
00303                       else
00304                           scatter_summary_level = 2;
00305                       break;
00306             case 's': seed = atoi(poptarg);
00307                       break;
00308             case 'v': prof.v_inf = atof(poptarg);
00309                       break;
00310             case 'V': debug = atoi(poptarg);
00311                       break;
00312             case 'x': prof.r1 = atof(poptarg);
00313                       break;
00314             case 'y': prof.r2 = atof(poptarg);
00315                       break;
00316             case 'z': prof.r3 = atof(poptarg);
00317                       break;
00318             case '?': params_to_usage(cerr, argv[0], param_string);
00319                       exit(1);
00320         }
00321 
00322     // Debugging:    |debug| = 0 ==> no debugging output (default)
00323     //               |debug| = 1 ==> output only at end
00324     //               |debug| = 2 ==> output at end of each top-level iteration
00325     //               |debug| = 3 ==> output after each sub-trial
00326     //
00327     //                debug  < 0 ==> simple statistics also (default: off)
00328 
00329     cpu_init();
00330 
00331     sigma_out out;
00332     int first_seed = srandinter(seed, n_rand);
00333 
00334     cerr << "random seed = " << first_seed << endl;
00335     print_profile(cerr, prof);
00336 
00337     initialize_stats();
00338 
00339     // Generic scattering experiment call, but an with additional hook to
00340     // transfer control to accumulate_stats (above) after each scattering.
00341 
00342     get_sigma3(max_trial_density, prof, out,
00343                debug, cpu_time_check,
00344                dt_snap, snap_cube_size,
00345                scatter_summary_level,
00346                accumulate_stats);
00347 
00348     // Out is the standard final descriptor of the scattering experiment.
00349     // It contains useful statistical and diagnostic information about
00350     // the details of what went on, but it may not be needed in all cases.
00351 
00352     // Standard "sigma3" output:
00353 
00354     print_sigma3(out, prof.v_inf * prof.v_inf);
00355     if (print_counts) print_all_sigma3_counts(out, cerr);
00356 
00357     // User-specific output:
00358 
00359     for (int i = 0; i < 72; i++) cerr << "-";
00360     cerr << endl;
00361 
00362     print_stats(prof, out);
00363 }