#include<stdio.h>
#include<stdlib.h>
#include<unistd.h>
#include<time.h>
#include<SDL.h>
#include<float.h>

#define DEBUGlog(f, args...) \
	;
	//fprintf(stderr, "in %s: "f, __FUNCTION__, ## args)

#define W_HEIGHT 480
#define W_WIDTH 640

/*
	N body problem
*/

static int NUM_BODY = 3;
//static float Gravitation = 6.67e-11 ;
//static float delta = 100;
//static float FIELD = 2e11;
static float Gravitation = 1.0f; // ?
static float delta = 0.05f; // 
static float FIELD = 400.0f; // -100 ~ 100
static const float eps = 1.5f;

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 moveCenter(int count, SDL_Surface *screen, int num);
__code CenteringVelocity(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)
{
	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);
	exit(ret);
}


__code initialize(int num)
{
	SDL_Surface *screen;

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

	/* 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);
	}
	screen = SDL_SetVideoMode(W_WIDTH, W_HEIGHT, 32, SDL_HWSURFACE | SDL_RESIZABLE); //Create a 640x480x32 resizable window
	atexit(SDL_Quit); //Now that we're enabled, make sure we cleanup

	goto starsInit(screen, num);
}

__code starsInit0(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;
#if 0
	/*  */
	stars_old[0].weight = 110;
	stars_old[0].v[0] = 0.0;
	stars_old[0].v[1] = -1.0;
	stars_old[0].v[2] = 0.0;
	stars_old[0].r[0] = 100.0;
	stars_old[0].r[1] = 0.0;
	stars_old[0].r[2] = 0.0;
	/*  */
	stars_old[1].weight = 110;
	stars_old[1].v[0] = 0.0;
	stars_old[1].v[1] = -1.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 = 110;
	stars_old[2].v[0] = -1.0;
	stars_old[2].v[1] = 0.0;
	stars_old[2].v[2] = 0.0;
	stars_old[2].r[0] = 0.0;
	stars_old[2].r[1] = 0.0;
	stars_old[2].r[2] = -70.0;
#else
	/*  */
	stars_old[0].weight = 100;
	stars_old[0].v[0] = -1.0;
	stars_old[0].v[1] = 0.0;
	stars_old[0].v[2] = 0.5;
	stars_old[0].r[0] = 50.0;
	stars_old[0].r[1] = 0.0;
	stars_old[0].r[2] = 0.0;
	/*  */
	stars_old[1].weight = 100;
	stars_old[1].v[0] = -0.5;
	stars_old[1].v[1] = 0.1;
	stars_old[1].v[2] = 1.732;
	stars_old[1].r[0] = -50.0;
	stars_old[1].r[1] = 0.0;
	stars_old[1].r[2] = 0.0;
	/*  */
	stars_old[2].weight = 100;
	stars_old[2].v[0] = 0.5;
	stars_old[2].v[1] = -0.1;
	stars_old[2].v[2] = -1.732;
	stars_old[2].r[0] = 0.0;
	stars_old[2].r[1] = 0.0;
	stars_old[2].r[2] = 86.60;
#endif

	for( i=0; i<num; i++){
		stars_new[i].weight = stars_old[i].weight;
		stars_new[i].rect.h = 5, stars_new[i].rect.w = 5;
		stars_old[i].rect.h = 5, stars_old[i].rect.w = 5;
	}

	goto loop(0, screen, num);
}

__code starsInit1(SDL_Surface *screen, int num)
{
	int i;
	/* Sun */
	stars_old[0].weight = 1.9891e30; // 1.9891*10^30
	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;
	/* Venus  */
	stars_old[1].weight = 4.869e24;  // 4.869*10^24
	stars_old[1].v[0] = 0.0;
	stars_old[1].v[1] = 3.50214e4;  // 35.0214 km/s
	stars_old[1].v[2] = 0.0;
	stars_old[1].r[0] = 1.08e11; // 108,208,930 km
	stars_old[1].r[1] = 0.0;
	stars_old[1].r[2] = 0.0;
	/* Earth  */
	stars_old[2].weight = 5.9742e24;  // 5.9742*10^24
	stars_old[2].v[0] = 0.0;
	stars_old[2].v[1] = 2.97859e4;  // 29.7859 km/s
	stars_old[2].v[2] = 0.0;
	stars_old[2].r[0] = 1.49e11; // 149,597,870km
	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;
		stars_new[i].rect.h = 5, stars_new[i].rect.w = 5;
		stars_old[i].rect.h = 5, stars_old[i].rect.w = 5;
	}

	goto loop(0, screen, num);
}

