delta_a.C Source File

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