Mercurial > hg > Members > kent > N-BodyProblem

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