out3.C Source File

00001 
00002 // out3.C: Perform a series of three-body scattering experiments
00003 //         and print out the result of each.
00004 
00005 // Starlab application:  scatter3.
00006 
00007 #include "scatter3.h"
00008 
00009 local void print(sdyn3* b)
00010     {
00011     fprintf(stdout, "%.5f ", b->get_time() + b->get_time_offset());
00012     for_all_daughters(sdyn3, b, bb)
00013         for (int k = 0; k < 3; k++) fprintf(stdout, "%.5f ", bb->get_pos()[k]);
00014     fprintf(stdout, "\n");
00015     }
00016 
00017 main(int argc, char **argv)
00018     {
00019 
00020     initial_state3 init;
00021     make_standard_init(init);
00022 
00023     int  seed       = 0;        // seed for random number generator
00024     int n_rand      = 0;        // number of times to invoke the generator
00025                                 // before starting for real
00026     int  n_experiments = 1;     // default: only one run
00027     real dt_out     =           // output time interval
00028           VERY_LARGE_NUMBER;
00029     real dt_snap    =           // output time interval
00030           VERY_LARGE_NUMBER;
00031     real dt_print    =          // print time interval
00032           VERY_LARGE_NUMBER;
00033 
00034     real cpu_time_check = 3600;
00035     real snap_cube_size = 10;
00036 
00037     int planar_flag = 0;
00038 
00039     bool  b_flag = FALSE;
00040     bool  q_flag = FALSE;
00041     bool  Q_flag = FALSE;
00042 
00043     extern char *poptarg;
00044     int c;
00045     char* param_string = "A:bc:C:d:D:e:L:m:M:n:N:pPqQr:R:s:S:U:v:x:y:z:";
00046 
00047     while ((c = pgetopt(argc, argv, param_string)) != -1)
00048         switch(c)
00049             {
00050             case 'A': init.eta = atof(poptarg);
00051                       break;
00052             case 'b': b_flag = 1 - b_flag;
00053                       break;
00054             case 'c': cpu_time_check = 3600*atof(poptarg);// (Specify in hours)
00055                       break;
00056             case 'C': snap_cube_size = atof(poptarg);
00057                       break;
00058             case 'd': dt_out = atof(poptarg);
00059                       break;
00060             case 'D': dt_print = atof(poptarg);
00061                       break;
00062             case 'e': init.ecc = atof(poptarg);
00063                       break;
00064             case 'L': init.r_init_min = atof(poptarg);
00065                       break;
00066             case 'm': init.m2 = atof(poptarg);
00067                       break;
00068             case 'M': init.m3 = atof(poptarg);
00069                       break;
00070             case 'n': n_experiments = atoi(poptarg);
00071                       break;
00072             case 'N': n_rand = atoi(poptarg);
00073                       break;
00074             case 'p': planar_flag = 1;
00075                       break;
00076             case 'P': planar_flag = -1;
00077                       break;
00078             case 'q': q_flag = 1 - q_flag;
00079                       break;
00080             case 'Q': Q_flag = 1 - Q_flag;
00081                       break;
00082             case 'r': init.rho = atof(poptarg);
00083                       break;
00084             case 'R': init.r_stop = init.r_init_min
00085                                   = init.r_init_max
00086                                   = atof(poptarg);
00087                       break;
00088             case 's': seed = atoi(poptarg);
00089                       break;
00090             case 'S': init.r_stop = atof(poptarg);
00091                       break;
00092             case 'U': init.r_init_max = atof(poptarg);
00093                       break;
00094             case 'v': init.v_inf = atof(poptarg);
00095                       break;
00096             case 'x': init.r1 = atof(poptarg);
00097                       break;
00098             case 'y': init.r2 = atof(poptarg);
00099                       break;
00100             case 'z': init.r3 = atof(poptarg);
00101                       break;
00102             case '?': params_to_usage(cerr, argv[0], param_string);
00103                       exit(1);
00104             }            
00105 
00106     if (Q_flag) q_flag = TRUE;
00107 
00108     if (init.m2 > 1)
00109         {
00110         cerr << "sigma3: init.m2 = " << init.m2 << " > 1" << endl;
00111         exit(1);
00112         }
00113 
00114     cpu_init();
00115     int random_seed = srandinter(seed, n_rand);
00116 
00117     while (n_experiments--)
00118         {
00119         cerr << "Random seed = " << get_initial_seed()
00120              << "  n_rand = " << get_n_rand() << flush;
00121 
00122         randomize_angles(init.phase);
00123 
00124         if (planar_flag == 1)
00125             init.phase.cos_theta = 1;   // Planar prograde
00126         else if (planar_flag == -1)
00127             init.phase.cos_theta = -1;  // Planar retrograde
00128 
00129         intermediate_state3 inter;
00130         final_state3 final;
00131 
00132         real cpu = cpu_time();  
00133         scatter3(init, inter, final, cpu_time_check,
00134                  dt_out, dt_snap, snap_cube_size,
00135                  dt_print, print);
00136         cpu = cpu_time() - cpu;
00137 
00138         cerr << ":  ";
00139         print_scatter3_outcome(inter, final, cerr);
00140 
00141         if (Q_flag) print_scatter3_summary(inter, final, cpu, cerr);
00142 
00143         if (!q_flag) print_scatter3_report(init, inter, final, b_flag, cerr);
00144 
00145         }
00146     }