Mercurial > hg > Members > kent > N-BodyProblem

changeset 0:249965d0a68f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
GL is no finish.
author kent
date Thu, 29 May 2008 19:35:24 +0900
parents
children 09e774f4433f
files Makefile main.cbc main_GL.cbc
diffstat 3 files changed, 883 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Makefile	Thu May 29 19:35:24 2008 +0900
@@ -0,0 +1,27 @@
+
+CbCC = ~/work/cvs_tmp/CbC_project/build-dir/INSTALL-DIR/bin/gcc
+.SUFFIXES: .cbc .o
+
+CFLAGS = -Wall -g -O2
+INCLUDES = `sdl-config --cflags`
+LIBS = `sdl-config --libs`
+
+GL_LIBS = -Wl,-framework,OpenGL
+
+TARGETS = main main_GL
+
+%.o: %.c
+
+all: $(TARGETS)
+
+.cbc.o:
+	$(CbCC) -c $(INCLUDES) $(CFLAGS) $^ -o $@
+
+main: main.o
+	$(CbCC) $(LIBS) $(CFLAGS) $^ -o $@
+main_GL: main_GL.o
+	$(CbCC) $(LIBS) $(GL_LIBS) $(CFLAGS) $^ -o $@
+
+clean:
+	rm -f *.o *.s $(TARGETS)
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main.cbc	Thu May 29 19:35:24 2008 +0900
@@ -0,0 +1,406 @@
+#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      
+       +--------+      	    
+*/		
+		
+		
+		
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main_GL.cbc	Thu May 29 19:35:24 2008 +0900
@@ -0,0 +1,450 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include<unistd.h>
+#include<time.h>
+#include<float.h>
+#include<SDL.h>
+#include<OpenGL/gl.h>
+#include<OpenGL/glu.h>
+
+#define DEBUGlog(f, args...) \
+	fprintf(stderr, "in %s: "f, __FUNCTION__, ## args)
+
+#define W_HEIGHT 480
+#define W_WIDTH 640
+
+/*
+	N body problem
+*/
+
+/* parameters  */
+static int NUM_BODY = 3;
+//static float Gravitation = 6.67e-11 ;
+//static float delta = 100;
+//static float FIELD = 2e11;
+static float Gravitation = 0.2f; // ?
+static float delta = 100.0f; // 0.01 ~ 100  ?
+static float FIELD = 400.0f; // -100 ~ 100
+static const float eps = 0.0f;
+
+/* for OpenGL Utility.  */
+GLUquadricObj **sphere;
+
+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);
+__code drawStars(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( 500.0,500.0,500.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);
+
+	//DEBUGlog("scr->w=%d\n", screen->w);
+	//DEBUGlog("scr->h=%d\n", screen->h);
+	for( i=0; i<num; i++){
+		sphere[i] = gluNewQuadric();
+		gluQuadricDrawStyle(sphere[i], GLU_FILL);
+	}
+
+	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;
+	/*  */
+	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;
+
+	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, num);
+				//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);
+				else if (event.key.keysym.sym==SDLK_l && event.key.state==SDL_PRESSED)
+					goto AdjustLook(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 AdjustLook(int count, SDL_Surface *screen, int num)
+{
+	int i;
+	float c0,c1,c2;
+	float v0,v1,v2,v;
+	c0=c1=c2=v0=v1=v2=0;
+
+	for (i=0; i<num; i++){
+		c0 += stars_new[i].r[0];
+		c1 += stars_new[i].r[1];
+		c2 += stars_new[i].r[2];
+	}
+	c0/=3.0; c1/=3.0; c2/=3.0; 
+	for (i=0; i<num; i++){
+		v0 += (stars_new[i].r[0]-c0)*(stars_new[i].r[0]-c0);
+		v1 += (stars_new[i].r[1]-c1)*(stars_new[i].r[1]-c1);
+		v2 += (stars_new[i].r[2]-c2)*(stars_new[i].r[2]-c2);
+	}
+	v = v0+v1+v2;
+	v0 = sqrt(v0); v1 = sqrt(v1); v2 = sqrt(v2);
+	v = sqrt(v);
+
+	glMatrixMode(GL_PROJECTION);
+	glLoadIdentity();
+	gluPerspective( 60.0, (float)screen->w/(float)screen->h, 1.0, v);
+	gluLookAt( v,v,v, 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);
+}
+
+__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 into 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 );
+		/* 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*/;
+
+	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(++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(count, screen, num);
+}
+
+__code drawStars(int count, SDL_Surface *screen, int num)
+{
+	int i;
+	DEBUGlog("draw\n");
+	glClear( GL_COLOR_BUFFER_BIT );
+	for (i=0; i<num; i++){
+		glPushMatrix();
+		glTranslatef( stars_new[i].r[0], stars_new[i].r[1], stars_new[i].r[2]);
+		gluSphere( sphere[i], 200.0, 8.0, 8.0 );
+		glPopMatrix();
+	}
+
+	glColor3d( 1.0, 1.0, 1.0);
+	glBegin(GL_LINES);
+	glVertex3d( 0.0, 0.0, 0.0);
+	glVertex3d( 10.0, 0.0, 0.0);
+	glEnd();
+	glBegin(GL_LINES);
+	glVertex3d( 0.0, 0.0, 0.0);
+	glVertex3d( 0.0, 10.0, 0.0);
+	glEnd();
+	glBegin(GL_LINES);
+	glVertex3d( 0.0, 0.0, 0.0);
+	glVertex3d( 0.0, 0.0, 10.0);
+	glEnd();
+
+	SDL_GL_SwapBuffers();
+	goto loop(0, screen, num);
+}
+