rendered paste body// ALMABENCH 1.0.1// Standard C++ version//// A number-crunching benchmark designed for cross-language and vendor// comparisons.//// Written by Scott Robert Ladd.//// No rights reserved. This is public domain software, for use by anyone.//// 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)//// 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 (and oft-updated) 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.// Added Timer() code#include <cstdlib>#include <ctime>#include <cmath>#include <iostream>using namespace std;static const double PI = 3.14159265358979323846;static const double J2000 = 2451545.0;static const double JCENTURY = 36525.0;static const double JMILLENIA = 365250.0;static const double TWOPI = 2.0 * PI;static const double A2R = PI / 648000.0;static const double R2H = 12.0 / PI;static const double R2D = 180.0 / PI;static const double GAUSSK = 0.01720209895;// number of days to include in teststatic const int TEST_LOOPS = 20;static const int TEST_LENGTH = 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 (true) { 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-2100int main(){ clock_t start=clock(); int i, n, p; double jd[2]; double pv[2][3]; double position[3]; //#pragma omp parallel for private(i, n, p, jd, pv, position) for (i = 0; i < TEST_LOOPS; ++i) { jd[0] = J2000; jd[1] = 0.0; for (n = 0; n < TEST_LENGTH; ++n) { jd[0] += 1.0; //#pragma omp critical //std::cout << "JD = " << jd[0] << "\n"; for (p = 0; p < 8; ++p) { planetpv(jd,p,pv); radecdist(pv,position); //#pragma omp critical //std::cout << p << " " << position[0] << " " << position[1] << "\n"; } } } clock_t endt=clock(); std::cout <<"Time: " << double(endt-start)/CLOCKS_PER_SEC * 1000 << " ms\n"; //doesn't work on vc++ without this cout << position[0] << jd[0] << pv[0][0]; return 0;}