#include<stdio.h>
#include<stdlib.h>
#include<unistd.h>
#include<time.h>
#include<float.h>
#include<math.h>
#include<SDL.h>
#include<OpenGL/gl.h>
#include<OpenGL/glu.h>

#ifdef DEBUG
  #define DEBUGlog(f, args...) \
	fprintf(stderr, "in %s: "f, __FUNCTION__, ## args)
#else
  #define DEBUGlog(f, args...) ;
#endif

#define W_HEIGHT 480
#define W_WIDTH 640

/*
	N body problem
*/

/* parameters  */
static int NUM_BODY = 3;
static float Gravitation = 100.0f; // ?
static float delta = 0.05f; // 0.01 ~ 100  ?
static const float eps = 15.0f;
/*
 * だいたい、平均距離をn倍にするんだったら、Gravitationoはn^2倍、
 * epsはn倍
 */

/* for OpenGL Utility.  */
GLUquadricObj **sphere;
float DefDistance = 300.0f;

typedef struct
{
	/* star's parameter.  */
	float weight;
	float a[3]; /* acceleration */
	float v[3]; /* velocity */
	float r[3]; /* position */
	/* for SDL.  */
	//SDL_Rect rect;
} body;
body *stars_old, *stars_new;


void start(void);
__code initialize(int num);
__code randomInit(SDL_Surface *screen, int num);
__code starsInit(SDL_Surface *screen, int num);
__code loop(int count, SDL_Surface *screen, int num);
__code compute(int count, SDL_Surface *screen, int num);
__code nextTurn(int count, SDL_Surface *screen, int num);
__code drawStars(SDL_Surface *screen, int num);
__code AdjustLooking(float distance, int count, SDL_Surface *screen, int num);

int main(int argc, char **argv)
{
	int ch;
	while ((ch = getopt(argc, argv, "s:g:")) != -1) {
		switch (ch) {
			case 's':
				delta = atof(optarg);
				break;
			case 'g':
				Gravitation = atof(optarg);
				break;
			case '?':
			default:
				break;
		}
	}
	start();
	return 0;
}
void start()
{
	goto initialize(NUM_BODY);
}

__code finish(int ret, int num)
{
	int i;
	DEBUGlog("Gravity = %e\n", Gravitation);
	DEBUGlog("FLT_MAX = %e\n", FLT_MAX);
	DEBUGlog("FLT_MAX_EXP = %d\n", FLT_MAX_EXP);
	DEBUGlog("FLT_MIN = %e\n", FLT_MIN);
	DEBUGlog("FLT_MIN_EXP = %d\n", FLT_MIN_EXP);
	DEBUGlog("FLT_EPSILON = %e\n", FLT_EPSILON);
	free(stars_old);
	free(stars_new);
	for( i=0; i<NUM_BODY; i++){
		gluDeleteQuadric(sphere[i]);
	}
	free(sphere);
	exit(ret);
}


__code initialize(int num)
{
	SDL_Surface *screen;
	int i;

	/* malloc.  */
	stars_old = malloc( sizeof(body)*num );
	stars_new = malloc( sizeof(body)*num );
	sphere = malloc( sizeof(GLUquadricObj*)*num );
	if (stars_old==NULL||stars_new==NULL||sphere==NULL){
		perror("malloc");
		goto finish(1, num);
	}

	/* SDL initialization.  */
	if(SDL_Init(SDL_INIT_VIDEO) < 0){ //Could we start SDL_VIDEO?
		fprintf(stderr,"Couldn't init SDL"); //Nope, output to stderr and quit 
		exit(1);
	}
    SDL_GL_SetAttribute( SDL_GL_RED_SIZE, 5);
    SDL_GL_SetAttribute( SDL_GL_GREEN_SIZE, 5);
    SDL_GL_SetAttribute( SDL_GL_BLUE_SIZE, 5);
    SDL_GL_SetAttribute( SDL_GL_DEPTH_SIZE, 16);
    SDL_GL_SetAttribute( SDL_GL_DOUBLEBUFFER, 1);

	screen = SDL_SetVideoMode(W_WIDTH, W_HEIGHT, 32, SDL_HWSURFACE | SDL_RESIZABLE | SDL_OPENGL ); //Create a 640x480x32 resizable window
	atexit(SDL_Quit); //Now that we're enabled, make sure we cleanup

	/* initialize OpenGL.  */
	glViewport(0,0,screen->w,screen->h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective( 60.0, (float)screen->w/(float)screen->h, 1.0, 1000.0);
	gluLookAt( 300.0,300.0,300.0, 0.0,0.0,0.0, 1.0,0.0,0.0);
	glClearColor(0.0, 0.0, 0.0, 0.0);
	glMatrixMode(GL_MODELVIEW);

	for( i=0; i<num; i++){
		sphere[i] = gluNewQuadric();
		gluQuadricDrawStyle(sphere[i], GLU_FILL);
	}

	goto starsInit(screen, num);
}

__code starsInitRandom(SDL_Surface *screen, int num)
{
	int i;
	srandom(time(NULL));
	for (i=0; i<num; i++){ // this loop should be split into few code segment..
		stars_old[i].weight = random()/(RAND_MAX+1.0)*5+5;
		stars_old[i].v[0] = random()/(RAND_MAX+1.0)*5+5;
		stars_old[i].v[1] = random()/(RAND_MAX+1.0)*5+5;
		stars_old[i].v[2] = random()/(RAND_MAX+1.0)*5+5;
		stars_old[i].r[0] = random()/(RAND_MAX+1.0)*5+5;
		stars_old[i].r[1] = random()/(RAND_MAX+1.0)*5+5;
		stars_old[i].r[2] = random()/(RAND_MAX+1.0)*5+5;
		stars_new[i].weight = stars_old[i].weight;
	}

	goto loop(0, screen, num);
}
__code starsInit(SDL_Surface *screen, int num)
{
	int i;
	/*  */
	stars_old[0].weight = 1000;
	stars_old[0].v[0] = 0.0;
	stars_old[0].v[1] = 0.0;
	stars_old[0].v[2] = 0.0;
	stars_old[0].r[0] = 0.0;
	stars_old[0].r[1] = 0.0;
	stars_old[0].r[2] = 0.0;
	/*  */
	stars_old[1].weight = 5;
	stars_old[1].v[0] = 0.1;
	stars_old[1].v[1] = 5.0;
	stars_old[1].v[2] = 0.0;
	stars_old[1].r[0] = 100.0;
	stars_old[1].r[1] = 0.0;
	stars_old[1].r[2] = 0.0;
	/*  */
	stars_old[2].weight = 5;
	stars_old[2].v[0] = 0.0;
	stars_old[2].v[1] = -5.0;
	stars_old[2].v[2] = 0.1;
	stars_old[2].r[0] = -100.0;
	stars_old[2].r[1] = 0.0;
	stars_old[2].r[2] = 0.0;

	for( i=0; i<num; i++){
		stars_new[i].weight = stars_old[i].weight;
	}

	goto loop(0, screen, num);
}


/* MAIN LOOP ROUTIN
 *
  +------+
  | loop |<-----------------------+
  +---+--+                        |
      +-----------+               |
      |count=num  | count<num     |
  +---v------+ +--v-------------+ |
  | nextTurn | | COMPUTE ROUTIN | |
  +---+------+ +-------+--------+ |
      |                |          |
  +---v---------+      v          |
  | DRAW ROUTIN +----->+----------+
  +-------------+    back to loop
 */
__code loop(int count, SDL_Surface *screen, int num)
{
	SDL_Event event;

	usleep(1000);
	/* check SDL event.  */
	while(SDL_PollEvent(&event)){ //Poll events
		switch(event.type){ //Check event type
			case SDL_QUIT: //User hit the X (or equivelent)
				goto finish(1, num);
			case SDL_KEYDOWN:
				if (event.key.state==SDL_PRESSED){
					switch(event.key.keysym.sym){
						case SDLK_UP:
							goto AdjustLooking(DefDistance/=1.2, count, screen, num);
						case SDLK_DOWN:
							goto AdjustLooking(DefDistance*=1.2, count, screen, num);
						case SDLK_l:
							goto AdjustLooking(DefDistance, count, screen, num);
						case SDLK_q:
							goto finish(0, num);
						default:
							break;
					}
				}
				break;
			default:
				break;
		} //Finished with current event
	} //Done with all events for now

	if ( count<num ){
		DEBUGlog("count %d, goto commpute().\n", count);
		goto compute(count, screen, num);
	}else{
		DEBUGlog("count %d, goto nextTurn()\n", count);
		goto nextTurn(count, screen, num);
	}
}

__code AdjustLooking(float distance, int count, SDL_Surface *screen, int num)
{
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective( 60.0, (float)screen->w/(float)screen->h, 1.0, distance*2.5);
	gluLookAt( distance,distance,distance, 0.0,0.0,0.0, 1.0,0.0,0.0);
	glClearColor(0.0, 0.0, 0.0, 0.0);
	glMatrixMode(GL_MODELVIEW);

	goto loop(0, screen, num);
}

/*
 * COMPUTING ROUTINs
 *
 *  next x = x+dx
 *  dx = v*dt  =>  dx/dt = v
 *  dv = a*dt  =>  dv/dt = a
 *  a = F/m;
 *  F = F1 + F2 + ... + Fn
 *  Fi = G m1m2/r^2
        |
  +-----v-----+
  |  compute  |      
  +-----+-----+      
        |<----------+
  +-----v---------+ | i++
  |  computeForce +-+
  +-----+---------+ i<num
        | i==num
  +-----+---------+
  |  com|ute_out  |
  +-----+---------+
        v to loop
 */
__code compute_out(int count, SDL_Surface *screen, int num);
__code computeForce(int i, int count, SDL_Surface *screen, int num);
__code compute(int count, SDL_Surface *screen, int num)
{
	/* a is accel this planet receive now.  */
	stars_old[count].a[0]=stars_old[count].a[1]=stars_old[count].a[2]=0.0f;
	DEBUGlog("count=%d\n", count);

	goto computeForce( 0, count, screen, num);
}
__code computeForce(int i, int count, SDL_Surface *screen, int num)
{
	//float F;
	float d0, d1, d2, d;
	float a;
	//body *o, *m;
	//o = &stars_old[i]; m = &stars_old[count];
	
	if (i>=num) goto compute_out(count, screen, num);
	/* skip compute with itself.  */
	else if (i==count) goto computeForce(i+1, count, screen, num);

	/* compute distance between two i-th planet and itself.  */
	d0 = stars_old[i].r[0] - stars_old[count].r[0];
	d1 = stars_old[i].r[1] - stars_old[count].r[1];
	d2 = stars_old[i].r[2] - stars_old[count].r[2];
	d  = ( d0*d0+d1*d1+d2*d2+eps*eps );
	/* compute force it receive from i-th planet.  */
	//F  = Gravitation * stars_old[i].weight * stars_old[count].weight / d;
	/* and accel.  */
	//a = F/stars_old[count].weight;
	a  = Gravitation/d * stars_old[i].weight ;
	stars_old[count].a[0] += a*d0/sqrtf(d);
	stars_old[count].a[1] += a*d1/sqrtf(d);
	stars_old[count].a[2] += a*d2/sqrtf(d);
	DEBUGlog("a=%e, d=%e, d0=%e\n", a, d, d0);
	DEBUGlog("g*w=%e\n", Gravitation*stars_old[i].weight);

	goto computeForce( i+1, count, screen, num);
}
__code compute_out(int count, SDL_Surface *screen, int num)
{

	stars_new[count].v[0] = stars_old[count].v[0]+stars_old[count].a[0]*delta;
	stars_new[count].v[1] = stars_old[count].v[1]+stars_old[count].a[1]*delta;
	stars_new[count].v[2] = stars_old[count].v[2]+stars_old[count].a[2]*delta;
	stars_new[count].r[0] = stars_old[count].r[0]+stars_new[count].v[0]*delta/*+stars_old[count].a[0]*delta*delta/2.0*/;
	stars_new[count].r[1] = stars_old[count].r[1]+stars_new[count].v[1]*delta/*+stars_old[count].a[1]*delta*delta/2.0*/;
	stars_new[count].r[2] = stars_old[count].r[2]+stars_new[count].v[2]*delta/*+stars_old[count].a[2]*delta*delta/2.0*/;

	DEBUGlog("x = %e\n", stars_new[count].r[0]);
	DEBUGlog("y = %e\n", stars_new[count].r[1]);
	DEBUGlog("z = %e\n", stars_new[count].r[2]);
	DEBUGlog("vx = %e\n", stars_new[count].v[0]);
	DEBUGlog("vy = %e\n", stars_new[count].v[1]);
	DEBUGlog("vz = %e\n", stars_new[count].v[2]);
	DEBUGlog("a0=%e,a1=%e,a2=%e\n", stars_new[count].a[0], stars_new[count].a[1], stars_new[count].a[2]);

	goto loop(1+count, screen, num);
}

/* it is executed after all bodies parameter is computed.  */
__code nextTurn(int count, SDL_Surface *screen, int num)
{
	body *tmp;
	tmp = stars_old;
	stars_old = stars_new;
	stars_new = tmp;
	//SDL_UpdateRect(screen, 0, 0, 0, 0);

	goto drawStars(screen, num);
}



/*
 * DRAWING ROUTINs
 *
       |
  +----+------+  
  | drawStars |   
  +----+------+   
       |<--------+
  +----v-----+   | count++
  | drawStar +---+
  +----+-----+  count<num
       | count==num
  +----v----------+
  | drawStars_out |
  +----+----------+
       | to loop
 */
__code drawStar(int i, SDL_Surface *screen, int num);
__code drawStars_out(SDL_Surface *screen, int num);

__code drawStars(SDL_Surface *screen, int num)
{
	DEBUGlog("draw\n");
	glClear( GL_COLOR_BUFFER_BIT );
	glColor3d( 1.0, 1.0, 1.0);
	glBegin(GL_LINES);
	glVertex3d( 0.0, 0.0, 0.0);
	glVertex3d( 200.0, 0.0, 0.0);
	glEnd();
	glBegin(GL_LINES);
	glVertex3d( 0.0, 0.0, 0.0);
	glVertex3d( 0.0, 200.0, 0.0);
	glEnd();
	glBegin(GL_LINES);
	glVertex3d( 0.0, 0.0, 0.0);
	glVertex3d( 0.0, 0.0, 200.0);
	glEnd();

	goto drawStar(0, screen, num);
}

__code drawStars_out(SDL_Surface *screen, int num)
{
	DEBUGlog("draw axis: screen=%p, num=%d\n", screen, num);
	SDL_GL_SwapBuffers();
	DEBUGlog("go loop\n");
	goto loop(0, screen, num);
}

__code drawStar(int i, SDL_Surface *screen, int num)
{
	if (i==num)
		goto drawStars_out(screen, num);

	DEBUGlog("i=%d, num=%d\n", i, num);
	glPushMatrix();
	glTranslatef( stars_new[i].r[0], stars_new[i].r[1], stars_new[i].r[2]);
	gluSphere( sphere[i], 2.0, 32.0, 32.0 );
	glPopMatrix();

	DEBUGlog("to myself\n");
	goto drawStar(i+1, screen, num);
}