Mercurial > hg > Members > kent > N-BodyProblem
comparison main_GL.cbc @ 0:249965d0a68f
GL is no finish.
author | kent |
---|---|
date | Thu, 29 May 2008 19:35:24 +0900 |
parents | |
children | 09e774f4433f |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:249965d0a68f |
---|---|
1 #include<stdio.h> | |
2 #include<stdlib.h> | |
3 #include<unistd.h> | |
4 #include<time.h> | |
5 #include<float.h> | |
6 #include<SDL.h> | |
7 #include<OpenGL/gl.h> | |
8 #include<OpenGL/glu.h> | |
9 | |
10 #define DEBUGlog(f, args...) \ | |
11 fprintf(stderr, "in %s: "f, __FUNCTION__, ## args) | |
12 | |
13 #define W_HEIGHT 480 | |
14 #define W_WIDTH 640 | |
15 | |
16 /* | |
17 N body problem | |
18 */ | |
19 | |
20 /* parameters */ | |
21 static int NUM_BODY = 3; | |
22 //static float Gravitation = 6.67e-11 ; | |
23 //static float delta = 100; | |
24 //static float FIELD = 2e11; | |
25 static float Gravitation = 0.2f; // ? | |
26 static float delta = 100.0f; // 0.01 ~ 100 ? | |
27 static float FIELD = 400.0f; // -100 ~ 100 | |
28 static const float eps = 0.0f; | |
29 | |
30 /* for OpenGL Utility. */ | |
31 GLUquadricObj **sphere; | |
32 | |
33 typedef struct | |
34 { | |
35 /* star's parameter. */ | |
36 float weight; | |
37 float a[3]; /* acceleration */ | |
38 float v[3]; /* velocity */ | |
39 float r[3]; /* position */ | |
40 /* for SDL. */ | |
41 //SDL_Rect rect; | |
42 } body; | |
43 body *stars_old, *stars_new; | |
44 | |
45 | |
46 void start(void); | |
47 __code initialize(int num); | |
48 __code randomInit(SDL_Surface *screen, int num); | |
49 __code starsInit(SDL_Surface *screen, int num); | |
50 __code loop(int count, SDL_Surface *screen, int num); | |
51 __code compute(int count, SDL_Surface *screen, int num); | |
52 __code nextTurn(int count, SDL_Surface *screen, int num); | |
53 __code moveCenter(int count, SDL_Surface *screen, int num); | |
54 __code CenteringVelocity(int count, SDL_Surface *screen, int num); | |
55 __code drawStars(int count, SDL_Surface *screen, int num); | |
56 | |
57 int main(int argc, char **argv) | |
58 { | |
59 int ch; | |
60 while ((ch = getopt(argc, argv, "s:g:")) != -1) { | |
61 switch (ch) { | |
62 case 's': | |
63 delta = atof(optarg); | |
64 break; | |
65 case 'g': | |
66 Gravitation = atof(optarg); | |
67 break; | |
68 case '?': | |
69 default: | |
70 break; | |
71 } | |
72 } | |
73 start(); | |
74 return 0; | |
75 } | |
76 void start() | |
77 { | |
78 goto initialize(NUM_BODY); | |
79 } | |
80 | |
81 __code finish(int ret, int num) | |
82 { | |
83 int i; | |
84 DEBUGlog("Gravity = %e\n", Gravitation); | |
85 DEBUGlog("FLT_MAX = %e\n", FLT_MAX); | |
86 DEBUGlog("FLT_MAX_EXP = %d\n", FLT_MAX_EXP); | |
87 DEBUGlog("FLT_MIN = %e\n", FLT_MIN); | |
88 DEBUGlog("FLT_MIN_EXP = %d\n", FLT_MIN_EXP); | |
89 DEBUGlog("FLT_EPSILON = %e\n", FLT_EPSILON); | |
90 free(stars_old); | |
91 free(stars_new); | |
92 for( i=0; i<NUM_BODY; i++){ | |
93 gluDeleteQuadric(sphere[i]); | |
94 } | |
95 free(sphere); | |
96 exit(ret); | |
97 } | |
98 | |
99 | |
100 __code initialize(int num) | |
101 { | |
102 SDL_Surface *screen; | |
103 int i; | |
104 | |
105 /* malloc. */ | |
106 stars_old = malloc( sizeof(body)*num ); | |
107 stars_new = malloc( sizeof(body)*num ); | |
108 sphere = malloc( sizeof(GLUquadricObj*)*num ); | |
109 if (stars_old==NULL||stars_new==NULL||sphere==NULL){ | |
110 perror("malloc"); | |
111 goto finish(1, num); | |
112 } | |
113 | |
114 /* SDL initialization. */ | |
115 if(SDL_Init(SDL_INIT_VIDEO) < 0){ //Could we start SDL_VIDEO? | |
116 fprintf(stderr,"Couldn't init SDL"); //Nope, output to stderr and quit | |
117 exit(1); | |
118 } | |
119 SDL_GL_SetAttribute( SDL_GL_RED_SIZE, 5); | |
120 SDL_GL_SetAttribute( SDL_GL_GREEN_SIZE, 5); | |
121 SDL_GL_SetAttribute( SDL_GL_BLUE_SIZE, 5); | |
122 SDL_GL_SetAttribute( SDL_GL_DEPTH_SIZE, 16); | |
123 SDL_GL_SetAttribute( SDL_GL_DOUBLEBUFFER, 1); | |
124 | |
125 screen = SDL_SetVideoMode(W_WIDTH, W_HEIGHT, 32, SDL_HWSURFACE | SDL_RESIZABLE | SDL_OPENGL ); //Create a 640x480x32 resizable window | |
126 atexit(SDL_Quit); //Now that we're enabled, make sure we cleanup | |
127 | |
128 /* initialize OpenGL. */ | |
129 glViewport(0,0,screen->w,screen->h); | |
130 glMatrixMode(GL_PROJECTION); | |
131 glLoadIdentity(); | |
132 gluPerspective( 60.0, (float)screen->w/(float)screen->h, 1.0, 1000.0); | |
133 gluLookAt( 500.0,500.0,500.0, 0.0,0.0,0.0, 1.0,0.0,0.0); | |
134 glClearColor(0.0, 0.0, 0.0, 0.0); | |
135 glMatrixMode(GL_MODELVIEW); | |
136 | |
137 //DEBUGlog("scr->w=%d\n", screen->w); | |
138 //DEBUGlog("scr->h=%d\n", screen->h); | |
139 for( i=0; i<num; i++){ | |
140 sphere[i] = gluNewQuadric(); | |
141 gluQuadricDrawStyle(sphere[i], GLU_FILL); | |
142 } | |
143 | |
144 goto starsInit(screen, num); | |
145 } | |
146 | |
147 __code starsInit0(SDL_Surface *screen, int num) | |
148 { | |
149 int i; | |
150 srandom(time(NULL)); | |
151 for (i=0; i<num; i++){ // this loop should be split into few code segment.. | |
152 stars_old[i].weight = random()/(RAND_MAX+1.0)*5+5; | |
153 stars_old[i].v[0] = random()/(RAND_MAX+1.0)*5+5; | |
154 stars_old[i].v[1] = random()/(RAND_MAX+1.0)*5+5; | |
155 stars_old[i].v[2] = random()/(RAND_MAX+1.0)*5+5; | |
156 stars_old[i].r[0] = random()/(RAND_MAX+1.0)*5+5; | |
157 stars_old[i].r[1] = random()/(RAND_MAX+1.0)*5+5; | |
158 stars_old[i].r[2] = random()/(RAND_MAX+1.0)*5+5; | |
159 stars_new[i].weight = stars_old[i].weight; | |
160 } | |
161 | |
162 goto loop(0, screen, num); | |
163 } | |
164 __code starsInit(SDL_Surface *screen, int num) | |
165 { | |
166 int i; | |
167 /* */ | |
168 stars_old[0].weight = 110; | |
169 stars_old[0].v[0] = 0.0; | |
170 stars_old[0].v[1] = -1.0; | |
171 stars_old[0].v[2] = 0.0; | |
172 stars_old[0].r[0] = 100.0; | |
173 stars_old[0].r[1] = 0.0; | |
174 stars_old[0].r[2] = 0.0; | |
175 /* */ | |
176 stars_old[1].weight = 110; | |
177 stars_old[1].v[0] = 0.0; | |
178 stars_old[1].v[1] = -1.0; | |
179 stars_old[1].v[2] = 0.0; | |
180 stars_old[1].r[0] = -100.0; | |
181 stars_old[1].r[1] = 0.0; | |
182 stars_old[1].r[2] = 0.0; | |
183 /* */ | |
184 stars_old[2].weight = 110; | |
185 stars_old[2].v[0] = -1.0; | |
186 stars_old[2].v[1] = 0.0; | |
187 stars_old[2].v[2] = 0.0; | |
188 stars_old[2].r[0] = 0.0; | |
189 stars_old[2].r[1] = 0.0; | |
190 stars_old[2].r[2] = -70.0; | |
191 | |
192 for( i=0; i<num; i++){ | |
193 stars_new[i].weight = stars_old[i].weight; | |
194 //stars_new[i].rect.h = 5, stars_new[i].rect.w = 5; | |
195 //stars_old[i].rect.h = 5, stars_old[i].rect.w = 5; | |
196 } | |
197 | |
198 goto loop(0, screen, num); | |
199 } | |
200 | |
201 __code starsInit1(SDL_Surface *screen, int num) | |
202 { | |
203 int i; | |
204 /* Sun */ | |
205 stars_old[0].weight = 1.9891e30; // 1.9891*10^30 | |
206 stars_old[0].v[0] = 0.0; | |
207 stars_old[0].v[1] = 0.0; | |
208 stars_old[0].v[2] = 0.0; | |
209 stars_old[0].r[0] = 0.0; | |
210 stars_old[0].r[1] = 0.0; | |
211 stars_old[0].r[2] = 0.0; | |
212 /* Venus */ | |
213 stars_old[1].weight = 4.869e24; // 4.869*10^24 | |
214 stars_old[1].v[0] = 0.0; | |
215 stars_old[1].v[1] = 3.50214e4; // 35.0214 km/s | |
216 stars_old[1].v[2] = 0.0; | |
217 stars_old[1].r[0] = 1.08e11; // 108,208,930 km | |
218 stars_old[1].r[1] = 0.0; | |
219 stars_old[1].r[2] = 0.0; | |
220 /* Earth */ | |
221 stars_old[2].weight = 5.9742e24; // 5.9742*10^24 | |
222 stars_old[2].v[0] = 0.0; | |
223 stars_old[2].v[1] = 2.97859e4; // 29.7859 km/s | |
224 stars_old[2].v[2] = 0.0; | |
225 stars_old[2].r[0] = 1.49e11; // 149,597,870km | |
226 stars_old[2].r[1] = 0.0; | |
227 stars_old[2].r[2] = 0.0; | |
228 | |
229 for( i=0; i<num; i++){ | |
230 stars_new[i].weight = stars_old[i].weight; | |
231 //stars_new[i].rect.h = 5, stars_new[i].rect.w = 5; | |
232 //stars_old[i].rect.h = 5, stars_old[i].rect.w = 5; | |
233 } | |
234 | |
235 goto loop(0, screen, num); | |
236 } | |
237 | |
238 __code loop(int count, SDL_Surface *screen, int num) | |
239 { | |
240 SDL_Event event; | |
241 | |
242 /* check SDL event. */ | |
243 while(SDL_PollEvent(&event)){ //Poll events | |
244 switch(event.type){ //Check event type | |
245 case SDL_QUIT: //User hit the X (or equivelent) | |
246 goto finish(1, num); | |
247 //case SDL_VIDEORESIZE: //User resized window | |
248 //screen = SDL_SetVideoMode(event.resize.w, event.resize.h, 32, | |
249 //SDL_HWSURFACE | SDL_RESIZABLE); // Create new window | |
250 //break; //Event handled, fetch next :) | |
251 case SDL_KEYDOWN: | |
252 if (event.key.keysym.sym==SDLK_UP && event.key.state==SDL_PRESSED) | |
253 FIELD *= 2.0; | |
254 else if (event.key.keysym.sym==SDLK_DOWN && event.key.state==SDL_PRESSED) | |
255 FIELD /= 2.0; | |
256 else if (event.key.keysym.sym==SDLK_r && event.key.state==SDL_PRESSED) | |
257 goto moveCenter(count, screen, num); | |
258 else if (event.key.keysym.sym==SDLK_v && event.key.state==SDL_PRESSED) | |
259 goto CenteringVelocity(count, screen, num); | |
260 else if (event.key.keysym.sym==SDLK_l && event.key.state==SDL_PRESSED) | |
261 goto AdjustLook(count, screen, num); | |
262 break; | |
263 default: | |
264 break; | |
265 } //Finished with current event | |
266 } //Done with all events for now | |
267 | |
268 if ( count<num ){ | |
269 DEBUGlog("count %d, goto commpute().\n", count); | |
270 goto compute(count, screen, num); | |
271 }else{ | |
272 DEBUGlog("count %d, goto nextTurn()\n", count); | |
273 goto nextTurn(count, screen, num); | |
274 } | |
275 } | |
276 | |
277 __code AdjustLook(int count, SDL_Surface *screen, int num) | |
278 { | |
279 int i; | |
280 float c0,c1,c2; | |
281 float v0,v1,v2,v; | |
282 c0=c1=c2=v0=v1=v2=0; | |
283 | |
284 for (i=0; i<num; i++){ | |
285 c0 += stars_new[i].r[0]; | |
286 c1 += stars_new[i].r[1]; | |
287 c2 += stars_new[i].r[2]; | |
288 } | |
289 c0/=3.0; c1/=3.0; c2/=3.0; | |
290 for (i=0; i<num; i++){ | |
291 v0 += (stars_new[i].r[0]-c0)*(stars_new[i].r[0]-c0); | |
292 v1 += (stars_new[i].r[1]-c1)*(stars_new[i].r[1]-c1); | |
293 v2 += (stars_new[i].r[2]-c2)*(stars_new[i].r[2]-c2); | |
294 } | |
295 v = v0+v1+v2; | |
296 v0 = sqrt(v0); v1 = sqrt(v1); v2 = sqrt(v2); | |
297 v = sqrt(v); | |
298 | |
299 glMatrixMode(GL_PROJECTION); | |
300 glLoadIdentity(); | |
301 gluPerspective( 60.0, (float)screen->w/(float)screen->h, 1.0, v); | |
302 gluLookAt( v,v,v, 0.0,0.0,0.0, 1.0,0.0,0.0); | |
303 glClearColor(0.0, 0.0, 0.0, 0.0); | |
304 glMatrixMode(GL_MODELVIEW); | |
305 | |
306 goto loop(0, screen, num); | |
307 } | |
308 | |
309 __code CenteringVelocity(int count, SDL_Surface *screen, int num) | |
310 { | |
311 int i; | |
312 float v0,v1,v2; | |
313 v0=v1=v2=0.0; | |
314 | |
315 for (i=0; i<num; i++){ | |
316 v0 += stars_old[i].v[0]; | |
317 v1 += stars_old[i].v[1]; | |
318 v2 += stars_old[i].v[2]; | |
319 } | |
320 v0/=(float)num; v1/=(float)num; v2/=(float)num; | |
321 for (i=0; i<num; i++){ | |
322 stars_old[i].v[0] -= v0; | |
323 stars_old[i].v[1] -= v1; | |
324 stars_old[i].v[2] -= v2; | |
325 } | |
326 | |
327 goto loop(0, screen, num); | |
328 } | |
329 __code moveCenter(int count, SDL_Surface *screen, int num) | |
330 { | |
331 int i; | |
332 float m0,m1,m2; | |
333 m0=m1=m2=0.0; | |
334 | |
335 for (i=0; i<num; i++){ | |
336 m0 += stars_old[i].r[0]; | |
337 m1 += stars_old[i].r[1]; | |
338 m2 += stars_old[i].r[2]; | |
339 } | |
340 m0/=(float)num; m1/=(float)num; m2/=(float)num; | |
341 for (i=0; i<num; i++){ | |
342 stars_old[i].r[0] -= m0; | |
343 stars_old[i].r[1] -= m1; | |
344 stars_old[i].r[2] -= m2; | |
345 } | |
346 | |
347 goto loop(0, screen, num); | |
348 } | |
349 | |
350 /* | |
351 * next x = x+dx | |
352 * dx = v*dt => dx/dt = v | |
353 * dv = a*dt => dv/dt = a | |
354 * a = F/m; | |
355 * F = F1 + F2 + ... + Fn | |
356 * Fi = G m1m2/r^2 | |
357 * | |
358 */ | |
359 __code compute(int count, SDL_Surface *screen, int num) | |
360 { | |
361 int i; | |
362 /* a is accel this planet receive now. */ | |
363 stars_old[count].a[0]=stars_old[count].a[1]=stars_old[count].a[2]=0.0f; | |
364 DEBUGlog("count=%d\n", count); | |
365 for (i=0; i<num; i++){ // this loop should be split into few code segment.. | |
366 //float F; | |
367 float d0, d1, d2, d; | |
368 float a; | |
369 //body *o, *m; | |
370 //o = &stars_old[i]; m = &stars_old[count]; | |
371 | |
372 /* skip compute with itself. */ | |
373 if ( i==count ) continue; | |
374 /* compute distance between two i-th planet and itself. */ | |
375 d0 = stars_old[i].r[0] - stars_old[count].r[0]; | |
376 d1 = stars_old[i].r[1] - stars_old[count].r[1]; | |
377 d2 = stars_old[i].r[2] - stars_old[count].r[2]; | |
378 d = ( d0*d0+d1*d1+d2*d2 ); | |
379 /* compute force it receive from i-th planet. */ | |
380 //F = Gravitation * stars_old[i].weight * stars_old[count].weight / d; | |
381 /* and accel. */ | |
382 //a = F/stars_old[count].weight; | |
383 a = Gravitation/d * stars_old[i].weight ; | |
384 stars_old[count].a[0] += a*d0/sqrt(d); | |
385 stars_old[count].a[1] += a*d1/sqrt(d); | |
386 stars_old[count].a[2] += a*d2/sqrt(d); | |
387 DEBUGlog("a=%e, d=%e, d0=%e\n", a, d, d0); | |
388 DEBUGlog("g*w=%e\n", Gravitation*stars_old[i].weight); | |
389 } | |
390 | |
391 stars_new[count].v[0] = stars_old[count].v[0]+stars_old[count].a[0]*delta; | |
392 stars_new[count].v[1] = stars_old[count].v[1]+stars_old[count].a[1]*delta; | |
393 stars_new[count].v[2] = stars_old[count].v[2]+stars_old[count].a[2]*delta; | |
394 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*/; | |
395 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*/; | |
396 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*/; | |
397 | |
398 DEBUGlog("x = %e\n", stars_new[count].r[0]); | |
399 DEBUGlog("y = %e\n", stars_new[count].r[1]); | |
400 DEBUGlog("z = %e\n", stars_new[count].r[2]); | |
401 DEBUGlog("vx = %e\n", stars_new[count].v[0]); | |
402 DEBUGlog("vy = %e\n", stars_new[count].v[1]); | |
403 DEBUGlog("vz = %e\n", stars_new[count].v[2]); | |
404 DEBUGlog("a0=%e,a1=%e,a2=%e\n", stars_new[count].a[0], stars_new[count].a[1], stars_new[count].a[2]); | |
405 | |
406 goto loop(++count, screen, num); | |
407 } | |
408 | |
409 /* it is executed after all bodies parameter is computed. */ | |
410 __code nextTurn(int count, SDL_Surface *screen, int num) | |
411 { | |
412 body *tmp; | |
413 tmp = stars_old; | |
414 stars_old = stars_new; | |
415 stars_new = tmp; | |
416 SDL_UpdateRect(screen, 0, 0, 0, 0); | |
417 | |
418 goto drawStars(count, screen, num); | |
419 } | |
420 | |
421 __code drawStars(int count, SDL_Surface *screen, int num) | |
422 { | |
423 int i; | |
424 DEBUGlog("draw\n"); | |
425 glClear( GL_COLOR_BUFFER_BIT ); | |
426 for (i=0; i<num; i++){ | |
427 glPushMatrix(); | |
428 glTranslatef( stars_new[i].r[0], stars_new[i].r[1], stars_new[i].r[2]); | |
429 gluSphere( sphere[i], 200.0, 8.0, 8.0 ); | |
430 glPopMatrix(); | |
431 } | |
432 | |
433 glColor3d( 1.0, 1.0, 1.0); | |
434 glBegin(GL_LINES); | |
435 glVertex3d( 0.0, 0.0, 0.0); | |
436 glVertex3d( 10.0, 0.0, 0.0); | |
437 glEnd(); | |
438 glBegin(GL_LINES); | |
439 glVertex3d( 0.0, 0.0, 0.0); | |
440 glVertex3d( 0.0, 10.0, 0.0); | |
441 glEnd(); | |
442 glBegin(GL_LINES); | |
443 glVertex3d( 0.0, 0.0, 0.0); | |
444 glVertex3d( 0.0, 0.0, 10.0); | |
445 glEnd(); | |
446 | |
447 SDL_GL_SwapBuffers(); | |
448 goto loop(0, screen, num); | |
449 } | |
450 |