home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / de118i.zip / FINDCENT.C < prev    next >
C/C++ Source or Header  |  1992-05-24  |  5KB  |  249 lines
  1. #ifndef NOINCS
  2. #include "mconf.h"
  3. #include "prec.h"
  4. #include "ssystem.h"
  5. #include "const.h"
  6. #endif
  7.  
  8. #define FIDEBUG 0
  9.  
  10. extern DOUBLE C;
  11. extern DOUBLE Rij[NTOTAL][NTOTAL];
  12. extern DOUBLE GMs[];
  13. extern DOUBLE EMRAT;
  14. extern DOUBLE yw[];
  15. extern DOUBLE ymoon[];
  16.  
  17.  
  18.  
  19. /* This routine adds a constant vector to move the
  20.  * barycenter to the origin.
  21.  */
  22. extern DOUBLE mustar[NTOTAL];
  23.  
  24. findcenter( t, yin )
  25. DOUBLE t;
  26. DOUBLE yin[];
  27. {
  28. int i, j, k, ii, jj;
  29. DOUBLE xx, yy, zz, vx, vy, vz, csqi, s, mu;
  30.  
  31. #if MOON
  32. fromemb( yin, yw );
  33. #endif
  34. #if AROIDS
  35. aroids( t, yw );
  36. #endif
  37. csqi = Half / (C*C);
  38. for( k=0; k<3; k++ )
  39.     { /* Iterate to find solution of implicit expressions. */
  40.  
  41.     distances(yw);
  42. #if DOREL
  43. /* Relativistic GM of each body */
  44.     ii = 6*FMASS;
  45.     for( i=FMASS; i<NTOTAL; i++ )
  46.         {
  47.         vx = yw[ii]; /* velocity */
  48.         vy = yw[ii+2];
  49.         vz = yw[ii+4];
  50.         s = vx * vx + vy * vy + vz * vz; /* square of velocity */
  51.         for( j=FMASS; j<NTOTAL; j++ )
  52.             {
  53.             if( j == i )
  54.                 continue;
  55.             s -= GMs[j]*Rij[i][j];
  56.             }
  57.         mustar[i] = GMs[i] * (One + csqi * s);
  58.         ii += 6;
  59.         }
  60. #endif
  61.  
  62. /* Compute relativistic (or Newtonian) center of mass. */
  63.     mu = Zero; /* total mass of the system */
  64.     xx = Zero;
  65.     yy = Zero;
  66.     zz = Zero;
  67.     vx = Zero;
  68.     vy = Zero;
  69.     vz = Zero;
  70.     ii = 6*FMASS;
  71.     for( i=FMASS; i<NTOTAL; i++ )
  72.         {
  73. #if DOREL
  74.         s = mustar[i];
  75. #else
  76.         s = GMs[i];
  77. #endif
  78.         mu += s;
  79.         xx += yw[ii+1] * s; /* position */
  80.         yy += yw[ii+3] * s;
  81.         zz += yw[ii+5] * s;
  82.         vx += yw[ii] * s; /* velocity */
  83.         vy += yw[ii+2] * s;
  84.         vz += yw[ii+4] * s;
  85.         ii += 6;
  86.         }
  87. /* Offset the coordinates of all the objects equally. */
  88.     xx /= mu;
  89.     yy /= mu;
  90.     zz /= mu;
  91.     vx /= mu;
  92.     vy /= mu;
  93.     vz /= mu;
  94.     jj = 6*FMASS;
  95.     for( i=FMASS; i<NTOTAL; i++ )
  96.         {
  97.         yw[jj+1] -= xx;
  98.         yw[jj+3] -= yy;
  99.         yw[jj+5] -= zz;
  100.         yw[jj] -= vx;
  101.         yw[jj+2] -= vy;
  102.         yw[jj+4] -= vz;
  103.         jj += 6;
  104.         }
  105. #if FIDEBUG
  106.     if( k == 0 )
  107.         {
  108.         xx = xx*xx + yy*yy + zz*zz;
  109.         vx = vx*vx + vy*vy + vz*vz;
  110.         printf( "Barycenter offset ^2 = %.5Le au, %.5Le au/d\n",
  111.             xx, vx );
  112.         }
  113. #endif
  114.     }
  115.  
  116. #if MOON
  117. jj = 6*FMASS;
  118. for( i=FMASS; i<NTOTAL; i++ )
  119.     {
  120.     if( i == IMOON )
  121.         jj += 6;
  122.     else
  123.         {
  124.         for( j=0; j<6; j++ )
  125.             {
  126.             yin[jj] = yw[jj];
  127.             jj += 1;
  128.             }
  129.         }
  130.     }
  131. /* Convert barycentric Earth and Moon to output EMB variables. */
  132. ii = 6*IEARTH;
  133. jj = 6*IMOON;
  134. for( i=0; i<6; i++ )
  135.     {
  136.     xx = yw[ii+i]; /* SS barycentric Earth */
  137.     yy = yw[jj+i]; /* SS barycentric Moon */
  138.     yin[ii+i] = (EMRAT * xx + yy)/(EMRAT+One); /* EMB */
  139. /*    yin[jj+i] = yy - xx; */ /* M = Moon - Earth */
  140.     yin[jj+i] = ymoon[i];
  141.     }
  142. #endif
  143.  
  144. #if FIDEBUG
  145. printf( "Sun's state:\n" );
  146. for( i=0; i<6; i++ )
  147.     {
  148.     xx = yw[6*ISUN+i];
  149.     printf( "%27.18Le,\n", xx );
  150.     }
  151. #endif
  152. }
  153.  
  154.  
  155. #if FIDEBUG
  156. /* Initial heliocentric accelerations, from DE118
  157.  * JED 2440400.5
  158.  */
  159. DOUBLE iniaccel[] = {
  160. /* Libs */
  161. -2.33840106202781591e-6L,
  162.  3.76106870289525089e-6L,
  163. -8.36366071885551190e-7L,
  164. /* Mercury */
  165. -1.93862993974248726e-3L,
  166.  5.20550192322063647e-4L,
  167.  4.77930057349665243e-4L,
  168.  
  169. -4.62592548867567062e-4L,
  170.  2.72876752608682010e-4L,
  171.  1.52182594409640823e-4L,
  172.  
  173. -2.92050891691249361e-5L,
  174.  2.61318585392418325e-4L,
  175.  1.13316577018344841e-4L,
  176.  
  177.  1.24770701161065107e-5L,
  178.  1.24988885615999961e-4L,
  179.  5.70375473684709056e-5L,
  180.  
  181.  9.85358573091297212e-6L,
  182.  1.40990437369746227e-6L,
  183.  3.64122840684491628e-7L,
  184.  
  185. -2.95364938480744481e-6L,
  186. -1.67622531774450414e-6L,
  187. -5.65095232380545449e-7L,
  188.  
  189.  8.90424390777303382e-7L,
  190.  4.87032002931296781e-8L,
  191.  8.78848553654263136e-9L,
  192.  
  193.  1.80994923472092847e-7L,
  194.  2.55021868335872720e-7L,
  195.  1.00082692642110886e-7L,
  196.  
  197.  2.88404461789641877e-7L,
  198.  7.37133065237483562e-9L,
  199. -8.27405630423931700e-8L,
  200. /* Moon */
  201.  5.40548681996079620e-5L,
  202.  1.26326595513576372e-4L,
  203.  6.91101026282511274e-5L,
  204. /* Sun */
  205. 0.0e0L,
  206. 0.0e0L,
  207. 0.0e0L,
  208. };
  209.  
  210. chkacc(v)
  211. DOUBLE v[];
  212. {
  213. DOUBLE asunx, asuny, asunz, dx, dy, dz;
  214. double px, py, pz;
  215. int i, ii, jj;
  216.  
  217. /* estimate the initial acceleration of the Sun
  218.  * from the given heliocentric acceleration of Pluto
  219.  */
  220. ii = 6*(IMOON-1);
  221. jj = 3*(IMOON-1);
  222. asunx = v[ii] - iniaccel[jj];
  223. asuny = v[ii+2] - iniaccel[jj+1];
  224. asunz = v[ii+4] - iniaccel[jj+2];
  225.  
  226. printf( "\n" );
  227. ii = 0;
  228. jj = 0;
  229. for( i=0; i<11; i++ )
  230.     {
  231.     dx = v[ii]   - iniaccel[jj];
  232.     dy = v[ii+2] - iniaccel[jj+1];
  233.     dz = v[ii+4] - iniaccel[jj+2];
  234.     if( (i != 0) && (i != IMOON) )
  235.         {
  236.         dx -= asunx;
  237.         dy -= asuny;
  238.         dz -= asunz;
  239.         }
  240.     px = dx;
  241.     py = dy;
  242.     pz = dz;
  243.     printf( "%11.3e %11.3e %11.3e\n", px, py, pz );
  244.     ii += 6;
  245.     jj += 3;
  246.     }
  247. }
  248. #endif
  249.