10 #if defined(SORT_STAR) || defined (SORT_ENERGYSTAR)
24 safeSumSortedArray(double *arr, int N) {
29 arr[i>>1] = arr[i]+arr[i+1];
39 initGalaxy(Galaxy *galaxy, u_int32_t nstars) {
40 memset(galaxy,0,sizeof(Galaxy));
44 galaxy->nstars = nstars;
46 galaxy->stars=(Star*)t0malloc( nstars*sizeof(Star) );
47 galaxy->forces=(Vector*)t0malloc( CNT_FARRAY(nstars)*sizeof(Vector) );
51 freeGalaxy(Galaxy *galaxy) {
55 tfree(galaxy->forces);
56 memset(galaxy,0,sizeof(Galaxy));
60 cntForce(Vector *force, Star *istar, Star *jstar) {
62 R2 = (jstar->c.x - istar->c.x)*(jstar->c.x - istar->c.x) +
63 (jstar->c.y - istar->c.y)*(jstar->c.y - istar->c.y) +
64 (jstar->c.z - istar->c.z)*(jstar->c.z - istar->c.z);
66 F = G*istar->mass*jstar->mass / (R2*sqrt(R2));
68 force->x = F * fabs(jstar->c.x - istar->c.x);
69 force->y = F * fabs(jstar->c.y - istar->c.y);
70 force->z = F * fabs(jstar->c.z - istar->c.z);
74 forceGalaxy(Galaxy *galaxy) {
75 Vector *force = galaxy->forces;
76 Star *jstar, *istar = galaxy->stars;
78 while( istar < galaxy->stars + galaxy->nstars - 1 ) {
80 while( jstar < galaxy->stars + galaxy->nstars ) {
81 cntForce(force, istar, jstar);
90 #define SUMFORCES(AXIS) \
93 cmpAcc_##AXIS(const void *a, const void *b) { \
94 return ( fabs((*(Vector**)a)->AXIS) > fabs((*(Vector**)b)->AXIS) ) ? 1 : -1; \
98 sumAxis_##AXIS( Vector **a, int n ) { \
99 int N = maxN2(n), i; \
102 qsort(a, n, sizeof(Vector*), cmpAcc_##AXIS); \
104 memset(arr, 0, sizeof(double)*(N-n)); \
106 arr[i+N-n] = a[i]->AXIS; \
108 return safeSumSortedArray(arr, N); \
119 accelerationGalaxy(Galaxy *galaxy) {
121 Star *istar = galaxy->stars,*jstar;
123 Vector* ptrforces[galaxy->nstars];
128 for(i=0; i<galaxy->nstars; i++ ) {
129 istar = galaxy->stars+i;
130 jstar = galaxy->stars;
132 memcpy( &(istar->pa), &(istar->a), sizeof(Vector) );
135 for(j=0; j<galaxy->nstars; j++ ) {
137 ptrforces[j] = GET_FORCE(galaxy, j, i);
138 ptrforces[j]->x = ( jstar->c.x>istar->c.x ) ? fabs(ptrforces[j]->x) : -fabs(ptrforces[j]->x);
139 ptrforces[j]->y = ( jstar->c.y>istar->c.y ) ? fabs(ptrforces[j]->y) : -fabs(ptrforces[j]->y);
140 ptrforces[j]->z = ( jstar->c.z>istar->c.z ) ? fabs(ptrforces[j]->z) : -fabs(ptrforces[j]->z);
142 ptrforces[j-1] = GET_FORCE(galaxy, i, j);
143 cntForce( ptrforces[j-1], istar, jstar );
144 ptrforces[j-1]->x = ( jstar->c.x>istar->c.x ) ? ptrforces[j-1]->x : -ptrforces[j-1]->x;
145 ptrforces[j-1]->y = ( jstar->c.y>istar->c.y ) ? ptrforces[j-1]->y : -ptrforces[j-1]->y;
146 ptrforces[j-1]->z = ( jstar->c.z>istar->c.z ) ? ptrforces[j-1]->z : -ptrforces[j-1]->z;
151 istar->a.x = sumAxis_x(ptrforces, galaxy->nstars-1);
152 istar->a.y = sumAxis_y(ptrforces, galaxy->nstars-1);
153 istar->a.z = sumAxis_z(ptrforces, galaxy->nstars-1);
156 memset(&(istar->a), 0, sizeof(Vector));
158 force=GET_FORCE(galaxy, 0, i);
160 istar->a.x += ( jstar->c.x>istar->c.x ) ? force->x : -force->x;
161 istar->a.y += ( jstar->c.y>istar->c.y ) ? force->y : -force->y;
162 istar->a.z += ( jstar->c.z>istar->c.z ) ? force->z : -force->z;
166 for(j=i+1;j<galaxy->nstars;j++) {
167 force = GET_FORCE(galaxy, i, j);
168 cntForce(force, istar, jstar);
170 istar->a.x += ( jstar->c.x>istar->c.x ) ? force->x : -force->x;
171 istar->a.y += ( jstar->c.y>istar->c.y ) ? force->y : -force->y;
172 istar->a.z += ( jstar->c.z>istar->c.z ) ? force->z : -force->z;
177 istar->a.x /= istar->mass;
178 istar->a.y /= istar->mass;
179 istar->a.z /= istar->mass;
182 if ( galaxy->firsttime )
183 memcpy( &(istar->pa), &(istar->a), sizeof(Vector) );
187 galaxy->firsttime = 0;
192 stepGalaxy(Galaxy *galaxy, double delta) {
193 Star *star = galaxy->stars;
195 while( star - galaxy->stars < galaxy->nstars ) {
197 star->c.x += delta*( delta*( (star->a.x - star->pa.x)/3.0 + star->a.x )/2.0 + star->v.x );
198 star->v.x += delta*( star->a.x + (star->a.x - star->pa.x)/2.0 );
199 star->c.y += delta*( delta*( (star->a.y - star->pa.y)/3.0 + star->a.y )/2.0 + star->v.y );
200 star->v.y += delta*( star->a.y + (star->a.y - star->pa.y)/2.0 );
201 star->c.z += delta*( delta*( (star->a.z - star->pa.z)/3.0 + star->a.z )/2.0 + star->v.z );
202 star->v.z += delta*( star->a.z + (star->a.z - star->pa.z)/2.0 );
204 star->c.x += delta*(star->v.x + star->a.x*delta/2.0);
205 star->v.x += star->a.x * delta;
206 star->c.y += delta*(star->v.y + star->a.y*delta/2.0);
207 star->v.y += star->a.y * delta;
208 star->c.z += delta*(star->v.z + star->a.z*delta/2.0);
209 star->v.z += star->a.z * delta;
217 unstepGalaxy(Galaxy *galaxy, double delta) {
218 Star *star = galaxy->stars;
220 while( star - galaxy->stars < galaxy->nstars ) {
222 star->v.x -= delta*( star->a.x + (star->a.x - star->pa.x)/2.0 );
223 star->c.x -= delta*( delta*( (star->a.x - star->pa.x)/3.0 + star->a.x )/2.0 + star->v.x );
224 star->v.y -= delta*( star->a.y + (star->a.y - star->pa.y)/2.0 );
225 star->c.y -= delta*( delta*( (star->a.y - star->pa.y)/3.0 + star->a.y )/2.0 + star->v.y );
226 star->v.z -= delta*( star->a.z + (star->a.z - star->pa.z)/2.0 );
227 star->c.z -= delta*( delta*( (star->a.z - star->pa.z)/3.0 + star->a.z )/2.0 + star->v.z );
229 star->v.x -= star->a.x * delta;
230 star->c.x -= delta*(star->v.x + star->a.x*delta/2.0);
231 star->v.y -= star->a.y * delta;
232 star->c.y -= delta*(star->v.y + star->a.y*delta/2.0);
233 star->v.z -= star->a.z * delta;
234 star->c.z -= delta*(star->v.z + star->a.z*delta/2.0);
241 #ifdef SORT_ENERGYSTAR
243 getImpulseX(Star *star) {
244 return star->mass * star->v.x;
248 getImpulseY(Star *star) {
249 return star->mass * star->v.y;
253 getImpulseZ(Star *star) {
254 return star->mass * star->v.z;
258 getMomentX(Star *star) {
259 return star->mass * ( star->v.y * star->c.z - star->c.y * star->v.z );
263 getMomentY(Star *star) {
264 return star->mass * ( star->v.z * star->c.x - star->c.z * star->v.x );
268 getMomentZ(Star *star) {
269 return star->mass * ( star->v.x * star->c.y - star->c.x * star->v.y );
273 getKinetic(Star *star) {
274 return star->mass * (star->v.x * star->v.x + star->v.y * star->v.y + star->v.z * star->v.z) / 2.0;
278 cmpDouble(const void *a, const void *b) {
279 return ( fabs(*(double*)a) > fabs(*(double*)b) ) ? 1 : -1;
283 safeSumArray(double *arr, int N) {
284 qsort(arr, N, sizeof(double), cmpDouble);
285 return safeSumSortedArray(arr,N);
289 safeSumStar(Star *stars, int n, double (*extractor)(Star*)) {
294 arr[i] = (*extractor)(stars+i);
295 memset( arr+n, 0, sizeof(double)*(N-n) );
297 return safeSumArray(arr, N);
303 EnergyImpulseGalaxy(Galaxy *galaxy) {
304 Star *istar=galaxy->stars, *jstar;
305 double potential=0.0;
307 double impulseX=0.0, impulseY=0.0, impulseZ=0.0, momentX=0.0, momentY=0.0, momentZ=0.0;
308 #ifdef SORT_ENERGYSTAR
309 double pe[ maxN2(CNT_FARRAY(galaxy->nstars)) ], *ptrp;
312 memset(pe, 0, sizeof(double)*maxN2(CNT_FARRAY(galaxy->nstars)));
315 while( istar - galaxy->stars < galaxy->nstars ) {
317 #ifndef SORT_ENERGYSTAR
318 impulseX = impulseX + istar->mass * istar->v.x;
319 impulseY = impulseY + istar->mass * istar->v.y;
320 impulseZ = impulseZ + istar->mass * istar->v.z;
322 momentX += istar->mass * ( istar->v.y * istar->c.z - istar->c.y * istar->v.z );
323 momentY += istar->mass * ( istar->v.z * istar->c.x - istar->c.z * istar->v.x );
324 momentZ += istar->mass * ( istar->v.x * istar->c.y - istar->c.x * istar->v.y );
326 kinetic += istar->mass * (istar->v.x * istar->v.x + istar->v.y * istar->v.y + istar->v.z * istar->v.z) / 2.0;
330 while( jstar - galaxy->stars < galaxy->nstars ) {
331 #ifdef SORT_ENERGYSTAR
336 G * istar->mass * jstar->mass / sqrt(
337 (jstar->c.x - istar->c.x)*(jstar->c.x - istar->c.x) +
338 (jstar->c.y - istar->c.y)*(jstar->c.y - istar->c.y) +
339 (jstar->c.z - istar->c.z)*(jstar->c.z - istar->c.z)
346 #ifdef SORT_ENERGYSTAR
347 potential = safeSumArray(pe, maxN2(ptrp-pe));
348 impulseX = safeSumStar(galaxy->stars, galaxy->nstars, getImpulseX);
349 impulseY = safeSumStar(galaxy->stars, galaxy->nstars, getImpulseY);
350 impulseZ = safeSumStar(galaxy->stars, galaxy->nstars, getImpulseZ);
351 momentX = safeSumStar(galaxy->stars, galaxy->nstars, getMomentX);
352 momentY = safeSumStar(galaxy->stars, galaxy->nstars, getMomentY);
353 momentZ = safeSumStar(galaxy->stars, galaxy->nstars, getMomentZ);
354 kinetic = safeSumStar(galaxy->stars, galaxy->nstars, getKinetic);
356 galaxy->Impulse = sqrt(impulseX*impulseX + impulseY*impulseY + impulseZ*impulseZ);
357 galaxy->Moment = sqrt(momentX *momentX + momentY *momentY + momentZ *momentZ );
358 galaxy->kineticEnergy = kinetic;
359 if ( !finite(potential) ) potential=0;
360 galaxy->potentialEnergy = -potential;
364 initLiveGalaxy(Galaxy *galaxy) {
365 if ( galaxy->errorLimit<=0) galaxy->errorLimit = 1e-8;
366 if ( galaxy->desiredDelta<=0 ) galaxy->desiredDelta=3600;
368 galaxy->error = galaxy->errorLimit;
369 galaxy->delta=galaxy->desiredDelta;
370 galaxy->elapsedTime=0.0;
375 EnergyImpulseGalaxy(galaxy);
379 #define iszero(x) ( fpclassify(x) & FP_ZERO )
380 #define CHK_ERRORVAL( what, err, val, prevval ) \
381 if ( finite(val) && finite(prevval) ) { \
382 err = fabs( (val)/(prevval) - 1.0 ); \
383 if ( !finite(err) ) \
384 err = galaxy->error; \
386 fprintf(stderr, "%s is not finite: val:%G prevval:%G\n", what, (val), (prevval)); \
387 err = galaxy->error; \
391 liveGalaxy(Galaxy *galaxy, pthread_mutex_t *mutex) {
392 double prevEnergy = galaxy->kineticEnergy + galaxy->potentialEnergy;
393 double prevImpulse = galaxy->Impulse;
394 double prevMoment = galaxy->Moment;
395 double errorI, errorE, errorM, error, prevError;
397 accelerationGalaxy(galaxy);
398 prevError = galaxy->error;
400 if ( mutex ) pthread_mutex_lock(mutex);
402 stepGalaxy(galaxy, galaxy->delta);
403 EnergyImpulseGalaxy(galaxy);
405 CHK_ERRORVAL("Energy", errorE, galaxy->kineticEnergy + galaxy->potentialEnergy, prevEnergy);
406 CHK_ERRORVAL("Impulse",errorI, galaxy->Impulse, prevImpulse);
407 CHK_ERRORVAL("Moment", errorM, galaxy->Moment, prevMoment);
409 /*error=fmax(errorE, fmax(errorI, errorM));*/
412 if ( error > galaxy->errorLimit ) {
413 unstepGalaxy(galaxy, galaxy->delta);
414 galaxy->delta /= 2.0;
415 } else if ( error*1e2 < galaxy->errorLimit && prevError < galaxy->errorLimit ) {
416 unstepGalaxy(galaxy, galaxy->delta);
417 galaxy->delta *= 2.0;
419 galaxy->error = error;
424 if ( mutex ) pthread_mutex_unlock(mutex);
426 galaxy->elapsedTime += galaxy->delta;
430 printStar(Star *star) {
431 printf("\tMass: %e\n", star->mass);
432 printf("\tPosit: x:%e y:%e z:%e\n", star->c.x, star->c.y, star->c.z);
433 printf("\tSpeed: x:%e y:%e z:%e\n", star->v.x, star->v.y, star->v.z);