/* almabench.c 26 March 2004 Written by Scott Robert Ladd. No rights reserved. This is public domain software, for use by anyone. A number-crunching benchmark that can be used as a fitness test for evolving optimal compiler options via genetic algorithm. This program calculates the daily ephemeris (at noon) for the years 2000-2099 using an algorithm developed by J.L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touze, G. Francou and J. Laskar of the Bureau des Longitudes, Paris, France), as detailed in Astronomy & Astrophysics 282, 663 (1994) This program's performance is heavily predicated on the compiler's effectiveness in generating the sin, cos, and sqrt functions. In some cases, this benchmark's performance may be adversely affected by your C library's implementation of those functions. A good compiler, in "fastest math" mode, will generate inline calls to processor instructions for these fucntions. Note that the code herein is design for the purpose of testing computational performance; error handling and other such "niceties" is virtually non-existent. Actual benchmark results can be found at: http://www.coyotegulch.com Please do not use this information or algorithm in any way that might upset the balance of the universe or otherwise cause planets to impact upon one another. */ #include #include #include #include #include #define CALC_PI 3.14159265358979323846 static const double PI = CALC_PI; static const double J2000 = 2451545.0; static const double JCENTURY = 36525.0; static const double JMILLENIA = 365250.0; static const double TWOPI = 2.0 * CALC_PI; static const double A2R = CALC_PI / 648000.0; static const double R2H = 12.0 / CALC_PI; static const double R2D = 180.0 / CALC_PI; static const double GAUSSK = 0.01720209895; // number of loops to perform static const int TEST_LOOPS = 2; // number of days to include in test static const int TEST_DAYS = 36525; // sin and cos of j2000 mean obliquity (iau 1976) static const double sineps = 0.3977771559319137; static const double coseps = 0.9174820620691818; static const double amas [8] = { 6023600.0, 408523.5, 328900.5, 3098710.0, 1047.355, 3498.5, 22869.0, 19314.0 }; // tables giving the mean keplerian elements, limited to t**2 terms: // a semi-major axis (au) // dlm mean longitude (degree and arcsecond) // e eccentricity // pi longitude of the perihelion (degree and arcsecond) // dinc inclination (degree and arcsecond) // omega longitude of the ascending node (degree and arcsecond) static const double a [8][3] = { { 0.3870983098, 0, 0 }, { 0.7233298200, 0, 0 }, { 1.0000010178, 0, 0 }, { 1.5236793419, 3e-10, 0 }, { 5.2026032092, 19132e-10, -39e-10 }, { 9.5549091915, -0.0000213896, 444e-10 }, { 19.2184460618, -3716e-10, 979e-10 }, { 30.1103868694, -16635e-10, 686e-10 } }; static const double dlm[8][3] = { { 252.25090552, 5381016286.88982, -1.92789 }, { 181.97980085, 2106641364.33548, 0.59381 }, { 100.46645683, 1295977422.83429, -2.04411 }, { 355.43299958, 689050774.93988, 0.94264 }, { 34.35151874, 109256603.77991, -30.60378 }, { 50.07744430, 43996098.55732, 75.61614 }, { 314.05500511, 15424811.93933, -1.75083 }, { 304.34866548, 7865503.20744, 0.21103 } }; static const double e[8][3] = { { 0.2056317526, 0.0002040653, -28349e-10 }, { 0.0067719164, -0.0004776521, 98127e-10 }, { 0.0167086342, -0.0004203654, -0.0000126734 }, { 0.0934006477, 0.0009048438, -80641e-10 }, { 0.0484979255, 0.0016322542, -0.0000471366 }, { 0.0555481426, -0.0034664062, -0.0000643639 }, { 0.0463812221, -0.0002729293, 0.0000078913 }, { 0.0094557470, 0.0000603263, 0 } }; static const double pi[8][3] = { { 77.45611904, 5719.11590, -4.83016 }, { 131.56370300, 175.48640, -498.48184 }, { 102.93734808, 11612.35290, 53.27577 }, { 336.06023395, 15980.45908, -62.32800 }, { 14.33120687, 7758.75163, 259.95938 }, { 93.05723748, 20395.49439, 190.25952 }, { 173.00529106, 3215.56238, -34.09288 }, { 48.12027554, 1050.71912, 27.39717 } }; static const double dinc[8][3] = { { 7.00498625, -214.25629, 0.28977 }, { 3.39466189, -30.84437, -11.67836 }, { 0, 469.97289, -3.35053 }, { 1.84972648, -293.31722, -8.11830 }, { 1.30326698, -71.55890, 11.95297 }, { 2.48887878, 91.85195, -17.66225 }, { 0.77319689, -60.72723, 1.25759 }, { 1.76995259, 8.12333, 0.08135 } }; static const double omega[8][3] = { { 48.33089304, -4515.21727, -31.79892 }, { 76.67992019, -10008.48154, -51.32614 }, { 174.87317577, -8679.27034, 15.34191 }, { 49.55809321, -10620.90088, -230.57416 }, { 100.46440702, 6362.03561, 326.52178 }, { 113.66550252, -9240.19942, -66.23743 }, { 74.00595701, 2669.15033, 145.93964 }, { 131.78405702, -221.94322, -0.78728 } }; // tables for trigonometric terms to be added to the mean elements // of the semi-major axes. static const double kp[8][9] = { { 69613.0, 75645.0, 88306.0, 59899.0, 15746.0, 71087.0, 142173.0, 3086.0, 0.0 }, { 21863.0, 32794.0, 26934.0, 10931.0, 26250.0, 43725.0, 53867.0, 28939.0, 0.0 }, { 16002.0, 21863.0, 32004.0, 10931.0, 14529.0, 16368.0, 15318.0, 32794.0, 0.0 }, { 6345.0, 7818.0, 15636.0, 7077.0, 8184.0, 14163.0, 1107.0, 4872.0, 0.0 }, { 1760.0, 1454.0, 1167.0, 880.0, 287.0, 2640.0, 19.0, 2047.0, 1454.0 }, { 574.0, 0.0, 880.0, 287.0, 19.0, 1760.0, 1167.0, 306.0, 574.0 }, { 204.0, 0.0, 177.0, 1265.0, 4.0, 385.0, 200.0, 208.0, 204.0 }, { 0.0, 102.0, 106.0, 4.0, 98.0, 1367.0, 487.0, 204.0, 0.0 } }; static const double ca[8][9] = { { 4.0, -13.0, 11.0, -9.0, -9.0, -3.0, -1.0, 4.0, 0.0 }, { -156.0, 59.0, -42.0, 6.0, 19.0, -20.0, -10.0, -12.0, 0.0 }, { 64.0, -152.0, 62.0, -8.0, 32.0, -41.0, 19.0, -11.0, 0.0 }, { 124.0, 621.0, -145.0, 208.0, 54.0, -57.0, 30.0, 15.0, 0.0 }, { -23437.0, -2634.0, 6601.0, 6259.0, -1507.0, -1821.0, 2620.0, -2115.0,-1489.0 }, { 62911.0,-119919.0, 79336.0, 17814.0,-24241.0, 12068.0, 8306.0, -4893.0, 8902.0 }, { 389061.0,-262125.0,-44088.0, 8387.0,-22976.0, -2093.0, -615.0, -9720.0, 6633.0 }, { -412235.0,-157046.0,-31430.0, 37817.0, -9740.0, -13.0, -7449.0, 9644.0, 0.0 } }; static const double sa[8][9] = { { -29.0, -1.0, 9.0, 6.0, -6.0, 5.0, 4.0, 0.0, 0.0 }, { -48.0, -125.0, -26.0, -37.0, 18.0, -13.0, -20.0, -2.0, 0.0 }, { -150.0, -46.0, 68.0, 54.0, 14.0, 24.0, -28.0, 22.0, 0.0 }, { -621.0, 532.0, -694.0, -20.0, 192.0, -94.0, 71.0, -73.0, 0.0 }, { -14614.0,-19828.0, -5869.0, 1881.0, -4372.0, -2255.0, 782.0, 930.0, 913.0 }, { 139737.0, 0.0, 24667.0, 51123.0, -5102.0, 7429.0, -4095.0, -1976.0,-9566.0 }, { -138081.0, 0.0, 37205.0,-49039.0,-41901.0,-33872.0,-27037.0,-12474.0,18797.0 }, { 0.0, 28492.0,133236.0, 69654.0, 52322.0,-49577.0,-26430.0, -3593.0, 0.0 } }; // tables giving the trigonometric terms to be added to the mean // elements of the mean longitudes. static const double kq[8][10] = { { 3086.0, 15746.0, 69613.0, 59899.0, 75645.0, 88306.0, 12661.0, 2658.0, 0.0, 0.0 }, { 21863.0, 32794.0, 10931.0, 73.0, 4387.0, 26934.0, 1473.0, 2157.0, 0.0, 0.0 }, { 10.0, 16002.0, 21863.0, 10931.0, 1473.0, 32004.0, 4387.0, 73.0, 0.0, 0.0 }, { 10.0, 6345.0, 7818.0, 1107.0, 15636.0, 7077.0, 8184.0, 532.0, 10.0, 0.0 }, { 19.0, 1760.0, 1454.0, 287.0, 1167.0, 880.0, 574.0, 2640.0, 19.0,1454.0 }, { 19.0, 574.0, 287.0, 306.0, 1760.0, 12.0, 31.0, 38.0, 19.0, 574.0 }, { 4.0, 204.0, 177.0, 8.0, 31.0, 200.0, 1265.0, 102.0, 4.0, 204.0 }, { 4.0, 102.0, 106.0, 8.0, 98.0, 1367.0, 487.0, 204.0, 4.0, 102.0 } }; static const double cl[8][10] = { { 21.0, -95.0, -157.0, 41.0, -5.0, 42.0, 23.0, 30.0, 0.0, 0.0 }, { -160.0, -313.0, -235.0, 60.0, -74.0, -76.0, -27.0, 34.0, 0.0, 0.0 }, { -325.0, -322.0, -79.0, 232.0, -52.0, 97.0, 55.0, -41.0, 0.0, 0.0 }, { 2268.0, -979.0, 802.0, 602.0, -668.0, -33.0, 345.0, 201.0, -55.0, 0.0 }, { 7610.0, -4997.0,-7689.0,-5841.0,-2617.0, 1115.0, -748.0, -607.0, 6074.0, 354.0 }, { -18549.0, 30125.0,20012.0, -730.0, 824.0, 23.0, 1289.0, -352.0,-14767.0,-2062.0 }, { -135245.0,-14594.0, 4197.0,-4030.0,-5630.0,-2898.0, 2540.0, -306.0, 2939.0, 1986.0 }, { 89948.0, 2103.0, 8963.0, 2695.0, 3682.0, 1648.0, 866.0, -154.0, -1963.0, -283.0 } }; static const double sl[8][10] = { { -342.0, 136.0, -23.0, 62.0, 66.0, -52.0, -33.0, 17.0, 0.0, 0.0 }, { 524.0, -149.0, -35.0, 117.0, 151.0, 122.0, -71.0, -62.0, 0.0, 0.0 }, { -105.0, -137.0, 258.0, 35.0, -116.0, -88.0, -112.0, -80.0, 0.0, 0.0 }, { 854.0, -205.0, -936.0, -240.0, 140.0, -341.0, -97.0, -232.0, 536.0, 0.0 }, { -56980.0, 8016.0, 1012.0, 1448.0,-3024.0,-3710.0, 318.0, 503.0, 3767.0, 577.0 }, { 138606.0,-13478.0,-4964.0, 1441.0,-1319.0,-1482.0, 427.0, 1236.0, -9167.0,-1918.0 }, { 71234.0,-41116.0, 5334.0,-4935.0,-1848.0, 66.0, 434.0,-1748.0, 3780.0, -701.0 }, { -47645.0, 11647.0, 2166.0, 3194.0, 679.0, 0.0, -244.0, -419.0, -2531.0, 48.0 } }; //--------------------------------------------------------------------------- // Normalize angle into the range -pi <= A < +pi. double anpm (double a) { double w = fmod(a,TWOPI); if (fabs(w) >= PI) w = w - ((a < 0) ? -TWOPI : TWOPI); return w; } //--------------------------------------------------------------------------- // The reference frame is equatorial and is with respect to the // mean equator and equinox of epoch j2000. void planetpv (double epoch[2], int np, double pv[2][3]) { // working storage int i, j, k; double t, da, dl, de, dp, di, doh, dmu, arga, argl, am; double ae, dae, ae2, at, r, v, si2, xq, xp, tl, xsw; double xcw, xm2, xf, ci2, xms, xmc, xpxq2, x, y, z; // time: julian millennia since j2000. t = ((epoch[0] - J2000) + epoch[1]) / JMILLENIA; // compute the mean elements. da = a[np][0] + (a[np][1] + a[np][2] * t ) * t; dl = (3600.0 * dlm[np][0] + (dlm[np][1] + dlm[np][2] * t ) * t ) * A2R; de = e[np][0] + (e[np][1] + e[np][2] * t ) * t; dp = anpm((3600.0 * pi[np][0] + (pi[np][1] + pi[np][2] * t ) * t ) * A2R ); di = (3600.0 * dinc[np][0] + (dinc[np][1] + dinc[np][2] * t ) * t ) * A2R; doh = anpm((3600.0 * omega[np][0] + (omega[np][1] + omega[np][2] * t ) * t ) * A2R ); // apply the trigonometric terms. dmu = 0.35953620 * t; for (k = 0; k < 8; ++k) { arga = kp[np][k] * dmu; argl = kq[np][k] * dmu; da = da + (ca[np][k] * cos(arga) + sa[np][k] * sin(arga)) * 0.0000001; dl = dl + (cl[np][k] * cos(argl) + sl[np][k] * sin(argl)) * 0.0000001; } arga = kp[np][8] * dmu; da = da + t * (ca[np][8] * cos(arga) + sa[np][8] * sin(arga)) * 0.0000001; for (k = 8; k <= 9; ++k) { argl = kq[np][k] * dmu; dl = dl + t * ( cl[np][k] * cos(argl) + sl[np][k] * sin(argl) ) * 0.0000001; } dl = fmod(dl,TWOPI); // iterative solution of kepler's equation to get eccentric anomaly. am = dl - dp; ae = am + de * sin(am); k = 0; while (1) { dae = (am - ae + de * sin(ae)) / (1.0 - de * cos(ae)); ae = ae + dae; k = k + 1; if ((k >= 10) || (fabs(dae) < 1e-12)) break; } // true anomaly. ae2 = ae / 2.0; at = 2.0 * atan2(sqrt((1.0 + de)/(1.0 - de)) * sin(ae2), cos(ae2)); // distance (au) and speed (radians per day). r = da * (1.0 - de * cos(ae)); v = GAUSSK * sqrt((1.0 + 1.0 / amas[np] ) / (da * da * da)); si2 = sin(di / 2.0); xq = si2 * cos(doh); xp = si2 * sin(doh); tl = at + dp; xsw = sin(tl); xcw = cos(tl); xm2 = 2.0 * (xp * xcw - xq * xsw ); xf = da / sqrt(1.0 - de*de); ci2 = cos(di / 2.0); xms = (de * sin(dp) + xsw) * xf; xmc = (de * cos(dp) + xcw) * xf; xpxq2 = 2.0 * xp * xq; // position (j2000 ecliptic x,y,z in au). x = r * (xcw - xm2 * xp); y = r * (xsw + xm2 * xq); z = r * (-xm2 * ci2); // rotate to equatorial. pv[0][0] = x; pv[0][1] = y * coseps - z * sineps; pv[0][2] = y * sineps + z * coseps; // velocity (j2000 ecliptic xdot,ydot,zdot in au/d). x = v * ((-1.0 + 2.0 * xp * xp) * xms + xpxq2 * xmc); y = v * (( 1.0 - 2.0 * xq * xq ) * xmc - xpxq2 * xms); z = v * (2.0 * ci2 * (xp * xms + xq * xmc)); // rotate to equatorial. pv[1][0] = x; pv[1][1] = y * coseps - z * sineps; pv[1][2] = y * sineps + z * coseps; } //--------------------------------------------------------------------------- // Computes RA, Declination, and distance from a state vector returned by // planetpv. void radecdist(double state[2][3], double rdd[3]) { // distance rdd[2] = sqrt(state[0][0] * state[0][0] + state[0][1] * state[0][1] + state[0][2] * state[0][2]); // RA rdd[0] = atan2(state[0][1], state[0][0]) * R2H; if (rdd[0] < 0.0) rdd[0] += 24.0; // Declination rdd[1] = asin(state[0][2] / rdd[2]) * R2D; } //--------------------------------------------------------------------------- // Entry point // Calculate RA and Dec for noon on every day in 1900-2100 int main(int argc, char ** argv) { int i, n, p; double jd[2]; double pv[2][3]; double position[3]; bool ga_testing = false; // do we have verbose output? if (argc > 1) { for (i = 1; i < argc; ++i) { if (!strcmp(argv[1],"-ga")) { ga_testing = true; break; } } } // get starting time struct timespec start, stop; clock_gettime(CLOCK_REALTIME,&start); // main loop for (i = 0; i < TEST_LOOPS; ++i) { jd[0] = J2000; jd[1] = 0.0; for (n = 0; n < TEST_DAYS; ++n) { jd[0] += 1.0; for (p = 0; p < 8; ++p) { planetpv(jd,p,pv); radecdist(pv,position); } } } // get final time clock_gettime(CLOCK_REALTIME,&stop); double run_time = (stop.tv_sec - start.tv_sec) + (double)(stop.tv_nsec - start.tv_nsec) / 1000000000.0; // report runtime if (ga_testing) fprintf(stdout,"%f",run_time); else fprintf(stdout,"\nalmabench (Std. C) run time: %f\n\n",run_time); fflush(stdout); return 0; }