/* standard defines -- it's better to have them here than
  * having to s/X/Y/g the whole thing over and over again...
  * andrey
  */
#include <u.h>
#include <libc.h>
#include <draw.h>
#include <event.h>

#define NULL nil
#define XPoint Point
#define SINF sin
#define COSF cos
#define NRAND nrand
#define LRAND lrand
#define MAXRAND ((2<<31)-1)
#define FLOATRAND frand()
#define M_PI 3.141592653589793238462643383 /* should be enough, no? */
/* end of plan9-related defines */


#define MINSIZE       1
#define MINGALAXIES    2
#define MAX_STARS    1000
#define MAX_IDELTAT    50
/* These come originally from the Cluster-version */
#define DEFAULT_GALAXIES  3
#define DEFAULT_STARS    1000
#define DEFAULT_HITITERATIONS  7500
#define DEFAULT_IDELTAT    200 /* 0.02 */
#define EPSILON 0.00000001

#define sqrt_EPSILON 0.0001

#define DELTAT (MAX_IDELTAT * 0.0001)

#define GALAXYRANGESIZE  0.1
#define GALAXYMINSIZE  0.15
#define QCONS    0.001


#define COLORBASE  16
  /* Colors for stars start here */
#define COLORSTEP  (NUMCOLORS/COLORBASE)
/* NUMCOLORS / COLORBASE colors per galaxy */


typedef struct {
 double      pos[3], vel[3];
} Star;


typedef struct {
 int         mass;
 int         nstars;
 Star       *stars;
 XPoint     *oldpoints;
 XPoint     *newpoints;
 double      pos[3], vel[3];
 int         galcol;
} Galaxy;

typedef struct {
 double      mat[3][3]; /* Movement of stars(?) */
 double      scale; /* Scale */
 int         midx; /* Middle of screen, x */
 int         midy; /* Middle of screen, y */
 double      size; /* */
 double      diff[3]; /* */
 Galaxy     *galaxies; /* the Whole Universe */
 int         ngalaxies; /* # galaxies */
 int         f_hititerations; /* # iterations before restart */
 int         step; /* */
 double      rot_y; /* rotation of eye around center of universe, around
y-axis*/
 double      rot_x; /* rotation of eye around center of universe, around
x-axis */
} unistruct;

static unistruct *universes = NULL;

static void
free_galaxies(unistruct * gp)
{
 if (gp->galaxies != NULL) {
  int         i;

  for (i = 0; i < gp->ngalaxies; i++) {
   Galaxy     *gt = &gp->galaxies[i];

   if (gt->stars != NULL)
    (void) free((void *) gt->stars);
   if (gt->oldpoints != NULL)
       (void) free((void *) gt->oldpoints);
   if (gt->newpoints != NULL)
       (void) free((void *) gt->newpoints);
  }
  (void) free((void *) gp->galaxies);
  gp->galaxies = NULL;
 }
}

static void
startover(void)
{
 unistruct  *gp = &universes[0];
 int         i, j; /* more tmp */
 double      w1, w2; /* more tmp */
 double      d, v, w, h; /* yet more tmp */

 gp->step = 0;
 gp->rot_y = 0;
 gp->rot_x = 0;

  gp->ngalaxies = MINGALAXIES + (int)(FLOATRAND+0.5);
 if (gp->galaxies == NULL)
  gp->galaxies = (Galaxy *) calloc(gp->ngalaxies, sizeof (Galaxy));

 for (i = 0; i < gp->ngalaxies; ++i) {
  Galaxy     *gt = &gp->galaxies[i];
  double      sinw1, sinw2, cosw1, cosw2;

  gt->galcol = NRAND(COLORBASE - 2);
  if (gt->galcol > 1)
   gt->galcol += 2; /* Mult 8; 16..31 no green stars */
  /* Galaxies still may have some green stars but are not all green. */

 if (gt->stars != NULL) {
   (void) free((void *) gt->stars);
   gt->stars = NULL;
  }
  gt->nstars = (NRAND(MAX_STARS / 2)) + MAX_STARS / 2;
  gt->stars = (Star *) malloc(gt->nstars * sizeof (Star));
  gt->oldpoints = (XPoint *) malloc(gt->nstars * sizeof (XPoint));
  gt->newpoints = (XPoint *) malloc(gt->nstars * sizeof (XPoint));

  w1 = 2.0 * M_PI * FLOATRAND;
  w2 = 2.0 * M_PI * FLOATRAND;
  sinw1 = SINF(w1);
  sinw2 = SINF(w2);
  cosw1 = COSF(w1);
  cosw2 = COSF(w2);

  gp->mat[0][0] = cosw2;
  gp->mat[0][1] = -sinw1 * sinw2;
  gp->mat[0][2] = cosw1 * sinw2;
  gp->mat[1][0] = 0.0;
  gp->mat[1][1] = cosw1;
  gp->mat[1][2] = sinw1;
  gp->mat[2][0] = -sinw2;
  gp->mat[2][1] = -sinw1 * cosw2;
  gp->mat[2][2] = cosw1 * cosw2;

  gt->vel[0] = FLOATRAND * 2.0 - 1.0;
  gt->vel[1] = FLOATRAND * 2.0 - 1.0;
  gt->vel[2] = FLOATRAND * 2.0 - 1.0;
  gt->pos[0] = -gt->vel[0] * DELTAT * gp->f_hititerations + FLOATRAND -
0.5;
  gt->pos[1] = -gt->vel[1] * DELTAT * gp->f_hititerations + FLOATRAND -
0.5;
  gt->pos[2] = -gt->vel[2] * DELTAT * gp->f_hititerations + FLOATRAND -
0.5;

 gt->mass = (int) (FLOATRAND * 1000.0) + 1;

  gp->size = GALAXYRANGESIZE * FLOATRAND;
 gp->size  +=  GALAXYMINSIZE;
 for (j = 0; j < gt->nstars; ++j) {
   Star       *st = &gt->stars[j];
   XPoint     *oldp = &gt->oldpoints[j];
   XPoint     *newp = &gt->newpoints[j];

   double      sinw, cosw;

   w = 2.0 * M_PI * FLOATRAND;
   sinw = SINF(w);
   cosw = COSF(w);
   d = FLOATRAND * gp->size;
   h = FLOATRAND * exp(-2.0 * (d / gp->size)) / 5.0 * gp->size;
   if (FLOATRAND < 0.5)
    h = -h;
   st->pos[0] = gp->mat[0][0] * d * cosw + gp->mat[1][0] * d * sinw +
gp->mat[2][0] * h + gt->pos[0];
   st->pos[1] = gp->mat[0][1] * d * cosw + gp->mat[1][1] * d * sinw +
gp->mat[2][1] * h + gt->pos[1];
   st->pos[2] = gp->mat[0][2] * d * cosw + gp->mat[1][2] * d * sinw +
gp->mat[2][2] * h + gt->pos[2];

   v = sqrt(gt->mass * QCONS / sqrt(d * d + h * h));
   st->vel[0] = -gp->mat[0][0] * v * sinw + gp->mat[1][0] * v * cosw +
gt->vel[0];
   st->vel[1] = -gp->mat[0][1] * v * sinw + gp->mat[1][1] * v * cosw +
gt->vel[1];
   st->vel[2] = -gp->mat[0][2] * v * sinw + gp->mat[1][2] * v * cosw +
gt->vel[2];

   st->vel[0] *= DELTAT;
   st->vel[1] *= DELTAT;
   st->vel[2] *= DELTAT;

   oldp->x = 0;
   oldp->y = 0;
   newp->x = 0;
   newp->y = 0;
  }

 }

 draw(screen, screen->r, display->black, nil, ZP);

}

