plot_stars.C Source File

00001 
00002 // plot_stars:  plot part of an N-body system.
00003 
00004 //.............................................................................
00005 //    version 1:  Steve, Jun 2000, adapted from "starplot"
00006 //.............................................................................
00007 //  Plot positions of some bodies in an N-body system, in a primitive way:
00008 //  the screen is divided in character positions, and the representation is:
00009 //  every character stands for one body; every * for more than one body.
00010 //  Very primitive, but independent of terminal or plotting software.
00011 //.............................................................................
00012 
00013 #include "dyn.h"
00014 
00015 #ifndef TOOLBOX
00016 
00017 local real find_neighbors(dyn *bi, dyn **list, int n)   // note: this n is n+1!
00018 {
00019     dyn *root = bi->get_root();
00020 
00021     real *dr2 = new real[n];
00022     for (int i = 0; i < n; i++) {
00023         list[i] = NULL;
00024         dr2[i] = VERY_LARGE_NUMBER;
00025     }
00026 
00027     for_all_daughters(dyn, root, bj) {
00028 
00029         real sep2 = square(bj->get_pos() - bi->get_pos());
00030         int j;
00031 
00032         for (j = n-1; j >= 0; j--)
00033             if (sep2 > dr2[j]) break;
00034 
00035         // bj should be in location j+1
00036 
00037         if (j < n-1) {
00038             for (int k = n-1; k > j+1; k--) {
00039                 list[k] = list[k-1];
00040                 dr2[k] = dr2[k-1];
00041             }
00042             list[j+1] = bj;
00043             dr2[j+1] = sep2;
00044         }
00045     }
00046 
00047     for (int i = n-1; i >= 0; i--)
00048         if (list[i]) return sqrt(dr2[i]);
00049 
00050     return 0;
00051 }
00052 
00053 #define  HBINS  41                      // number of horizontal bins
00054 #define  VBINS  21                      // number of vertical bins
00055 
00056 void plot_stars(dyn * bi,
00057                 int n,                  // default = 5
00058                 int k)                  // default = 3: plot (x, y)
00059 
00060 // Draw the n stars closest to bi.
00061 
00062 {
00063     if (n <= 0) return;
00064     bi = bi->get_top_level_node();      // deal only with top-level nodes
00065 
00066     // Make a list of the n (actual, not projected) nearest neighbors of bi.
00067     // Convenient to keep bi at location 0.
00068 
00069     dyn **list = new dynptr[n+1];
00070     real scale = find_neighbors(bi, list, n+1);
00071     if (scale <= 0) return;
00072 
00073     // Clear the display and make a box.
00074 
00075     char disp[HBINS][VBINS];            // pixels on the display
00076 
00077     for (int j = 0; j < VBINS; j++)
00078         for (int i = 0; i < HBINS; i++)
00079             disp[i][j] = ' ';
00080 
00081     for (int j = 0; j < VBINS; j++)
00082         disp[0][j] = disp[HBINS-1][j] = '|';
00083 
00084     for (int i = 0; i < HBINS; i++)
00085         disp[i][0] = disp[i][VBINS-1] = '-';
00086 
00087     disp[0][0] = disp[0][VBINS-1]
00088         = disp[HBINS-1][0] = disp[HBINS-1][VBINS-1] = '+';
00089 
00090     // Add the stars.
00091 
00092     int  kx, ky;                        // projected coordinates
00093 
00094     switch (k) {
00095         case 1: kx = 1; ky = 2; break;  // project along x
00096         case 2: kx = 2; ky = 0; break;  // project along y
00097         default:
00098         case 3: kx = 0; ky = 1; break;  // project along z
00099     }
00100 
00101     // Refine scale before proceeding.
00102 
00103     scale = 0;
00104     for (int m = 0; m <= n; m++) {
00105         if (list[m]) {
00106             real x = list[m]->get_pos()[kx] - list[0]->get_pos()[kx];
00107             real y = list[m]->get_pos()[ky] - list[0]->get_pos()[ky];
00108             if (abs(x) > scale) scale = abs(x);
00109             if (abs(y) > scale) scale = abs(y);
00110         }
00111     }
00112 
00113     real scale2 = 256;
00114     while (scale2 > scale) scale2 /= sqrt(2.0);
00115     scale2 *= sqrt(2.0);
00116     scale = scale2;
00117 
00118     // Create the plot.
00119 
00120     for (int m = 0; m <= n; m++) {
00121         if (list[m]) {
00122 
00123             // Center of mass.
00124 
00125             real x = list[m]->get_pos()[kx] - list[0]->get_pos()[kx];
00126             real y = list[m]->get_pos()[ky] - list[0]->get_pos()[ky];
00127             int i = (int)(0.5*(0.999999 + x/scale) * HBINS);
00128             int j = (int)(0.5*(0.999999 + y/scale) * VBINS);
00129             if (disp[i][j] == ' '
00130                 || i == 0 || i == HBINS-1
00131                 || j == 0 || j == VBINS-1) {
00132                 if (m < 10)
00133                     disp[i][j] = '0' + m;
00134                 else
00135                     disp[i][j] = 'a' + m - 10;
00136             } else
00137                 disp[i][j] = '*';
00138 
00139             // Add next-level components, if any.
00140 
00141             dyn *od = list[m]->get_oldest_daughter();
00142             if (od) {
00143                 x += od->get_pos()[kx];
00144                 y += od->get_pos()[ky];
00145                 int icm = i, jcm = j;
00146                 i = (int)(0.5*(0.999999 + x/scale) * HBINS);
00147                 j = (int)(0.5*(0.999999 + y/scale) * VBINS);
00148                 if (disp[i][j] == ' '
00149                      || (disp[i][j] != '*' && i == icm && j == jcm))
00150                     disp[i][j] = 'X';
00151                 else
00152                     disp[i][j] = '*';
00153 
00154                 dyn *yd = od->get_younger_sister();
00155                 x += yd->get_pos()[kx] - od->get_pos()[kx];
00156                 y += yd->get_pos()[ky] - od->get_pos()[ky];
00157                 i = (int)(0.5*(0.999999 + x/scale) * HBINS);
00158                 j = (int)(0.5*(0.999999 + y/scale) * VBINS);
00159                 if (disp[i][j] == ' ')
00160                     disp[i][j] = 'x';
00161                 else
00162                     disp[i][j] = '*';
00163             }
00164         }
00165     }
00166 
00167     // Print the results.
00168 
00169     int off = 10;
00170     cerr << endl;
00171     PRI(off+1); cerr << "system_time = " << bi->get_system_time() << endl;
00172     PRI(off+1); cerr << "plot scale  = " << 2*scale << endl << endl;
00173 
00174     int m = 0;
00175     for (int j = VBINS-1; j >= 0; j--) {
00176         PRI(off);
00177         for (int i = 0; i < HBINS; i++) cerr << disp[i][j];
00178         if (j < VBINS-1 && m <= n) {
00179             fprintf(stderr, "   %3d: %7s %10.2e",
00180                     m, list[m]->format_label(), list[m++]->get_mass()); 
00181             // cerr << "     " << m << ":  " << list[m]->format_label()
00182             //      << list[m++]->get_mass();
00183         }
00184         cerr << endl;
00185     }
00186     cerr << endl;
00187 
00188     // Add an array of all interparticle separations if not too large.
00189 
00190     if (n <= 5) {
00191 
00192         // Make a list of all nodes and leaves involved.
00193 
00194         int nnodes = 1000;
00195 
00196         while (nnodes > 10) {
00197           nnodes = 0;
00198           for (m = 0; m <= n; m++) {
00199             if (list[m]) {
00200               for_all_nodes(dyn, list[m], bb) nnodes++;
00201             }
00202           }
00203           if (nnodes > 10) n--;
00204           if (n <= 0) break;
00205         }
00206 
00207         if (n > 0) {
00208 
00209           PRL(nnodes);
00210 
00211           dyn ** nodes = new dynptr[nnodes];
00212 
00213           int i = 0;
00214           for (m = 0; m <= n; m++) {
00215             if (list[m]) {
00216               for_all_nodes(dyn, list[m], bb) nodes[i++] = bb;
00217             }
00218           }
00219 
00220           PRI(14);
00221           for (i = 0; i < nnodes; i++) {
00222 
00223             // Label placement is Steve's aesthetic judgement (6/00)...
00224 
00225             int len = strlen(nodes[i]->format_label());
00226             if (len > 9) len = 9;
00227             int ind = 7-len;
00228             if (len > 3) ind = 6 - (1+len)/2;
00229 
00230             PRI(ind);
00231             fprintf(stderr, "%*s", len, nodes[i]->format_label());
00232             PRI(10-ind-len);
00233           }
00234           cerr << endl;
00235 
00236           for (int j = 0; j < nnodes; j++) {
00237             fprintf(stderr, "%13s ", nodes[j]->format_label());
00238             for (i = 0; i < nnodes; i++)
00239                 if (i == j)
00240                     fprintf(stderr, "          ");
00241                 else
00242                     fprintf(stderr, "%10.2e",
00243                       abs(something_relative_to_root(nodes[i], &dyn::get_pos) -
00244                           something_relative_to_root(nodes[j], &dyn::get_pos)));
00245             cerr << endl;
00246           }
00247           cerr << endl;
00248         }
00249     }
00250 }
00251 
00252 #else
00253 
00254 main(int argc, char ** argv)
00255 {
00256     int  k = 3;
00257     int n = 5;
00258     int seed = 0;
00259     dyn *b;
00260 
00261     check_help();
00262 
00263     extern char *poptarg;
00264     int c;
00265     char* param_string = "a:n:s:";
00266 
00267     while ((c = pgetopt(argc, argv, param_string)) != -1)
00268         switch(c) {
00269             case 'a': k = atoi(poptarg);
00270                       break;
00271             case 'n': n = atoi(poptarg);
00272                       break;
00273             case 's': seed = atoi(poptarg);
00274                       break;
00275             case '?': params_to_usage(cerr, argv[0], param_string);
00276                       get_help();
00277                       exit(1);
00278     }            
00279 
00280     srandinter(seed);
00281 
00282     while (b = get_dyn(cin)) {
00283         dyn *bi = b->get_oldest_daughter();
00284         int nd = (int)(randinter(0, 1) * b->n_daughters());
00285         for (int i = 0; i < nd-1; i++) bi = bi->get_younger_sister();
00286         plot_stars(bi, n, k);
00287         rmtree(b);
00288     }
00289 }
00290 
00291 #endif