kira_sample_debug.C Source File

00001 
00002 local void dump_node(hdyn* b)
00003 {
00004     // Incomplete dump of hdyn data.
00005 
00006     cerr << endl;
00007     int p = cerr.precision(HIGH_PRECISION);
00008 
00009     PRL(b->format_label());
00010     PRL(b->get_time());
00011     PRL(b->get_timestep());
00012     PRL(b->get_unperturbed_timestep());
00013     PRL(b->get_pos());
00014     PRL(b->get_vel());
00015     PRL(b->get_acc());
00016     PRL(b->get_jerk());
00017     if (b->get_kepler())
00018         b->get_kepler()->print_all(cerr);
00019 
00020     cerr.precision(p);
00021 }
00022 
00023 local void dump_binaries(hdyn* b)
00024 {
00025     for_all_daughters(hdyn, b, bb)
00026         if (bb->get_oldest_daughter())
00027             for_all_nodes(hdyn, bb, bbb)
00028                 dump_node(bbb);
00029 }
00030 
00031 
00032 // Sample debugging lines to precede integrate_list() in evolve_system.
00033 
00034 
00035 #if 0
00036         for (int i = 0; i < n_next; i++) {
00037             hdyn* n = next_nodes[i];
00038             if (n->is_parent() || n->is_low_level_leaf()) {
00039                 PRL(t);
00040                 n->print_label(cerr);
00041                 dump_node(n);
00042             }
00043         }
00044 #endif
00045 
00046 #if 0
00047         for (int i = 0; i < n_next; i++) {
00048             if (next_nodes[i]) {
00049                 hdyn* bi = next_nodes[i];
00050                 if (bi->get_top_level_node()->is_parent()) {
00051                     cerr << endl;
00052                     cerr << "evolve_system: integration list contains "
00053                          << bi->format_label() << endl;
00054                 }
00055             }
00056         }
00057 #endif
00058 
00059 #if 0
00060         cerr << endl << "time = " << b->get_system_time() << endl;
00061 
00062         for (int j = 0; j < n_next; j++) {
00063             hdyn* n = next_nodes[j];
00064             cerr << n->format_label() << "  ";
00065         }
00066         cerr << endl;
00067 #endif
00068 
00069 //------------------------------------------------------------------------
00070 
00071 // Sample debugging lines to follow integrate_list() in evolve_system.
00072 
00073 #if 0
00074         for (int j = 0; j < n_next; j++) {
00075             hdyn* n = next_nodes[j];
00076             if (n && n->is_valid()) {
00077                 if (node_contains(n->get_top_level_node(), "13a")) {
00078                     if (n->is_top_level_node())
00079                         count0++;
00080                     hdyn* p = n->get_parent();
00081                     if (p->is_top_level_node())
00082                         count1++;
00083                     hdyn* pp = p->get_parent();
00084                     if (pp && pp->is_top_level_node())
00085                         count2++;
00086                     if ((count0+count1+count2)%1000 == 0) {
00087                         cerr << "counts:  " << count0 << "  " << count1
00088                              << "  " << count2 << "  " << count_ << "  ";
00089                         PRC(steps); PRL(cpu_time());
00090                     }
00091                 } else
00092                     count_++;
00093             }
00094         }
00095 #endif
00096 
00097 
00098 #if 0
00099         if (t > 18.999) {
00100             cerr << endl; PRC(t); PRL(n_next);
00101             for (int j = 0; j < n_next; j++) {
00102                 hdyn* n = next_nodes[j];
00103                 if (n && n->is_valid()) {
00104                     cerr << n->format_label() << ":  ";
00105                     if (n->is_low_level_node())
00106                         n->print_perturber_list(cerr);
00107                     pp3(n, cerr);
00108                 }
00109             }
00110             print_recalculated_energies(b, eps * eps);
00111         }
00112         if (t > 19.1) exit(0);
00113 #endif
00114 
00115 
00116 #if 0
00117         if (t > 0.46) {
00118 
00119             cerr << endl, PRC(t);
00120             cerr << "evolve_system: after integrate_list:";
00121 
00122             for (int i = 0; i < n_next; i++) {
00123                 hdyn* bi = next_nodes[i];
00124                 cerr << endl << "    ";
00125                 cerr << bi->format_label();
00126                 if (bi->is_top_level_node())
00127                     cerr << " (t";
00128                 else
00129                     cerr << " (l";
00130                 if (bi->is_leaf())
00131                     cerr << "l";
00132                 else
00133                     cerr << "n";
00134                 if (bi->get_kepler())
00135                     cerr << "u";
00136                 else
00137                     cerr << "p";
00138                 cerr << "), nn = ";
00139                 if (bi->get_nn())
00140                     cerr << bi->get_nn()->format_label();
00141                 else
00142                     cerr << "(nil)";
00143                 if (bi->get_top_level_node()->is_parent()) {
00144                     cerr << endl;
00145                     pp3(bi->get_top_level_node(), cerr);
00146                 }
00147             }
00148             cerr << endl;
00149             print_recalculated_energies(b, eps * eps);
00150 
00151         } else {
00152 
00153             for (int i = 0; i < n_next; i++) {
00154                 if (next_nodes[i]) {
00155                     hdyn* bi = next_nodes[i];
00156                     if (bi->get_top_level_node()->is_parent()) {
00157                         cerr << endl;
00158                         cerr << "evolve_system: pp3 triggered by "
00159                              << bi->format_label() << endl;
00160                         pp3(bi->get_top_level_node(), cerr);
00161                     }
00162                 }
00163             }
00164 
00165         }
00166 #endif
00167 
00168 
00169 #if 0
00170         for (int i = 0; i < n_next; i++) {
00171             if (next_nodes[i] != NULL) {
00172                 char c[128];
00173                 real dt = next_nodes[i]->get_timestep();
00174                 int idt = (int)(log(1.0/dt)/log(2.0) + .1);
00175                 sprintf(c, "time = %.14f, new timestep = %.14f = 2^-%d",
00176                         t, dt, idt);
00177                 next_nodes[i]->log_comment(c);
00178             }
00179         }
00180 #endif
00181 
00182 #if 0
00183         if (t > 2.61) {
00184             for_all_nodes(hdyn, b, bb) {
00185                 if (bb != b && bb->get_oldest_daughter()
00186                     && !bb->get_kepler()) {
00187                     if (!bb->get_valid_perturbers()) {
00188                         PRC(t); cerr << "  no valid perturbers for "
00189                                      << bb->format_label() << endl;
00190                     } else if (bb->get_n_perturbers() > 0) {
00191                         for (int ip = 0; ip < bb->get_n_perturbers(); ip++) {
00192                             hdyn* p = bb->get_perturber_list()[ip];
00193                             if (!p) {
00194                                 PRC(t); cerr << "  NULL perturber #" << ip
00195                                     << " for " << bb->format_label() << endl;
00196                             } else if (!p->is_valid()) {
00197                                 PRC(t); cerr << "  invalid perturber #" << ip
00198                                     << " for " << bb->format_label() << endl;
00199                             }
00200                         }                       
00201                     }
00202                 }
00203             }
00204         }
00205 
00206         if (t > 28.1679) {
00207             cerr << endl;
00208             PRL(t);
00209             for (int j = 0; j < n_next; j++) {
00210                 hdyn* n = next_nodes[j];
00211                 PRC(j);
00212                 if (n) pp3(n);
00213                 else cerr << endl;
00214             }
00215         }
00216 
00217         if (t > 2.67) {
00218             bool print_energy = false;
00219             for (int j = 0; j < n_next; j++) {
00220                 hdyn* n = next_nodes[j];
00221                 if (n && n->is_valid()
00222                     && n->is_top_level_node())
00223                     print_energy = true;
00224             }
00225             if (print_energy)
00226                 print_recalculated_energies(b, eps * eps);
00227         }
00228 
00229 
00230         if (t > 1.3 && t < 1.375) {
00231             for (int j = 0; j < n_next; j++) {
00232                 hdyn* n = next_nodes[j];
00233                 if (n && n->is_valid()
00234                     && n->name_is("(1376b,1376a)")) {
00235                     cerr << "node " << n->format_label() << "  ";
00236                     PRL(t);
00237                     dump_story(n->get_dyn_story());
00238                 }
00239             }
00240         }
00241 
00242 
00243 
00244         if (t > 566.246119)
00245             pp3((hdyn*)node_with_name("(1840,2812)", b->get_root()));
00246 
00247         if (fmod(count, 100) == 0) {
00248 
00249             // Primitive checksum...
00250 
00251             real checkx = 0, checkv = 0, checka = 0, checkj = 0;
00252             for_all_nodes(hdyn, b, bb) {
00253                 checkx += abs(bb->get_pos());
00254                 checkv += abs(bb->get_vel());
00255                 checka += abs(bb->get_acc());
00256                 checkj += abs(bb->get_jerk());
00257             }
00258             int p = cerr.precision(HIGH_PRECISION);
00259             cerr << endl; PRL(t);
00260             PRC(checkx); PRL(checkv);
00261             PRC(checka); PRL(checkj);
00262             cerr.precision(p);
00263         }
00264 
00265         if (t > 566.246127) exit(0);
00266                 
00267 
00268 //      if (b->get_system_time() > 566.8)
00269             for (int j = 0; j < n_next; j++) {
00270                 hdyn* n = next_nodes[j];
00271                 if (n && n->is_valid()
00272                     && n->is_low_level_node()
00273                     && fmod(n->get_steps(), 100.0) == 0
00274                     && n->name_is("1139")
00275                     && n->get_perturbation_squared() < 0.0001) {
00276 
00277                     pp3(n);
00278 
00279                     real dist = abs(n->get_pos()
00280                                     - n->get_younger_sister()->get_pos());
00281                     real dtff = sqrt(dist*dist*dist
00282                                      / (2*n->get_parent()->get_mass()));
00283                     real dtv  = sqrt(square(n->get_pos())
00284                                      / square(n->get_vel()));
00285 
00286                     real kepstep = eta*min(dtff, dtv);
00287                     real ratio = n->get_timestep()/kepstep;
00288                     PRC(kepstep); PRL(ratio);
00289                 }
00290             }
00291 #endif