pscat.C Source File

00001 
00002 #include "scatter3.h"
00003 #include "pvm_scatt.h"
00004 #include "pvm3.h"
00005 
00006 main(int argc, char **argv)
00007 {
00008     initial_state3 init;
00009     make_standard_init(init);
00010 
00011     int  seed       = 0;        // seed for random number generator
00012     int n_rand      = 0;        // number of times to invoke the generator
00013                                 // before starting for real
00014     int  n_experiments = 1;     // default: only one run
00015     real dt_out     =           // output time interval
00016           VERY_LARGE_NUMBER;
00017     real dt_snap    =           // output time interval
00018           VERY_LARGE_NUMBER;
00019 
00020     real cpu_time_check = 3600;
00021     real snap_cube_size = 10;
00022 
00023     int planar_flag = 0;
00024     int psi_flag = 0;
00025     real psi = 0;
00026 
00027     bool  b_flag = FALSE;
00028     bool  q_flag = FALSE;
00029     bool  Q_flag = FALSE;
00030 
00031     extern char *poptarg;
00032     int c;
00033     char* param_string = "A:bc:C:d:D:e:g:L:m:M:n:N:o:pPqQr:R:s:S:U:v:x:y:z:";
00034 
00035     while ((c = pgetopt(argc, argv, param_string)) != -1)
00036         switch(c) {
00037 
00038             case 'A': init.eta = atof(poptarg);
00039                       break;
00040             case 'b': b_flag = 1 - b_flag;
00041                       break;
00042             case 'c': cpu_time_check = 3600*atof(poptarg);// (Specify in hours)
00043                       break;
00044             case 'C': snap_cube_size = atof(poptarg);
00045                       break;
00046             case 'd': dt_out = atof(poptarg);
00047                       break;
00048             case 'D': dt_snap = atof(poptarg);
00049                       break;
00050             case 'e': init.ecc = atof(poptarg);
00051                       break;
00052             case 'g': init.tidal_tol_factor = atof(poptarg);
00053                       break;
00054             case 'L': init.r_init_min = atof(poptarg);
00055                       break;
00056             case 'm': init.m2 = atof(poptarg);
00057                       break;
00058             case 'M': init.m3 = atof(poptarg);
00059                       break;
00060             case 'n': n_experiments = atoi(poptarg);
00061                       break;
00062             case 'N': n_rand = atoi(poptarg);
00063                       break;
00064             case 'o': psi = atof(poptarg);
00065                       psi_flag = 1;
00066                       break;
00067             case 'p': planar_flag = 1;
00068                       break;
00069             case 'P': planar_flag = -1;
00070                       break;
00071             case 'q': q_flag = 1 - q_flag;
00072                       break;
00073             case 'Q': Q_flag = 1 - Q_flag;
00074                       break;
00075             case 'r': init.rho = atof(poptarg);
00076                       break;
00077             case 'R': init.r_stop = atof(poptarg);
00078                       init.r_init_min = init.r_init_max = abs(init.r_stop);
00079                       break;
00080             case 's': seed = atoi(poptarg);
00081                       break;
00082             case 'S': init.r_stop = atof(poptarg);
00083                       break;
00084             case 'U': init.r_init_max = atof(poptarg);
00085                       break;
00086             case 'v': init.v_inf = atof(poptarg);
00087                       break;
00088             case 'x': init.r1 = atof(poptarg);
00089                       break;
00090             case 'y': init.r2 = atof(poptarg);
00091                       break;
00092             case 'z': init.r3 = atof(poptarg);
00093                       break;
00094             case '?': params_to_usage(cerr, argv[0], param_string);
00095                       exit(1);
00096         }            
00097 
00098     if (Q_flag) q_flag = TRUE;
00099 
00100     if (init.m2 > 1) {
00101         cerr << "scatter3: init.m2 = " << init.m2 << " > 1" << endl;
00102         exit(1);
00103     }
00104 
00105     cpu_init();
00106     int random_seed = srandinter(seed, n_rand);
00107 
00108     for (int i = 0; i < n_experiments; i++) {
00109 
00110         if (n_experiments > 1) cerr << i+1 << ": ";
00111 
00112         cerr << "Random seed = " << get_initial_seed()
00113              << "  n_rand = " << get_n_rand() << endl << flush;
00114 
00115         // Normally, we just want random angles and phases.
00116         // However, in the planar case, it may be desireable to
00117         // control the orientation of the outer orbit.  The
00118         // angle between the inner and outer orbital axes in
00119         // the planar case is init.phase.psi.  Specify psi
00120         // in *degrees*, with psi = 0 meaning that the outer
00121         // periastron lies along the positive x-axis.
00122 
00123         randomize_angles(init.phase);
00124 
00125         if (planar_flag == 1) {
00126             init.phase.cos_theta = 1;   // Planar prograde
00127             if (psi_flag) init.phase.psi = psi * M_PI / 180.0;
00128         } else if (planar_flag == -1) {
00129             init.phase.cos_theta = -1;  // Planar retrograde
00130             if (psi_flag) init.phase.psi = psi * M_PI / 180.0;
00131         }
00132 
00133         intermediate_state3 inter;
00134         final_state3 final;
00135 
00136         real cpu = cpu_time();
00137 
00138         // scatter3(init, inter, final, cpu_time_check,
00139         //          dt_out, dt_snap, snap_cube_size);
00140 
00141         //-----------------------------------------------------------------
00142 
00143         // Perform scattering on a separate processor.
00144 
00145         char exe[100];
00146 
00147         strcpy(exe, getenv("STARLAB_PATH"));    // ASSUME that these
00148         if (exe[0] == '\0')                     // environment variables
00149             strcpy(exe, getenv("PWD"));         // are defined
00150         else
00151             strcat(exe, "/bin");
00152         strcat(exe, "/scatter3_slave.pvm");
00153 
00154         int tid;
00155         if (pvm_spawn(exe, (char**)0, 33, ".", 1, &tid) != 1) {
00156             cerr << "Unable to spawn slave task\n";
00157             exit(1);
00158         }
00159 
00160         if (pvm_recv(tid, HANDSHAKE_MSG) < 0) {
00161             cerr << "Error getting slave host ID\n";
00162             exit(1);
00163         }
00164 
00165         char name[100];
00166         pvm_upkstr(name);
00167 
00168         cerr << "initialized slave node " << name << endl;
00169         cerr << "tid = " << tid << endl;
00170 
00171         pvm_notify(PvmTaskExit, FAILURE_MSG, 1, &tid);
00172 
00173         // Pack the data.
00174 
00175         pvm_initsend(PvmDataRaw);
00176 
00177         int n_rand = get_initial_seed();
00178         pvm_pkint(&n_rand, 1, 1);
00179         pack(init);
00180 
00181         pvm_pkdouble(&cpu_time_check, 1, 1);
00182         pvm_pkdouble(&dt_snap, 1, 1);
00183         pvm_pkdouble(&snap_cube_size, 1, 1);
00184 
00185         // Send the data as message SEND_DATA_MSG.
00186 
00187         pvm_send(tid, SEND_DATA_MSG);
00188 
00189         struct timeval timeout;
00190         timeout.tv_sec = 3600;          // Give up after an hour.
00191         int bufid, length, msgid, tid_r;
00192 
00193         bufid = pvm_trecv(-1, -1, &timeout);
00194         pvm_bufinfo(bufid, &length, &msgid, &tid_r);
00195 
00196         cerr << "returned msgid = " << msgid << endl;
00197         if (msgid != RETURN_DATA_MSG) exit(0);
00198 
00199         pvm_upkint(&n_rand, 1, 1);
00200 
00201         // Note: n_rand should = random_seed[task]
00202 
00203         unpack(init);
00204         unpack(inter);
00205         unpack(final);
00206 
00207         int dum;
00208         pvm_upkint(&dum, 1, 1);
00209         pvm_upkdouble(&cpu, 1, 1);
00210 
00211         //-----------------------------------------------------------------
00212 
00213         // cpu = cpu_time() - cpu;
00214 
00215         cerr << ":  ";
00216         print_scatter3_outcome(inter, final, cerr);
00217 
00218         if (Q_flag) print_scatter3_summary(inter, final, cpu, cerr);
00219 
00220         if (!q_flag) print_scatter3_report(init, inter, final,
00221                                            cpu, b_flag, cerr);
00222     }
00223 }