home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ode / ssystem.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  14.5 KB  |  702 lines
  1. /* Main program to perform a point-mass integration of the
  2.  * Solar system, with relativistic corrections.  The initial
  3.  * conditions are those of the JPL DE200 numerical integration.
  4.  * The final solution, given in the vector yn1[], is taken from
  5.  * the DE200 for 100 days later.  A test run of 400 steps in double
  6.  * precison IEEE arithmetic with step size = 1/4 day and 11th
  7.  * order Adams integration yielded agreement to about 10^-10 au
  8.  * for the inner planets, 10^-12 au for Pluto, and 10^-8 au
  9.  * for the Moon.  The DE200, of course, included Earth-Moon figure
  10.  * effects and five of the asteroids, and its arithmetic was
  11.  * probably higher than double precision.
  12.  *
  13.  * If the Adams order is set to N, the first N+1 steps are
  14.  * performed by the Runge-Kutta subroutine rungek.c.  These
  15.  * take about 7 times longer than the subsequent ones done
  16.  * by the Adams-Bashforth-Moulton routine adams.c.
  17.  *
  18.  * Neither the step size nor the approximation order has been
  19.  * made self-adjusting.  Some experimentation is recommended, to
  20.  * achieve the desired balance of speed and accuracy.  The
  21.  * example above was chosen to give about the same theoretical
  22.  * error per day as specified in the DE102 reference below.
  23.  *
  24.  * -Steve Moshier
  25.  */
  26.  
  27. #include "int.h"
  28.  
  29. /* Compute relativity corrections, or not: */
  30. #define DOREL 1
  31.  
  32. /* Include the Moon as a separate body, or not: */
  33. #define MOON 1
  34.  
  35. /* Constrain the relativistic barycenter of the solar
  36.  * system to stay at the origin.  Note, with the asteroids
  37.  * (10^-9 solar masses) missing, this step is possibly useless
  38.  * and was omitted in the test example.  The code for it is
  39.  * supplied at the end of this file.
  40.  */
  41. #define BARYC 0
  42.  
  43. double sqrt();
  44.  
  45. #if MOON
  46. #define NBODY 11
  47. #define IMOON 9
  48. #else
  49. #define NBODY 10
  50. #endif
  51. #define ISUN (NBODY-1)
  52. #define NTOTAL (NBODY)
  53. #define NEQ (6*NBODY)
  54. #define MAXORD 13
  55.  
  56. #if DOREL
  57. #if BARYC
  58. #define NEQC (NEQ-6)
  59. #else
  60. #define NEQC NEQ
  61. #endif
  62. #else
  63. #define NEQC (NEQ)
  64. #endif
  65.  
  66. static double vnewt[NEQ];
  67. static double awork[ (NEQ*(MAXORD+2))+(MAXORD+2)+(6*NEQ) ];
  68. static double rwork[(MAXRUNG+2)*NEQ];
  69. static double Rij[NTOTAL][NTOTAL];
  70.  
  71. /* Speed of light. Note, all dimensions are in astronomical
  72.  * units (au) and days (86,400 seconds).
  73.  */
  74. #define C 173.1446328
  75.  
  76. /* Solar system barycentric velocity and position
  77.  * state of Mercury, Venus, EMB, ..., Pluto, MOON, Sun
  78.  * in the order dx/dt, x, dy/dt, y, dz/dt, z for each object.
  79.  *
  80.  * EMB is the arithmetic average of Earth and Moon, weighted by
  81.  * their masses. The coordinates of the MOON variable are the difference
  82.  * between the solar system barycentric coordinates of the Moon
  83.  * minus those of the Earth.
  84.  */
  85.  
  86. /* Julian date of the initial state yn[]
  87.  * June 28.0, 1969:
  88.  */
  89. #define JD0 2440400.50
  90.  
  91. static double yn[] = {
  92.  3.367492758813184e-003, /* Mercury */
  93.  3.617650084026384e-001,
  94.   2.489452047536652e-002,
  95.  -9.077580081092690e-002,
  96.   1.294630071041128e-002,
  97.  -8.571250013597834e-002,
  98.  
  99.  1.095206748972985e-002, /* Venus */
  100.  6.127542336221157e-001,
  101.   1.561768454932838e-002,
  102.  -3.483591895730089e-001,
  103.   6.331105550175031e-003,
  104.  -1.952758215816601e-001,
  105.  
  106.  1.681126759384403e-002, /* EMB */
  107.  1.205197129002194e-001,
  108.   1.748308899617618e-003,
  109.  -9.258323035784091e-001,
  110.   7.582041930175932e-004,
  111.  -4.015377611192672e-001,
  112.  
  113.  1.448165239083299e-002, /* Mars */
  114. -1.101837696063662e-001,
  115.   2.424628123494023e-004,
  116.  -1.327593278911103e+000,
  117.  -2.815196436092114e-004,
  118.  -6.058866834707231e-001,
  119.  
  120.  1.092015004926191e-003, /* Jupiter */
  121. -5.379703845550703e+000,
  122.  -6.518116265389094e-003,
  123.  -8.304767428789135e-001,
  124.  -2.820782997470719e-003,
  125.  -2.248295252614748e-001,
  126.  
  127. -3.217557586964646e-003, /* Saturn */
  128.  7.894394222587950e+000,
  129.   4.335809727271375e-003,
  130.   4.596483181197301e+000,
  131.   1.928645932522133e-003,
  132.   1.558697670343092e+000,
  133.  
  134.  2.211920155513113e-004, /* Uranus */
  135. -1.826538634212984e+001,
  136.  -3.762477278855801e-003,
  137.  -1.161959795438843e+000,
  138.  -1.651014777644869e-003,
  139.  -2.501080003472951e-001,
  140.  
  141.  2.642773506889557e-003, /* Neptune */
  142. -1.605499953831724e+001,
  143.  -1.498309226803948e-003,
  144.  -2.394216805154041e+001,
  145.  -6.790394433396178e-004,
  146.  -9.400159225576589e+000,
  147.  
  148.  3.221893230538327e-004, /* Pluto */
  149. -3.048349217668479e+001,
  150.  -3.143582327761257e-003,
  151.  -8.724432817555443e-001,
  152.  -1.077956399281471e-003,
  153.   8.911620591248875e+000,
  154. #if MOON
  155.  6.010848314829119e-004,
  156. -8.081772358351251e-004,
  157.  -1.674454691500606e-004,
  158.  -1.994630037441996e-003,
  159.  -8.556208109904861e-005,
  160.  -1.087262721620868e-003,
  161. #endif
  162. -3.524457445683394e-007, /* Sun */
  163.  4.504795855674601e-003,
  164.   5.177637780672221e-006,
  165.   7.732544746890742e-004,
  166.   2.229113252400401e-006,
  167.   2.685039985573271e-004
  168. };
  169.  
  170.  
  171. /* Julian date of the state yn1[]: */
  172. #define JD1 2440500.50
  173.  
  174. static double yn1[] = {
  175. -2.457085900366332e-002,
  176.  2.381597794693489e-001,
  177.   1.851778146680485e-002,
  178.   1.996323209474686e-001,
  179.   1.244039124191612e-002,
  180.   8.216842968554781e-002,
  181.  
  182. -1.619513532964599e-002,
  183. -4.288838589850871e-001,
  184.  -1.160656795724704e-002,
  185.   5.131125475173164e-001,
  186.  -4.195286108659672e-003,
  187.   2.581368657564266e-001,
  188.  
  189. -4.148833933015173e-003,
  190.  9.784941438766318e-001,
  191.   1.532745170011609e-002,
  192.   2.074410311959717e-001,
  193.   6.646462310684277e-003,
  194.   8.989091702317244e-002,
  195.  
  196.  8.350422878366247e-003,
  197.  1.151108827619137e+000,
  198.   1.172723634265029e-002,
  199.  -6.894717355437423e-001,
  200.   5.152260844982779e-003,
  201.  -3.474143803147213e-001,
  202.  
  203.  2.061461649176091e-003,
  204. -5.221785914203378e+000,
  205.  -6.307486608326940e-003,
  206.  -1.472733759859347e+000,
  207.  -2.754121811344123e-003,
  208.  -5.039994787357370e-001,
  209.  
  210. -3.506509603502537e-003,
  211.  7.558106559699063e+000,
  212.   4.155902532833776e-003,
  213.   5.021213850462841e+000,
  214.   1.866781421893536e-003,
  215.   1.748532556433840e+000,
  216.  
  217.  3.093738240936483e-004,
  218. -1.823885644320201e+001,
  219.  -3.755962489549436e-003,
  220.  -1.537896840580550e+000,
  221.  -1.649410160709415e-003,
  222.  -4.151358646421778e-001,
  223.  
  224.  2.659694712570950e-003,
  225. -1.578987378978514e+001,
  226.  -1.472781001534497e-003,
  227.  -2.409072389641744e+001,
  228.  -6.690116571980626e-004,
  229.  -9.467562384490781e+000,
  230.  
  231.  3.504107115900905e-004,
  232. -3.044986307301521e+001,
  233.  -3.142629755658865e-003,
  234.  -1.186756347889178e+000,
  235.  -1.086164460786462e-003,
  236.   8.803414050220775e+000,
  237. #if MOON
  238. -4.181662799524538e-004,
  239. -1.749747307712349e-003,
  240.  -3.263107309648734e-004,
  241.   1.807047780153076e-003,
  242.  -1.847196922506988e-004,
  243.   9.595533885666276e-004,
  244. #endif
  245. -1.063378525544779e-006,
  246.  4.434825780648206e-003,
  247.   5.049061922605073e-006,
  248.   1.283313197720549e-003,
  249.   2.188756349489559e-006,
  250.   4.887462262250162e-004
  251. };
  252.  
  253.  
  254. /* Earth's mass divided by Moon's mass
  255.  */
  256. #define EMRAT 81.300587
  257.  
  258.  
  259. /* GM's of the solar system bodies
  260.  * These are scaled such that GMsun = k^2
  261.  * (k = Gaussian gravitational constant).
  262.  */
  263. double GMs[] = {
  264.  4.912547451450812e-011, /* Mercury */
  265.  7.243456209632766e-010, /* Venus */
  266. #if MOON
  267.  (8.997011658557308e-010*EMRAT)/(1.0+EMRAT), /* Earth */
  268. #else
  269.  8.997011658557308e-010, /* EMB */
  270. #endif
  271.  9.549528942224058e-011, /* Mars */
  272.  2.825342103445926e-007, /* Jupiter */
  273.  8.459468504830660e-008, /* Saturn */
  274.  1.288816238138035e-008, /* Uranus */
  275.  1.532112481284276e-008, /* Neptune */
  276.  2.276247751863699e-012, /* Pluto */
  277. #if MOON
  278.  (8.997011658557308e-010)/(1.0+EMRAT), /* Moon */
  279. #endif
  280.  2.959122082855911e-004 /* Sun */
  281. };
  282.  
  283.  
  284. main()
  285. {
  286. double t, e0, e, err, h, ccor;
  287. double *pdv;
  288. int i, ii, j, nsteps, ord;
  289. double adstep();
  290.  
  291.  
  292. orlup:
  293. printf( "Adams order ? " );
  294. scanf( "%d", &ord );
  295. if( (ord > MAXORD) || (ord < 1) )
  296.     {
  297.     printf( "order must be between 1 and %d\n", MAXORD );
  298.     goto orlup;
  299.     }
  300.  
  301. printf( "step size, days ? " );
  302. scanf( "%lf", &h );
  303.  
  304. printf( "Number of steps ? " );
  305. scanf( "%d", &nsteps );
  306.  
  307. for( i=0; i<((MAXRUNG+2)*NEQ); i++ )
  308.     rwork[i] = 0.0;
  309.  
  310. printf( "initializing...\n" );
  311. t = JD0;
  312. err = 0.0;
  313. rkstart( NEQ, rwork );
  314. adstart( h, yn, awork, NEQ, ord, t );
  315. printf( "initialized.\n" );
  316.  
  317.  
  318. for( j=1; j<=nsteps; j++ )
  319.     {
  320.     err += adstep( &t, yn, NEQC );
  321.     printf( "%5d %11.2lf %.3e\n", j, t, err );
  322.     }
  323.  
  324. printf( "Final x, y, z, positions and errors:\n" );
  325. ii = 0;
  326. for (i=0; i<NTOTAL; i++)
  327.     {
  328.     printf("%19.15lf  %19.15lf  %19.15lf\n",
  329.         yn[ii+1], yn[ii+3], yn[ii+5] );
  330.     printf("%19.3e    %19.3e    %19.3e\n",
  331.         yn[ii+1] - yn1[ii+1],
  332.         yn[ii+3] - yn1[ii+3],
  333.         yn[ii+5] - yn1[ii+5] );
  334.     ii += 6;
  335.     }
  336. }
  337.  
  338.  
  339.  
  340. /* Subroutine func() calculates the forces and accelerations.
  341.  * Data in the output vector v[] are in the order
  342.  * d^2x/dt^2, dx/dt, d^2y/dt^2, dy/dt, d^2z/dt^2, dz/dt
  343.  * for each object.  For this problem the velocities dx/dt, ...
  344.  * are simply copied over from the input array y[].
  345.  */
  346.  
  347. #if MOON
  348. #define yxx yin
  349. #else
  350. #define  yxx y
  351. #endif
  352.  
  353. func( t, yxx, v )
  354. double t; /* time */
  355. double yxx[]; /* input state: velocity and position */
  356. double v[]; /* output: calculated acceleration, copy of input velocity */
  357. {
  358. int i, j, k;
  359. int ii, jj;
  360. double xs, ys, zs;
  361. double xd, yd, zd, r, s, csqi;
  362. double temp, rc, vi;
  363. double xv, yv, zv;
  364.  
  365. /* Copy input to temp and unravel earth/moon
  366.  */
  367. #if MOON
  368. static double xm, ym, zm;
  369. static double y[NEQ];
  370.  
  371. for( i=0; i<NEQ; i++ )
  372.     y[i] = yin[i];
  373. /* Locally replace input variable EMB by barycentric earth
  374.  * and input variable M by barycentric Moon.
  375.  */
  376. xm = yin[55]; /* M */
  377. ym = yin[57];
  378. zm = yin[59];
  379. ii = 12;
  380. jj = 54;
  381. for( i=0; i<6; i++ )
  382.     {
  383.     zd = yin[ii+i];            /* EMB */
  384.     yd = yin[jj+i]/(1.0+EMRAT);    /* M */
  385.     y[ii+i] = zd -  yd;        /* r_e */
  386.     y[jj+i] = zd + EMRAT * yd;    /* r_m */
  387.     }
  388. #endif
  389.  
  390. /* Constrain the barycenter to stay at the origin.
  391.  * This is done by adjusting the Sun, which body is then
  392.  * not included with the variables that are integrated.
  393.  */
  394. #if DOREL
  395. #if BARYC
  396.     fixsun( y, v );
  397. #endif
  398. #endif
  399.  
  400. /* Table of distances between objects i and j */
  401. jj = 0;
  402. for (j=0; j<NTOTAL; j++)
  403.     {
  404.     ii = 0;
  405.     xv = y[jj+1];
  406.     yv = y[jj+3];
  407.     zv = y[jj+5];
  408.     for (i=0; i<j; i++)
  409.         {
  410.         xd = xv - y[ii+1];
  411.         yd = yv - y[ii+3];
  412.         zd = zv - y[ii+5];
  413.         r = sqrt(xd*xd + yd*yd + zd*zd);
  414.         Rij[j][i] = r;
  415.         Rij[i][j]  = r;
  416.         ii += 6;
  417.         }
  418.     jj += 6;
  419.     }
  420.  
  421. #if MOON
  422. /* Take the input M vector for distance from Earth to Moon. */
  423. r = sqrt( xm*xm + ym*ym + zm*zm );
  424. Rij[2][9] = r;
  425. Rij[9][2] = r;
  426. #endif
  427.  
  428. /* Compute Newtonian acceleration. */
  429. ii = 0;
  430. for( i=0; i<NTOTAL; i++ )
  431.     {
  432.     xs = 0.0;
  433.     ys = 0.0;
  434.     zs = 0.0;
  435.     xv = y[ii+1];
  436.     yv = y[ii+3];
  437.     zv = y[ii+5];
  438.     jj = 0;
  439.     for( j=0; j<NTOTAL; j++ )
  440.         {
  441.         if( j != i )
  442.             {
  443.             xd = y[jj+1] - xv;
  444.             yd = y[jj+3] - yv;
  445.             zd = y[jj+5] - zv;
  446.             r = Rij[i][j];
  447.             temp = GMs[j]/(r * r * r);
  448.             xs += xd * temp;
  449.             ys += yd * temp;
  450.             zs += zd * temp;
  451.             }
  452.         jj += 6;
  453.         }
  454.     vnewt[ii] = xs; /* acceleration */
  455.     vnewt[ii+2] = ys;
  456.     vnewt[ii+4] = zs;
  457.     vnewt[ii+1] = y[ii]; /* velocity */
  458.     vnewt[ii+3] = y[ii+2];
  459.     vnewt[ii+5] = y[ii+4];
  460.     ii += 6;
  461.     }
  462.  
  463. #if DOREL
  464.  
  465. /* Relativistic corrections.  Reference:
  466.  *
  467.  * Newhall, XX, E. M. Standish, and J. G. Williams, "DE102: a
  468.  * numerically integrated ephemeris of the Moon and planets
  469.  * spanning forty-four centuries," Astronomy and Astrophysics
  470.  * 125, 150-167 (1983).
  471.  */
  472.     
  473. csqi = 1.0/(C*C);
  474. ii = 0;
  475. for (i=0; i<NBODY; i++)
  476.     {
  477.     v[ii] = 0.0;
  478.     v[ii+2] = 0.0;
  479.     v[ii+4] = 0.0;
  480.  
  481.     xv = y[ii]; /* velocity */
  482.     yv = y[ii+2];
  483.     zv = y[ii+4];
  484.     vi = xv*xv + yv*yv + zv*zv;
  485.  
  486.     xs = 0.0;
  487.     ys = 0.0;
  488.     zs = 0.0;
  489.     jj = 0;
  490.     for (j=0; j<NTOTAL; j++)
  491.         {
  492.         if( j == i )
  493.             {
  494.             jj += 6;
  495.             continue; /* skip to next j if i = j */
  496.             }
  497.  
  498.         xd = y[jj+1] - y[ii+1];
  499.         yd = y[jj+3] - y[ii+3];
  500.         zd = y[jj+5] - y[ii+5];
  501.  
  502.         s = 0.0;
  503.         for (k=0; k<NTOTAL; k++)
  504.             {
  505.             if( k == i )
  506.                 continue;
  507.             s += GMs[k]/Rij[i][k];
  508.             }
  509.         rc = -4.0 * csqi * s;
  510.  
  511.         s = 0.0;
  512.         for (k=0; k<NTOTAL; k++)
  513.             {
  514.             if(k ==j )
  515.                 continue;
  516.             s += GMs[k]/Rij[j][k];
  517.             }
  518.         rc -= csqi * s;
  519.  
  520.         rc += vi * csqi;
  521.  
  522.         xv = y[jj];
  523.         yv = y[jj+2];
  524.         zv = y[jj+4];
  525.  
  526.         r = xv * xv + yv * yv + zv * zv;
  527.         rc += 2.0 * csqi * r;
  528.  
  529.         r = y[ii] * y[jj] + y[ii+2] * y[jj+2] + y[ii+4] * y[jj+4];
  530.         rc -= 4.0 * csqi * r;
  531.  
  532.         s = -xd * y[jj] - yd * y[jj+2] - zd * y[jj+4];
  533.         s /= Rij[i][j];
  534.         rc -= 1.5 * csqi * s * s;
  535.  
  536.         s = xd * vnewt[jj] + yd * vnewt[jj+2] + zd * vnewt[jj+4];
  537.         rc += 0.5 * csqi * s;
  538.  
  539.         r = Rij[i][j];
  540.         temp = GMs[j]/(r * r * r);
  541.  
  542.         rc = temp * (1.0 + rc );
  543.         xs += xd * rc;
  544.         ys += yd * rc;
  545.         zs += zd * rc;
  546.  
  547.         s = -xd * (4.0*y[ii] - 3.0*y[jj])
  548.             - yd * (4.0*y[ii+2] - 3.0*y[jj+2])
  549.             - zd * (4.0*y[ii+4] - 3.0*y[jj+4]);
  550.         s = s * temp * csqi;
  551.         xs += s * (y[ii] - y[jj]);
  552.         ys += s * (y[ii+2] - y[jj+2]);
  553.         zs += s * (y[ii+4] - y[jj+4]);
  554.  
  555.         temp = 3.5 * csqi * GMs[j] / r;
  556.         xs += temp * vnewt[jj];
  557.         ys += temp * vnewt[jj+2];
  558.         zs += temp * vnewt[jj+4];
  559.         jj += 6;
  560.         }
  561.     v[ii] = xs;
  562.     v[ii+2] = ys;
  563.     v[ii+4] = zs;
  564.     v[ii+1] = y[ii];
  565.     v[ii+3] = y[ii+2];
  566.     v[ii+5] = y[ii+4];
  567.     ii += 6;
  568.     }
  569. #else
  570. /* No relativity theory. Return the Newtonian accelerations. */
  571. ii = 0;
  572. for( i=0; i<NTOTAL; i++ )
  573.     {
  574.     v[ii] = vnewt[ii];
  575.     v[ii+2] = vnewt[ii+2];
  576.     v[ii+4] = vnewt[ii+4];
  577.     v[ii+1] = y[ii];
  578.     v[ii+3] = y[ii+2];
  579.     v[ii+5] = y[ii+4];
  580.     ii += 6;
  581.     }
  582. #endif
  583.  
  584. #if MOON
  585. /* Convert barycentric Earth and Moon to output EMB and M variables. */
  586. ii = 12;
  587. jj = 54;
  588. for( i=0; i<6; i += 2 )
  589.     {
  590.     xd = v[ii+i]; /* Earth */
  591.     yd = v[jj+i]; /* Moon */
  592.     v[ii+i] = (EMRAT * xd + yd)/(EMRAT+1.0); /* EMB */
  593.     v[jj+i] = yd - xd; /* M = Moon - Earth */
  594.     v[ii+i+1] = yin[ii+i];
  595.     v[jj+i+1] = yin[jj+i];
  596.     }
  597. #endif
  598. }
  599.  
  600. #if DOREL
  601. /* Constraint that the center of the relativistic masses = zero.
  602.  */
  603. static double ysun[6];
  604.  
  605. fixsun( y, v )
  606. double y[], v[];
  607. {
  608. double xx, yy, zz, vx, vy, vz, ax, ay, az;
  609. double csqi, s;
  610. double mustar[NTOTAL];
  611. int i, j, k, ii, jj;
  612.  
  613. for( ii=0; ii<6; ii++ )
  614.     ysun[ii] = y[ii+(6*ISUN)];
  615.  
  616. csqi = 0.5 / (C*C);
  617. for( k=0; k<2; k++ )
  618. { /* Iterate to find solution of implicit expressions. */
  619.  
  620. /* Table of distances between bodies i and j.
  621.  * Note, most of this work may already have been done by func().
  622.  */
  623. ii = 6;
  624. for( i=1; i<NTOTAL; i++ )
  625.     {
  626.     jj = 0;
  627.     vx = y[ii+1]; /* position */
  628.     vy = y[ii+3];
  629.     vz = y[ii+5];
  630.     for( j=0; j<i; j++ )
  631.         {
  632.         xx = vx - y[jj+1];
  633.         yy = vy - y[jj+3];
  634.         zz = vz - y[jj+5];
  635.         xx = sqrt( xx*xx + yy*yy + zz*zz );
  636.         Rij[i][j] = xx;
  637.         Rij[j][i] = xx;
  638.         jj += 6;
  639.         }
  640.     ii += 6;
  641.     }
  642. /* Relativistic GM of each body */
  643. ii = 0;
  644. for( i=0; i<NTOTAL; i++ )
  645.     {
  646.     vx = y[ii]; /* velocity */
  647.     vy = y[ii+2];
  648.     vz = y[ii+4];
  649.     s = vx * vx + vy * vy + vz * vz; /* square of velocity */
  650.     for( j=0; j<NTOTAL; j++ )
  651.         {
  652.         if( j == i )
  653.             continue;
  654.         s -= GMs[j]/Rij[i][j];
  655.         }
  656.     mustar[i] = GMs[i] * (1.0 + csqi * s);
  657.     ii += 6;
  658.     }
  659. /* Compute center of mass with Sun omitted. */
  660. xx = 0.0;
  661. yy = 0.0;
  662. zz = 0.0;
  663. vx = 0.0;
  664. vy = 0.0;
  665. vz = 0.0;
  666. ii = 0;
  667. for( i=0; i<ISUN; i++ )
  668.     {
  669.     s = mustar[i];
  670.     xx += y[ii+1] * s; /* position */
  671.     yy += y[ii+3] * s;
  672.     zz += y[ii+5] * s;
  673.     vx += y[ii] * s; /* velocity */
  674.     vy += y[ii+2] * s;
  675.     vz += y[ii+4] * s;
  676.     ii += 6;
  677.     }
  678. /* Evaluate the Sun so that the center = 0. */
  679. s = mustar[ISUN];
  680. jj = 6*ISUN;
  681. y[jj+1] = -xx/s;
  682. y[jj+3] = -yy/s;
  683. y[jj+5] = -zz/s;
  684. y[jj] = -vx/s;
  685. y[jj+2] = -vy/s;
  686. y[jj+4] = -vz/s;
  687. v[jj+1] = y[jj];
  688. v[jj+3] = y[jj+2];
  689. v[jj+5] = y[jj+4];
  690. }
  691.  
  692. /* debug
  693. for( ii=0; ii<6; ii++ )
  694.     {
  695.     xx = ysun[ii] - y[ii+6*ISUN];
  696.     printf( "%.1e\n", xx );
  697.     }
  698. */
  699. }
  700. #endif
  701.  
  702.