galaxy.c

   1 #include <math.h>
   2 #include <string.h>
   3
   4
   5 #include "tmalloc.h"
   6 #include "tlog.h"
   7
   8 #include "galaxy.h"
   9
  10 #if defined(SORT_STAR) || defined (SORT_ENERGYSTAR)
  11 static int
  12 maxN2(int n) {
  13         int i=0;
  14
  15         while(n) {
  16                 n >>= 1;
  17                 i++;
  18         }
  19
  20         return 1<<i;
  21 }
  22
  23 static double
  24 safeSumSortedArray(double *arr, int N) {
  25         int i;
  26
  27         while(N>1) {
  28                 for(i=0;i<N;i+=2)
  29                         arr[i>>1] = arr[i]+arr[i+1];
  30                 N>>=1;
  31         }
  32
  33         return arr[0];
  34 }
  35
  36 #endif
  37
  38 void
  39 initGalaxy(Galaxy *galaxy, u_int32_t nstars) {
  40         memset(galaxy,0,sizeof(Galaxy));
  41
  42         if ( nstars < 2 )
  43                 nstars=2;
  44         galaxy->nstars = nstars;
  45
  46         galaxy->stars=(Star*)t0malloc( nstars*sizeof(Star) );
  47         galaxy->forces=(Vector*)t0malloc( CNT_FARRAY(nstars)*sizeof(Vector) );
  48 }
  49
  50 void
  51 freeGalaxy(Galaxy *galaxy) {
  52         if ( galaxy->stars )
  53                 tfree(galaxy->stars);
  54         if ( galaxy->forces )
  55                 tfree(galaxy->forces);
  56         memset(galaxy,0,sizeof(Galaxy));
  57 }
  58
  59 static void
  60 cntForce(Vector *force, Star *istar, Star *jstar) {
  61         double          R2,F;
  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);
  65
  66         F = G*istar->mass*jstar->mass / (R2*sqrt(R2));
  67
  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);
  71 }
  72
  73 void
  74 forceGalaxy(Galaxy *galaxy) {
  75         Vector          *force = galaxy->forces;
  76         Star            *jstar, *istar = galaxy->stars;
  77
  78         while( istar < galaxy->stars + galaxy->nstars - 1 ) {
  79                 jstar = istar+1;
  80                 while( jstar < galaxy->stars + galaxy->nstars ) {
  81                         cntForce(force, istar, jstar);
  82                         force++; jstar++;
  83                 }
  84                 istar++;
  85         }
  86 }
  87
  88 #ifdef  SORT_STAR
  89
  90 #define SUMFORCES(AXIS)                                                                 \
  91                                                                                         \
  92 static int                                                                              \
  93 cmpAcc_##AXIS(const void *a, const void *b) {                                           \
  94         return ( fabs((*(Vector**)a)->AXIS) > fabs((*(Vector**)b)->AXIS) ) ? 1 : -1;    \
  95 }                                                                                       \
  96                                                                                         \
  97 static double                                                                           \
  98 sumAxis_##AXIS( Vector **a, int n ) {                                                   \
  99         int N = maxN2(n), i;                                                            \
 100         double  arr[N];                                                                 \
 101                                                                                         \
 102         qsort(a, n, sizeof(Vector*), cmpAcc_##AXIS);                                    \
 103                                                                                         \
 104         memset(arr, 0, sizeof(double)*(N-n));                                           \
 105         for(i=0;i<n;i++)                                                                \
 106                 arr[i+N-n] = a[i]->AXIS;                                                \
 107                                                                                         \
 108         return safeSumSortedArray(arr, N);                                              \
 109 }
 110
 111 SUMFORCES(x)
 112 SUMFORCES(y)
 113 SUMFORCES(z)
 114
 115 #endif
 116
 117
 118 void
 119 accelerationGalaxy(Galaxy *galaxy) {
 120         u_int32_t       i,j;
 121         Star            *istar = galaxy->stars,*jstar;
 122 #ifdef  SORT_STAR
 123         Vector*         ptrforces[galaxy->nstars];
 124 #else
 125         Vector          *force;
 126 #endif
 127
 128         for(i=0; i<galaxy->nstars; i++ ) {
 129                 istar = galaxy->stars+i;
 130                 jstar = galaxy->stars;
 131
 132                 memcpy( &(istar->pa), &(istar->a), sizeof(Vector) );
 133
 134 #ifdef  SORT_STAR
 135                 for(j=0; j<galaxy->nstars; j++ ) {
 136                         if ( j<i ) {
 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);
 141                         } else if ( j>i ) {
 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;
 147                         }
 148                         jstar++;
 149                 }
 150
 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);
 154
 155 #else
 156                 memset(&(istar->a), 0, sizeof(Vector));
 157
 158                 force=GET_FORCE(galaxy, 0, i);
 159                 for(j=0;j<i;j++) {
 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;
 163                         force++; jstar++;
 164                 }
 165                 jstar++;
 166                 for(j=i+1;j<galaxy->nstars;j++) {
 167                         force = GET_FORCE(galaxy, i, j);
 168                         cntForce(force, istar, jstar);
 169
 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;
 173                         jstar++;
 174                 }
 175 #endif
 176
 177                 istar->a.x /= istar->mass;
 178                 istar->a.y /= istar->mass;
 179                 istar->a.z /= istar->mass;
 180
 181 #ifdef TAYLOR3
 182                 if ( galaxy->firsttime )
 183                         memcpy( &(istar->pa), &(istar->a), sizeof(Vector) );
 184 #endif
 185         }
 186 #ifdef TAYLOR3
 187         galaxy->firsttime = 0;
 188 #endif
 189 }
 190
 191 void
 192 stepGalaxy(Galaxy *galaxy, double delta) {
 193         Star            *star = galaxy->stars;
 194
 195         while( star - galaxy->stars < galaxy->nstars ) {
 196 #ifdef TAYLOR3
 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 );
 203 #else
 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;
 210 #endif
 211
 212                 star++;
 213         }
 214 }
 215
 216 void
 217 unstepGalaxy(Galaxy *galaxy, double delta) {
 218         Star            *star = galaxy->stars;
 219
 220         while( star - galaxy->stars < galaxy->nstars ) {
 221 #ifdef TAYLOR3
 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 );
 228 #else
 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);
 235 #endif
 236
 237                 star++;
 238         }
 239 }
 240
 241 #ifdef SORT_ENERGYSTAR
 242 static double
 243 getImpulseX(Star *star) {
 244         return star->mass * star->v.x;
 245 }
 246
 247 static double
 248 getImpulseY(Star *star) {
 249         return star->mass * star->v.y;
 250 }
 251
 252 static double
 253 getImpulseZ(Star *star) {
 254         return star->mass * star->v.z;
 255 }
 256
 257 static double
 258 getMomentX(Star *star) {
 259         return star->mass * ( star->v.y * star->c.z - star->c.y * star->v.z );
 260 }
 261
 262 static double
 263 getMomentY(Star *star) {
 264         return star->mass * ( star->v.z * star->c.x - star->c.z * star->v.x );
 265 }
 266
 267 static double
 268 getMomentZ(Star *star) {
 269         return star->mass * ( star->v.x * star->c.y - star->c.x * star->v.y );
 270 }
 271
 272 static double
 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;
 275 }
 276
 277 static int
 278 cmpDouble(const void *a, const void *b) {
 279         return ( fabs(*(double*)a) > fabs(*(double*)b) ) ? 1 : -1;
 280 }
 281
 282 static double
 283 safeSumArray(double *arr, int N) {
 284         qsort(arr, N, sizeof(double), cmpDouble);
 285         return safeSumSortedArray(arr,N);
 286 }
 287
 288 static double
 289 safeSumStar(Star *stars, int n, double (*extractor)(Star*)) {
 290         int N=maxN2(n), i;
 291         double arr[N];
 292
 293         for(i=0;i<n;i++)
 294                 arr[i] = (*extractor)(stars+i);
 295         memset( arr+n, 0, sizeof(double)*(N-n) );
 296
 297         return safeSumArray(arr, N);
 298 }
 299
 300 #endif
 301
 302 void
 303 EnergyImpulseGalaxy(Galaxy *galaxy) {
 304         Star            *istar=galaxy->stars, *jstar;
 305         double potential=0.0;
 306         double kinetic=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;
 310
 311         ptrp=pe;
 312         memset(pe, 0, sizeof(double)*maxN2(CNT_FARRAY(galaxy->nstars)));
 313
 314 #endif
 315         while( istar - galaxy->stars < galaxy->nstars ) {
 316
 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;
 321
 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 );
 325
 326                 kinetic += istar->mass * (istar->v.x * istar->v.x + istar->v.y * istar->v.y + istar->v.z * istar->v.z) / 2.0;
 327 #endif
 328
 329                 jstar=istar+1;
 330                 while( jstar - galaxy->stars < galaxy->nstars ) {
 331 #ifdef SORT_ENERGYSTAR
 332                         *ptrp++ =
 333 #else
 334                         potential +=
 335 #endif
 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)
 340                         );
 341                         jstar++;
 342                 }
 343                 istar++;
 344         }
 345
 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);
 355 #endif
 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;
 361 }
 362
 363 void
 364 initLiveGalaxy(Galaxy *galaxy) {
 365         if ( galaxy->errorLimit<=0) galaxy->errorLimit = 1e-8;
 366         if ( galaxy->desiredDelta<=0 ) galaxy->desiredDelta=3600;
 367
 368         galaxy->error = galaxy->errorLimit;
 369         galaxy->delta=galaxy->desiredDelta;
 370         galaxy->elapsedTime=0.0;
 371
 372 #ifdef TAYLOR3
 373         galaxy->firsttime=1;
 374 #endif
 375         EnergyImpulseGalaxy(galaxy);
 376 }
 377
 378
 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;                                                           \
 385 } else {                                                                                        \
 386         fprintf(stderr, "%s is not finite: val:%G prevval:%G\n", what, (val), (prevval));       \
 387         err = galaxy->error;                                                                    \
 388 }
 389
 390 void
 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;
 396
 397         accelerationGalaxy(galaxy);
 398         prevError = galaxy->error;
 399
 400         if ( mutex ) pthread_mutex_lock(mutex);
 401         while(1) {
 402                 stepGalaxy(galaxy, galaxy->delta);
 403                 EnergyImpulseGalaxy(galaxy);
 404
 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);
 408
 409                 /*error=fmax(errorE, fmax(errorI, errorM));*/
 410                 error = errorE;
 411
 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;
 418                 } else {
 419                         galaxy->error = error;
 420                         break;
 421                 }
 422                 prevError=error;
 423         }
 424         if ( mutex ) pthread_mutex_unlock(mutex);
 425
 426         galaxy->elapsedTime += galaxy->delta;
 427 }
 428
 429 void
 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);
 434 }