- for (k = i + 1; k < gp->ngalaxies; ++k) {
- Galaxy *gtk = &gp->galaxies[k];
- double d0 = gtk->pos[0] - gt->pos[0];
- double d1 = gtk->pos[1] - gt->pos[1];
- double d2 = gtk->pos[2] - gt->pos[2];
-
- d = d0 * d0 + d1 * d1 + d2 * d2;
- if (d > EPSILON)
- d = gt->mass * gt->mass / (d * sqrt(d)) * DELTAT * QCONS;
- else
- d = gt->mass * gt->mass / (EPSILON * sqrt_EPSILON) * DELTAT * QCONS;
-
- d0 *= d;
- d1 *= d;
- d2 *= d;
- gt->vel[0] += d0 / gt->mass;
- gt->vel[1] += d1 / gt->mass;
- gt->vel[2] += d2 / gt->mass;
- gtk->vel[0] -= d0 / gtk->mass;
- gtk->vel[1] -= d1 / gtk->mass;
- gtk->vel[2] -= d2 / gtk->mass;
+ cox = COSF(gp->rot_y);
+ six = SINF(gp->rot_y);
+ cor = COSF(gp->rot_x);
+ sir = SINF(gp->rot_x);
+
+ eps = 1/(EPSILON * sqrt_EPSILON * DELTAT * DELTAT * QCONS);
+
+ for (i = 0; i < gp->ngalaxies; ++i) {
+ Galaxy *gt = &gp->galaxies[i];
+
+ for (j = 0; j < gp->galaxies[i].nstars; ++j) {
+ Star *st = >->stars[j];
+ XPoint *newp = >->newpoints[j];
+ double v0 = st->vel[0];
+ double v1 = st->vel[1];
+ double v2 = st->vel[2];
+
+ for (k = 0; k < gp->ngalaxies; ++k) {
+ Galaxy *gtk = &gp->galaxies[k];
+ double d0 = gtk->pos[0] - st->pos[0];
+ double d1 = gtk->pos[1] - st->pos[1];
+ double d2 = gtk->pos[2] - st->pos[2];
+
+ d = d0 * d0 + d1 * d1 + d2 * d2;
+ if (d > EPSILON)
+ d = gtk->mass / (d * sqrt(d)) * DELTAT * DELTAT * QCONS;
+ else
+ d = gtk->mass / (eps * sqrt(eps));
+ v0 += d0 * d;
+ v1 += d1 * d;
+ v2 += d2 * d;
+ }
+
+ st->vel[0] = v0;
+ st->vel[1] = v1;
+ st->vel[2] = v2;
+
+ st->pos[0] += v0;
+ st->pos[1] += v1;
+ st->pos[2] += v2;
+
+ newp->x = (short) (((cox * st->pos[0]) - (six * st->pos[2])) *
+ gp->scale) + gp->midx;
+ newp->y = (short) (((cor * st->pos[1]) - (sir * ((six * st->pos[0]) +
+ (cox * st->pos[2]))))
+ * gp->scale) + gp->midy;
+
+ }
+
+ for (k = i + 1; k < gp->ngalaxies; ++k) {
+ Galaxy *gtk = &gp->galaxies[k];
+ double d0 = gtk->pos[0] - gt->pos[0];
+ double d1 = gtk->pos[1] - gt->pos[1];
+ double d2 = gtk->pos[2] - gt->pos[2];
+
+ d = d0 * d0 + d1 * d1 + d2 * d2;
+ if (d > EPSILON)
+ d = 1 / (d * sqrt(d)) * DELTAT * QCONS;
+ else
+ d = 1 / (EPSILON * sqrt_EPSILON) * DELTAT * QCONS;
+
+ d0 *= d;
+ d1 *= d;
+ d2 *= d;
+ gt->vel[0] += d0 * gtk->mass;
+ gt->vel[1] += d1 * gtk->mass;
+ gt->vel[2] += d2 * gtk->mass;
+ gtk->vel[0] -= d0 * gt->mass;
+ gtk->vel[1] -= d1 * gt->mass;
+ gtk->vel[2] -= d2 * gt->mass;
+ }
+
+ gt->pos[0] += gt->vel[0] * DELTAT;
+ gt->pos[1] += gt->vel[1] * DELTAT;
+ gt->pos[2] += gt->vel[2] * DELTAT;
+
+ if (dbufp) {
+ XSetForeground(display, gc, MI_WIN_BLACK_PIXEL(mi));
+ XDrawPoints(display, window, gc, gt->oldpoints, gt->nstars,
+ CoordModeOrigin);
+ }
+ XSetForeground(display, gc, MI_PIXEL(mi, COLORSTEP * gt->galcol));
+ XDrawPoints(display, window, gc, gt->newpoints, gt->nstars,
+ CoordModeOrigin);
+
+ dummy = gt->oldpoints;
+ gt->oldpoints = gt->newpoints;
+ gt->newpoints = dummy;