diff options
Diffstat (limited to 'test/cminor/almabench.cmp')
| -rw-r--r-- | test/cminor/almabench.cmp | 159 |
1 files changed, 159 insertions, 0 deletions
diff --git a/test/cminor/almabench.cmp b/test/cminor/almabench.cmp new file mode 100644 index 0000000..e9e8392 --- /dev/null +++ b/test/cminor/almabench.cmp @@ -0,0 +1,159 @@ +#define PI 3.14159265358979323846 +#define J2000 2451545.0 +#define JCENTURY 36525.0 +#define JMILLENIA 365250.0 +#define TWOPI (2.0 *f PI) +#define A2R (PI /f 648000.0) +#define R2H (12.0 /f PI) +#define R2D (180.0 /f PI) +#define GAUSSK 0.01720209895 +#define TEST_LOOPS 20 +#define TEST_LENGTH 36525 +#define sineps 0.3977771559319137 +#define coseps 0.9174820620691818 + +/* Access to tables */ +#define amas(n) float64["amas" + (n) * 8] +#define a(x,y) float64["a" + ((x) * 24 + (y) * 8)] +#define dlm(x,y) float64["dlm" + ((x) * 24 + (y) * 8)] +#define e(x,y) float64["e" + ((x) * 24 + (y) * 8)] +#define pi(x,y) float64["pi" + ((x) * 24 + (y) * 8)] +#define dinc(x,y) float64["dinc" + ((x) * 24 + (y) * 8)] +#define omega(x,y) float64["omega" + ((x) * 24 + (y) * 8)] +#define kp(x,y) float64["kp" + ((x) * 72 + (y) * 8)] +#define ca(x,y) float64["ca" + ((x) * 72 + (y) * 8)] +#define sa(x,y) float64["sa" + ((x) * 72 + (y) * 8)] +#define kq(x,y) float64["kq" + ((x) * 80 + (y) * 8)] +#define cl(x,y) float64["cl" + ((x) * 80 + (y) * 8)] +#define sl(x,y) float64["sl" + ((x) * 80 + (y) * 8)] + +/* Function calls */ +#define cos(x) ("cos_static"(x): float -> float) +#define sin(x) ("sin_static"(x): float -> float) +#define atan2(x,y) ("atan2_static"(x,y): float -> float -> float) +#define asin(x) ("asin_static"(x): float -> float) +#define sqrt(x) ("sqrt_static"(x): float -> float) +#define fmod(x,y) ("fmod_static"(x,y): float -> float -> float) +#define anpm(x) ("anpm"(x) : float -> float) + +"anpm"(a): float -> float +{ + var w; + w = fmod(a,TWOPI); + if (absf(w) >=f PI) + w = w -f ((a <f 0.0) ? -f TWOPI : TWOPI); + return w; +} + +"planetpv" (epoch_, np, pv_): int -> int -> int -> void +{ +#define epoch(x) float64[epoch_ + (x) * 8] +#define pv(x,y) float64[pv_ + (x) * 24 + (y) * 8] + + var i, j, k; + var t, da, dl, de, dp, di; + var doh, dmu, arga, argl, am; + var ae, dae, ae2, at, r, v; + var si2, xq, xp, tl, xsw; + var xcw, xm2, xf, ci2, xms, xmc; + var xpxq2, x, y, z; + + t = ((epoch(0) -f J2000) +f epoch(1)) /f JMILLENIA; + + da = a(np,0) +f (a(np,1) +f a(np,2) *f t ) *f t; + dl = (3600.0 *f dlm(np,0) +f (dlm(np,1) +f dlm(np,2) *f t ) *f t ) *f A2R; + de = e(np,0) +f (e(np,1) +f e(np,2) *f t ) *f t; + dp = anpm((3600.0 *f pi(np,0) +f (pi(np,1) +f pi(np,2) *f t ) *f t ) *f A2R ); + di = (3600.0 *f dinc(np,0) +f (dinc(np,1) +f dinc(np,2) *f t ) *f t ) *f A2R; + doh = anpm((3600.0 *f omega(np,0) +f (omega(np,1) +f omega(np,2) *f t ) *f t ) *f A2R ); + + dmu = 0.35953620 *f t; + + k = 0; + {{ loop { + if (! (k < 8)) exit; + arga = kp(np,k) *f dmu; + argl = kq(np,k) *f dmu; + da = da +f (ca(np,k) *f cos(arga) +f sa(np,k) *f sin(arga)) *f 0.0000001; + dl = dl +f (cl(np,k) *f cos(argl) +f sl(np,k) *f sin(argl)) *f 0.0000001; + k = k + 1; + } }} + + arga = kp(np,8) *f dmu; + da = da +f t *f (ca(np,8) *f cos(arga) +f sa(np,8) *f sin(arga)) *f 0.0000001; + + k = 8; + {{ loop { + if (! (k <= 9)) exit; + argl = kq(np,k) *f dmu; + dl = dl +f t *f ( cl(np,k) *f cos(argl) +f sl(np,k) *f sin(argl) ) *f 0.0000001; + k = k + 1; + } }} + + dl = "fmod_static"(dl,TWOPI) : float -> float -> float; + + am = dl -f dp; + ae = am +f de *f sin(am); + k = 0; + + {{ loop { + dae = (am -f ae +f de *f sin(ae)) /f (1.0 -f de *f cos(ae)); + ae = ae +f dae; + k = k + 1; + + if ((k >= 10) || (absf(dae) <f 1e-12)) + exit; + } }} + + ae2 = ae /f 2.0; + at = 2.0 *f atan2(sqrt((1.0 +f de) /f (1.0 -f de)) *f sin(ae2), cos(ae2)); + + r = da *f (1.0 -f de *f cos(ae)); + v = GAUSSK *f sqrt((1.0 +f 1.0 /f amas(np) ) /f (da *f da *f da)); + + si2 = sin(di /f 2.0); + xq = si2 *f cos(doh); + xp = si2 *f sin(doh); + tl = at +f dp; + xsw = sin(tl); + xcw = cos(tl); + xm2 = 2.0 *f (xp *f xcw -f xq *f xsw ); + xf = da /f sqrt(1.0 -f de *f de); + ci2 = cos(di /f 2.0); + xms = (de *f sin(dp) +f xsw) *f xf; + xmc = (de *f cos(dp) +f xcw) *f xf; + xpxq2 = 2.0 *f xp *f xq; + + x = r *f (xcw -f xm2 *f xp); + y = r *f (xsw +f xm2 *f xq); + z = r *f (-f xm2 *f ci2); + + pv(0,0) = x; + pv(0,1) = y *f coseps -f z *f sineps; + pv(0,2) = y *f sineps +f z *f coseps; + + x = v *f ((-f 1.0 +f 2.0 *f xp *f xp) *f xms +f xpxq2 *f xmc); + y = v *f (( 1.0 -f 2.0 *f xq *f xq ) *f xmc -f xpxq2 *f xms); + z = v *f (2.0 *f ci2 *f (xp *f xms +f xq *f xmc)); + + pv(1,0) = x; + pv(1,1) = y *f coseps -f z *f sineps; + pv(1,2) = y *f sineps +f z *f coseps; + +#undef epoch +#undef pv +} + +"radecdist"(state_, rdd_): int -> int -> void +{ +#define state(x,y) float64[state_ + (x) * 24 + (y) * 8] +#define rdd(x) float64[rdd_ + (x) * 8] + + rdd(2) = sqrt(state(0,0) *f state(0,0) +f state(0,1) *f state(0,1) +f state(0,2) *f state(0,2)); + rdd(0) = atan2(state(0,1), state(0,0)) *f R2H; + if (rdd(0) <f 0.0) rdd(0) = rdd(0) +f 24.0; + rdd(1) = asin(state(0,2) /f rdd(2)) *f R2D; + +#undef state +#undef rdd +} |