Mercurial > hg > Members > kent > N-BodyProblem
view main_GL.cbc @ 0:249965d0a68f
GL is no finish.
author | kent |
---|---|
date | Thu, 29 May 2008 19:35:24 +0900 |
parents | |
children | 09e774f4433f |
line wrap: on
line source
#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); }