Mercurial > hg > Members > kent > N-BodyProblem

comparison main_GL.cbc @ 3:01290d71ef9c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
add figures into comment area of main_GL.cbc
author kent
date Fri, 30 May 2008 12:36:59 +0900
parents 3e543e31b6eb
children 5535e1d5e693
comparison
equal deleted inserted replaced
2:3e543e31b6eb 3:01290d71ef9c
1 #include<stdio.h> 1 #include<stdio.h>
2 #include<stdlib.h> 2 #include<stdlib.h>
3 #include<unistd.h> 3 #include<unistd.h>
4 #include<time.h> 4 #include<time.h>
5 #include<float.h> 5 #include<float.h>
6 #include<math.h>
6 #include<SDL.h> 7 #include<SDL.h>
7 #include<OpenGL/gl.h> 8 #include<OpenGL/gl.h>
8 #include<OpenGL/glu.h> 9 #include<OpenGL/glu.h>
9 10
10 #define DEBUGlog(f, args...) \ 11 #ifdef DEBUG
11 ; 12 #define DEBUGlog(f, args...) \
12 //fprintf(stderr, "in %s: "f, __FUNCTION__, ## args) 13 fprintf(stderr, "in %s: "f, __FUNCTION__, ## args)
14 #else
15 #define DEBUGlog(f, args...) ;
16 #endif
13 17
14 #define W_HEIGHT 480 18 #define W_HEIGHT 480
15 #define W_WIDTH 640 19 #define W_WIDTH 640
16 20
17 /* 21 /*
18 N body problem 22 N body problem
19 */ 23 */
20 24
21 /* parameters */ 25 /* parameters */
22 static int NUM_BODY = 3; 26 static int NUM_BODY = 3;
23 //static float Gravitation = 6.67e-11 ;
24 //static float delta = 100;
25 //static float FIELD = 2e11;
26 static float Gravitation = 100.0f; // ? 27 static float Gravitation = 100.0f; // ?
27 static float delta = 0.05f; // 0.01 ~ 100 ? 28 static float delta = 0.05f; // 0.01 ~ 100 ?
28 static float FIELD = 400.0f; // -100 ~ 100
29 static const float eps = 15.0f; 29 static const float eps = 15.0f;
30 /*
31 * だいたい、平均距離をn倍にするんだったら、Gravitationoはn^2倍、
32 * epsはn倍
33 */
30 34
31 /* for OpenGL Utility. */ 35 /* for OpenGL Utility. */
32 GLUquadricObj **sphere; 36 GLUquadricObj **sphere;
33 float DefDistance = 300.0f; 37 float DefDistance = 300.0f;
34 38
50 __code randomInit(SDL_Surface *screen, int num); 54 __code randomInit(SDL_Surface *screen, int num);
51 __code starsInit(SDL_Surface *screen, int num); 55 __code starsInit(SDL_Surface *screen, int num);
52 __code loop(int count, SDL_Surface *screen, int num); 56 __code loop(int count, SDL_Surface *screen, int num);
53 __code compute(int count, SDL_Surface *screen, int num); 57 __code compute(int count, SDL_Surface *screen, int num);
54 __code nextTurn(int count, SDL_Surface *screen, int num); 58 __code nextTurn(int count, SDL_Surface *screen, int num);
55 __code moveCenter(int count, SDL_Surface *screen, int num); 59 __code drawStars(SDL_Surface *screen, int num);
56 __code CenteringVelocity(int count, SDL_Surface *screen, int num);
57 __code drawStars(int count, SDL_Surface *screen, int num);
58 __code AdjustLooking(float distance, int count, SDL_Surface *screen, int num); 60 __code AdjustLooking(float distance, int count, SDL_Surface *screen, int num);
59 61
60 int main(int argc, char **argv) 62 int main(int argc, char **argv)
61 { 63 {
62 int ch; 64 int ch;
135 gluPerspective( 60.0, (float)screen->w/(float)screen->h, 1.0, 1000.0); 137 gluPerspective( 60.0, (float)screen->w/(float)screen->h, 1.0, 1000.0);
136 gluLookAt( 300.0,300.0,300.0, 0.0,0.0,0.0, 1.0,0.0,0.0); 138 gluLookAt( 300.0,300.0,300.0, 0.0,0.0,0.0, 1.0,0.0,0.0);
137 glClearColor(0.0, 0.0, 0.0, 0.0); 139 glClearColor(0.0, 0.0, 0.0, 0.0);
138 glMatrixMode(GL_MODELVIEW); 140 glMatrixMode(GL_MODELVIEW);
139 141
140 //DEBUGlog("scr->w=%d\n", screen->w);
141 //DEBUGlog("scr->h=%d\n", screen->h);
142 for( i=0; i<num; i++){ 142 for( i=0; i<num; i++){
143 sphere[i] = gluNewQuadric(); 143 sphere[i] = gluNewQuadric();
144 gluQuadricDrawStyle(sphere[i], GLU_FILL); 144 gluQuadricDrawStyle(sphere[i], GLU_FILL);
145 } 145 }
146 146
197 } 197 }
198 198
199 goto loop(0, screen, num); 199 goto loop(0, screen, num);
200 } 200 }
201 201
202
203 /* MAIN LOOP ROUTIN
204 *
205 +------+
206 | loop |<-----------------------+
207 +---+--+ |
208 +-----------+ |
209 |count=num | count<num |
210 +---v------+ +--v-------------+ |
211 | nextTurn | | COMPUTE ROUTIN | |
212 +---+------+ +-------+--------+ |
213 | | |
214 +---v---------+ v |
215 | DRAW ROUTIN +----->+----------+
216 +-------------+ back to loop
217 */
202 __code loop(int count, SDL_Surface *screen, int num) 218 __code loop(int count, SDL_Surface *screen, int num)
203 { 219 {
204 SDL_Event event; 220 SDL_Event event;
205 221
206 usleep(80); 222 usleep(1000);
207 /* check SDL event. */ 223 /* check SDL event. */
208 while(SDL_PollEvent(&event)){ //Poll events 224 while(SDL_PollEvent(&event)){ //Poll events
209 switch(event.type){ //Check event type 225 switch(event.type){ //Check event type
210 case SDL_QUIT: //User hit the X (or equivelent) 226 case SDL_QUIT: //User hit the X (or equivelent)
211 goto finish(1, num); 227 goto finish(1, num);
212 //case SDL_VIDEORESIZE: //User resized window
213 //screen = SDL_SetVideoMode(event.resize.w, event.resize.h, 32,
214 //SDL_HWSURFACE | SDL_RESIZABLE); // Create new window
215 //break; //Event handled, fetch next :)
216 case SDL_KEYDOWN: 228 case SDL_KEYDOWN:
217 if (event.key.state==SDL_PRESSED){ 229 if (event.key.state==SDL_PRESSED){
218 switch(event.key.keysym.sym){ 230 switch(event.key.keysym.sym){
219 case SDLK_UP: 231 case SDLK_UP:
220 goto AdjustLooking(DefDistance/=1.2, count, screen, num); 232 goto AdjustLooking(DefDistance/=1.2, count, screen, num);
253 glMatrixMode(GL_MODELVIEW); 265 glMatrixMode(GL_MODELVIEW);
254 266
255 goto loop(0, screen, num); 267 goto loop(0, screen, num);
256 } 268 }
257 269
258 __code CenteringVelocity(int count, SDL_Surface *screen, int num)
259 {
260 int i;
261 float v0,v1,v2;
262 v0=v1=v2=0.0;
263
264 for (i=0; i<num; i++){
265 v0 += stars_old[i].v[0];
266 v1 += stars_old[i].v[1];
267 v2 += stars_old[i].v[2];
268 }
269 v0/=(float)num; v1/=(float)num; v2/=(float)num;
270 for (i=0; i<num; i++){
271 stars_old[i].v[0] -= v0;
272 stars_old[i].v[1] -= v1;
273 stars_old[i].v[2] -= v2;
274 }
275
276 goto loop(0, screen, num);
277 }
278 __code moveCenter(int count, SDL_Surface *screen, int num)
279 {
280 int i;
281 float m0,m1,m2;
282 m0=m1=m2=0.0;
283
284 for (i=0; i<num; i++){
285 m0 += stars_old[i].r[0];
286 m1 += stars_old[i].r[1];
287 m2 += stars_old[i].r[2];
288 }
289 m0/=(float)num; m1/=(float)num; m2/=(float)num;
290 for (i=0; i<num; i++){
291 stars_old[i].r[0] -= m0;
292 stars_old[i].r[1] -= m1;
293 stars_old[i].r[2] -= m2;
294 }
295
296 goto loop(0, screen, num);
297 }
298
299 /* 270 /*
271 * COMPUTING ROUTINs
272 *
300 * next x = x+dx 273 * next x = x+dx
301 * dx = v*dt => dx/dt = v 274 * dx = v*dt => dx/dt = v
302 * dv = a*dt => dv/dt = a 275 * dv = a*dt => dv/dt = a
303 * a = F/m; 276 * a = F/m;
304 * F = F1 + F2 + ... + Fn 277 * F = F1 + F2 + ... + Fn
305 * Fi = G m1m2/r^2 278 * Fi = G m1m2/r^2
306 * 279 |
280 |
281 +-----v-----+
282 | compute |
283 +-----+-----+
284 |<----------+
285 +-----v---------+ | i++
286 | computeForce +-+
287 +-----+---------+ i<num
288 | i==num
289 +-----+---------+
290 | com|ute_out |
291 +-----+---------+
292 v to loop
307 */ 293 */
294 __code compute_out(int count, SDL_Surface *screen, int num);
295 __code computeForce(int i, int count, SDL_Surface *screen, int num);
308 __code compute(int count, SDL_Surface *screen, int num) 296 __code compute(int count, SDL_Surface *screen, int num)
309 { 297 {
310 int i;
311 /* a is accel this planet receive now. */ 298 /* a is accel this planet receive now. */
312 stars_old[count].a[0]=stars_old[count].a[1]=stars_old[count].a[2]=0.0f; 299 stars_old[count].a[0]=stars_old[count].a[1]=stars_old[count].a[2]=0.0f;
313 DEBUGlog("count=%d\n", count); 300 DEBUGlog("count=%d\n", count);
314 for (i=0; i<num; i++){ // this loop should be split into few code segment.. 301
315 //float F; 302 goto computeForce( 0, count, screen, num);
316 float d0, d1, d2, d; 303 }
317 float a; 304 __code computeForce(int i, int count, SDL_Surface *screen, int num)
318 //body *o, *m; 305 {
319 //o = &stars_old[i]; m = &stars_old[count]; 306 //float F;
320 307 float d0, d1, d2, d;
321 /* skip compute with itself. */ 308 float a;
322 if ( i==count ) continue; 309 //body *o, *m;
323 /* compute distance between two i-th planet and itself. */ 310 //o = &stars_old[i]; m = &stars_old[count];
324 d0 = stars_old[i].r[0] - stars_old[count].r[0]; 311
325 d1 = stars_old[i].r[1] - stars_old[count].r[1]; 312 if (i>=num) goto compute_out(count, screen, num);
326 d2 = stars_old[i].r[2] - stars_old[count].r[2]; 313 /* skip compute with itself. */
327 d = ( d0*d0+d1*d1+d2*d2+eps*eps ); 314 else if (i==count) goto computeForce(i+1, count, screen, num);
328 /* compute force it receive from i-th planet. */ 315
329 //F = Gravitation * stars_old[i].weight * stars_old[count].weight / d; 316 /* compute distance between two i-th planet and itself. */
330 /* and accel. */ 317 d0 = stars_old[i].r[0] - stars_old[count].r[0];
331 //a = F/stars_old[count].weight; 318 d1 = stars_old[i].r[1] - stars_old[count].r[1];
332 a = Gravitation/d * stars_old[i].weight ; 319 d2 = stars_old[i].r[2] - stars_old[count].r[2];
333 stars_old[count].a[0] += a*d0/sqrt(d); 320 d = ( d0*d0+d1*d1+d2*d2+eps*eps );
334 stars_old[count].a[1] += a*d1/sqrt(d); 321 /* compute force it receive from i-th planet. */
335 stars_old[count].a[2] += a*d2/sqrt(d); 322 //F = Gravitation * stars_old[i].weight * stars_old[count].weight / d;
336 DEBUGlog("a=%e, d=%e, d0=%e\n", a, d, d0); 323 /* and accel. */
337 DEBUGlog("g*w=%e\n", Gravitation*stars_old[i].weight); 324 //a = F/stars_old[count].weight;
338 } 325 a = Gravitation/d * stars_old[i].weight ;
326 stars_old[count].a[0] += a*d0/sqrtf(d);
327 stars_old[count].a[1] += a*d1/sqrtf(d);
328 stars_old[count].a[2] += a*d2/sqrtf(d);
329 DEBUGlog("a=%e, d=%e, d0=%e\n", a, d, d0);
330 DEBUGlog("g*w=%e\n", Gravitation*stars_old[i].weight);
331
332 goto computeForce( i+1, count, screen, num);
333 }
334 __code compute_out(int count, SDL_Surface *screen, int num)
335 {
339 336
340 stars_new[count].v[0] = stars_old[count].v[0]+stars_old[count].a[0]*delta; 337 stars_new[count].v[0] = stars_old[count].v[0]+stars_old[count].a[0]*delta;
341 stars_new[count].v[1] = stars_old[count].v[1]+stars_old[count].a[1]*delta; 338 stars_new[count].v[1] = stars_old[count].v[1]+stars_old[count].a[1]*delta;
342 stars_new[count].v[2] = stars_old[count].v[2]+stars_old[count].a[2]*delta; 339 stars_new[count].v[2] = stars_old[count].v[2]+stars_old[count].a[2]*delta;
343 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*/; 340 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*/;
350 DEBUGlog("vx = %e\n", stars_new[count].v[0]); 347 DEBUGlog("vx = %e\n", stars_new[count].v[0]);
351 DEBUGlog("vy = %e\n", stars_new[count].v[1]); 348 DEBUGlog("vy = %e\n", stars_new[count].v[1]);
352 DEBUGlog("vz = %e\n", stars_new[count].v[2]); 349 DEBUGlog("vz = %e\n", stars_new[count].v[2]);
353 DEBUGlog("a0=%e,a1=%e,a2=%e\n", stars_new[count].a[0], stars_new[count].a[1], stars_new[count].a[2]); 350 DEBUGlog("a0=%e,a1=%e,a2=%e\n", stars_new[count].a[0], stars_new[count].a[1], stars_new[count].a[2]);
354 351
355 goto loop(++count, screen, num); 352 goto loop(1+count, screen, num);
356 } 353 }
357 354
358 /* it is executed after all bodies parameter is computed. */ 355 /* it is executed after all bodies parameter is computed. */
359 __code nextTurn(int count, SDL_Surface *screen, int num) 356 __code nextTurn(int count, SDL_Surface *screen, int num)
360 { 357 {
361 body *tmp; 358 body *tmp;
362 tmp = stars_old; 359 tmp = stars_old;
363 stars_old = stars_new; 360 stars_old = stars_new;
364 stars_new = tmp; 361 stars_new = tmp;
365 SDL_UpdateRect(screen, 0, 0, 0, 0); 362 //SDL_UpdateRect(screen, 0, 0, 0, 0);
366 363
367 goto drawStars(count, screen, num); 364 goto drawStars(screen, num);
368 } 365 }
369 366
370 __code drawStars(int count, SDL_Surface *screen, int num) 367
371 { 368
372 int i; 369 /*
370 * DRAWING ROUTINs
371 *
372 |
373 +----+------+
374 | drawStars |
375 +----+------+
376 |<--------+
377 +----v-----+ | count++
378 | drawStar +---+
379 +----+-----+ count<num
380 | count==num
381 +----v----------+
382 | drawStars_out |
383 +----+----------+
384 | to loop
385 */
386 __code drawStar(int i, SDL_Surface *screen, int num);
387 __code drawStars_out(SDL_Surface *screen, int num);
388
389 __code drawStars(SDL_Surface *screen, int num)
390 {
373 DEBUGlog("draw\n"); 391 DEBUGlog("draw\n");
374 glClear( GL_COLOR_BUFFER_BIT ); 392 glClear( GL_COLOR_BUFFER_BIT );
375 for (i=0; i<num; i++){
376 glPushMatrix();
377 glTranslatef( stars_new[i].r[0], stars_new[i].r[1], stars_new[i].r[2]);
378 gluSphere( sphere[i], 2.0, 32.0, 32.0 );
379 glPopMatrix();
380 }
381
382 glColor3d( 1.0, 1.0, 1.0); 393 glColor3d( 1.0, 1.0, 1.0);
383 glBegin(GL_LINES); 394 glBegin(GL_LINES);
384 glVertex3d( 0.0, 0.0, 0.0); 395 glVertex3d( 0.0, 0.0, 0.0);
385 glVertex3d( 200.0, 0.0, 0.0); 396 glVertex3d( 200.0, 0.0, 0.0);
386 glEnd(); 397 glEnd();
391 glBegin(GL_LINES); 402 glBegin(GL_LINES);
392 glVertex3d( 0.0, 0.0, 0.0); 403 glVertex3d( 0.0, 0.0, 0.0);
393 glVertex3d( 0.0, 0.0, 200.0); 404 glVertex3d( 0.0, 0.0, 200.0);
394 glEnd(); 405 glEnd();
395 406
407 goto drawStar(0, screen, num);
408 }
409
410 __code drawStars_out(SDL_Surface *screen, int num)
411 {
412 DEBUGlog("draw axis: screen=%p, num=%d\n", screen, num);
396 SDL_GL_SwapBuffers(); 413 SDL_GL_SwapBuffers();
414 DEBUGlog("go loop\n");
397 goto loop(0, screen, num); 415 goto loop(0, screen, num);
398 } 416 }
399 417
418 __code drawStar(int i, SDL_Surface *screen, int num)
419 {
420 if (i==num)
421 goto drawStars_out(screen, num);
422
423 DEBUGlog("i=%d, num=%d\n", i, num);
424 glPushMatrix();
425 glTranslatef( stars_new[i].r[0], stars_new[i].r[1], stars_new[i].r[2]);
426 gluSphere( sphere[i], 2.0, 32.0, 32.0 );
427 glPopMatrix();
428
429 DEBUGlog("to myself\n");
430 goto drawStar(i+1, screen, num);
431 }
432