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