__code loop(int count, SDL_Surface *screen, int num)
{
	SDL_Event event;

	/* 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);
				//case SDL_VIDEORESIZE: //User resized window
				//screen = SDL_SetVideoMode(event.resize.w, event.resize.h, 32,
				//SDL_HWSURFACE | SDL_RESIZABLE); // Create new window
				//break; //Event handled, fetch next :)
			case SDL_KEYDOWN:
				if (event.key.keysym.sym==SDLK_UP && event.key.state==SDL_PRESSED)
					FIELD *= 2.0;
				else if (event.key.keysym.sym==SDLK_DOWN && event.key.state==SDL_PRESSED)
					FIELD /= 2.0;
				else if (event.key.keysym.sym==SDLK_r && event.key.state==SDL_PRESSED)
					goto moveCenter(count, screen, num);
				else if (event.key.keysym.sym==SDLK_v && event.key.state==SDL_PRESSED)
					goto CenteringVelocity(count, screen, num);
				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 CenteringVelocity(int count, SDL_Surface *screen, int num)
{		
	int i;	
	float v0,v1,v2;
	v0=v1=v2=0.0;
		
	for (i=0; i<num; i++){
		v0 += stars_old[i].v[0];
		v1 += stars_old[i].v[1];
		v2 += stars_old[i].v[2];
	}	
	v0/=(float)num; v1/=(float)num; v2/=(float)num;
	for (i=0; i<num; i++){
		stars_old[i].v[0] -= v0;
		stars_old[i].v[1] -= v1;
		stars_old[i].v[2] -= v2;
	}	
		
	goto loop(0, screen, num);
}		
__code moveCenter(int count, SDL_Surface *screen, int num)
{		
	int i;	
	float m0,m1,m2;
	m0=m1=m2=0.0;
		
	for (i=0; i<num; i++){
		m0 += stars_old[i].r[0];
		m1 += stars_old[i].r[1];
		m2 += stars_old[i].r[2];
	}	
	m0/=(float)num; m1/=(float)num; m2/=(float)num;
	for (i=0; i<num; i++){
		stars_old[i].r[0] -= m0;
		stars_old[i].r[1] -= m1;
		stars_old[i].r[2] -= m2;
	}	
		
	goto loop(0, screen, num);
}		
		
/*		
 *  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
 *		
 */		
__code compute(int count, SDL_Surface *screen, int num)
{		
	int i;	
	/* 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);
	for (i=0; i<num; i++){ // this loop should be split to few code segment..
		//float F;
		float d0, d1, d2, d;
		float a;
		//body *o, *m;
		//o = &stars_old[i]; m = &stars_old[count];
		
		/* skip compute with itself.  */
		if ( i==count ) continue;
		/* 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/sqrt(d);
		stars_old[count].a[1] += a*d1/sqrt(d);
		stars_old[count].a[2] += a*d2/sqrt(d);
		DEBUGlog("a=%e, d=%e, d0=%e\n", a, d, d0);
		DEBUGlog("g*w=%e\n", Gravitation*stars_old[i].weight);
	}	

	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*/;

	stars_new[count].rect.x = (stars_new[count].r[0]*(float)W_HEIGHT/FIELD/2)+(W_WIDTH/2);
	stars_new[count].rect.y = (stars_new[count].r[1]*(float)W_HEIGHT/FIELD/2)+(W_HEIGHT/2);
	DEBUGlog("x = %e, x=%d\n", stars_new[count].r[0], stars_new[count].rect.x);
	DEBUGlog("y = %e, y=%d\n", stars_new[count].r[1], stars_new[count].rect.y);
	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]);
	DEBUGlog("x=%d,y=%d,h=%d,w=%d\n",stars_old[count].rect.x,stars_old[count].rect.y,stars_old[count].rect.h,stars_old[count].rect.w);
	SDL_FillRect(screen, &stars_old[count].rect, SDL_MapRGB(screen->format, 0x00,0x00,0x00));
	SDL_FillRect(screen, &stars_new[count].rect, SDL_MapRGB(screen->format, 0xff,0x88,0x33));
	//SDL_BlitSurface(screen, &stars_old[count].rect, screen, &stars_new[count].rect);

	goto loop(++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 loop(0, screen, num);
}		
		
	       	       	     
	 		     
/*
       +--------+-----------+
       |        v           |
       |  +-----------+	    |
       |  | loop      |	    |
       |  +-----+-----+	    |
       |       	|      	    |
       |  +-----v-----+	    |
       |  + compute   +-----+
       |  +-----+-----+ count++
       |       	|
       |       	| count==num
       |  +-----v-----+	   
       |  | nextTurn  |	       
       |  +-----+-----+	       
       |        | count=0      
       +--------+      	    
*/