void
init_galaxy(void)
{
 unistruct  *gp;

 if (universes == NULL) {
  if ((universes = (unistruct *) calloc(1, sizeof (unistruct))) == NULL)
   return;
 }
 gp = &universes[0];

 gp->f_hititerations = 100;

 gp->scale = (double) (Dx(screen->r) + Dy(screen->r)) / 8.0;
 gp->midx =  screen->r.min.x + Dx(screen->r)  / 2;
 gp->midy =  screen->r.min.y +Dy(screen->r) / 2;
startover();
}

void
draw_galaxy(void)
{
 unistruct  *gp = &universes[0];
 double      d, eps, cox, six, cor, sir;  /* tmp */
 int         i, j, k; /* more tmp */
 int 	     i2;
    XPoint    *dummy = NULL;

//draw(screen, screen->r, display->black, nil, ZP);

gp->rot_y += 0.01;
 gp->rot_x += 0.004;

 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 = &gt->stars[j];
   XPoint     *newp = &gt->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 = gt->mass / (d * sqrt(d)) * DELTAT * DELTAT * QCONS;
    else
        d = gt->mass * 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) (((st->pos[0]) - (st->pos[2])) * gp->scale) + gp->midx;
   newp->y = (short) (((st->pos[1]) - (((st->pos[0]) + (st->pos[2])))) * gp->scale) + gp->midy;
//   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 = 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;
  }

  gt->pos[0] += gt->vel[0] * DELTAT;
  gt->pos[1] += gt->vel[1] * DELTAT;
  gt->pos[2] += gt->vel[2] * DELTAT;


/* original drawing routines */
//         	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);

		for(i2 = 0; i2 < gt->nstars; i2++) {
			Rectangle r;
			Point p = Pt(1, 1);

			p = gt->oldpoints[i2];
			r = Rpt(p, addpt(p, Pt(2, 2)));
			draw(screen, r, display->black, nil, ZP);

			p = gt->newpoints[i2];
			r = Rpt(p, addpt(p, Pt(2, 2)));
			draw(screen, r, display->white, nil, ZP);

		}
	dummy = gt->oldpoints;
	gt->oldpoints = gt->newpoints;
	gt->newpoints = dummy;
 }

 gp->step++;
 if (gp->step > gp->f_hititerations * 4)
  startover();
}

void
release_galaxy(void)
{
 if (universes != NULL) {

   free_galaxies(&universes[0]);
  (void) free((void *) universes);
  universes = NULL;
 }
}

void
refresh_galaxy(void)
{
 /* Do nothing, it will refresh by itself */
}

void
eresized(int new)
{

	if(new && getwindow(display, Refnone) < 0) {
		sysfatal("can't reattach to window");
	}

	draw(screen, screen->r, display->black, nil, ZP);

	if(universes != nil)
		release_galaxy();
	init_galaxy();
}

char *buttons[] =
{
	"exit",
	0
};

Menu menu =
{
	buttons
};

void
main(int argc, char **argv)
{
	Mouse m;
	vlong tbegin, tend;
	ulong frames;

	srand(time(0));

	if(initdraw(nil, nil, "bez") < 0)
		sysfatal("initdraw failed: %r");

	einit(Emouse);

	eresized(0);
	tbegin = nsec();
	for(;;) {
		if(ecanmouse()) {
			m = emouse();
			if(m.buttons&4)
				if(emenuhit(3, &m, &menu) == 0) {
					tend = nsec();
					print("fps: %uf\n", frames/((tend - tbegin)/1000000000.0));
					exits(0);
				}
		}
		draw_galaxy();
		frames++;
	}
}--upas-kojbavveocoghhvaxhnnsdccem--