diff options
Diffstat (limited to 'test/c/almabench.c')
| -rw-r--r-- | test/c/almabench.c | 190 |
1 files changed, 190 insertions, 0 deletions
diff --git a/test/c/almabench.c b/test/c/almabench.c new file mode 100644 index 0000000..ef92a59 --- /dev/null +++ b/test/c/almabench.c @@ -0,0 +1,190 @@ +/* + almabench + Standard C version + 17 April 2003 + + 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) + + 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 <math.h> +#include <string.h> + +#define PI 3.14159265358979323846 +#define J2000 2451545.0 +#define JCENTURY 36525.0 +#define JMILLENIA 365250.0 +#define TWOPI (2.0 * PI) +#define A2R (PI / 648000.0) +#define R2H (12.0 / PI) +#define R2D (180.0 / PI) +#define GAUSSK 0.01720209895 +#define TEST_LOOPS 20 +#define TEST_LENGTH 36525 +#define sineps 0.3977771559319137 +#define coseps 0.9174820620691818 + +extern double amas[8]; +extern double a[8][3], dlm[8][3], e[8][3], pi[8][3], dinc[8][3], omega[8][3]; +extern double kp[8][9], ca[8][9], sa[8][9], kq[8][10], cl[8][10], sl[8][10]; + +extern double cos_static(double), sin_static(double), + atan2_static(double, double), asin_static(double), sqrt_static(double), + fmod_static(double, double), fabs_static(double); + +#define cos cos_static +#define sin sin_static +#define atan2 atan2_static +#define asin asin_static +#define sqrt sqrt_static +#define fmod fmod_static +#define fabs fabs_static + +//--------------------------------------------------------------------------- +// 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; +} |