(* * ALMABENCH 1.0.1 * Objective Caml version * * A number-crunching benchmark designed for cross-language and vendor * comparisons. * * Written by Shawn Wagner, from Scott Robert Ladd's versions for * C++ and java. * * 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. *) let pic = 3.14159265358979323846 and j2000 = 2451545.0 and jcentury = 36525.0 and jmillenia = 365250.0 let twopi = 2.0 *. pic and a2r = pic /. 648000.0 and r2h = 12.0 /. pic and r2d = 180.0 /. pic and gaussk = 0.01720209895 (* number of days to include in test *) let test_loops = 5 (* was: 20 *) and test_length = 36525 (* sin and cos of j2000 mean obliquity (iau 1976) *) and sineps = 0.3977771559319137 and coseps = 0.9174820620691818 and amas = [| 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) *) and a = [| [| 0.3870983098; 0.0; 0.0 |]; [| 0.7233298200; 0.0; 0.0 |]; [| 1.0000010178; 0.0; 0.0 |]; [| 1.5236793419; 3e-10; 0.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 |] |] and dlm = [| [| 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 |] |] and e = [| [| 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.0 |] |] and pi = [| [| 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 |] |] and dinc = [| [| 7.00498625; -214.25629; 0.28977 |]; [| 3.39466189; -30.84437; -11.67836 |]; [| 0.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 |] |] and omega = [| [| 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. *) and kp = [| [| 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 |] |] and ca = [| [| 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 |] |] and sa = [| [| -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 . *) and kq = [| [| 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 |] |] and cl = [| [| 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 |] |] and sl = [| [| -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. *) let anpm a = let w = mod_float a twopi in if abs_float w >= pic then begin if a < 0.0 then w +. twopi else w -. twopi end else w (* The reference frame is equatorial and is with respect to the * mean equator and equinox of epoch j2000. *) let planetpv epoch np pv = (* time: julian millennia since j2000. *) let t = ((epoch.(0) -. j2000) +. epoch.(1)) /. jmillenia in (* compute the mean elements. *) let da = ref (a.(np).(0) +. (a.(np).(1) +. a.(np).(2) *. t ) *. t) and dl = ref ((3600.0 *. dlm.(np).(0) +. (dlm.(np).(1) +. dlm.(np).(2) *. t ) *. t) *. a2r) and de = e.(np).(0) +. (e.(np).(1) +. e.(np).(2) *. t ) *. t and dp = anpm ((3600.0 *. pi.(np).(0) +. (pi.(np).(1) +. pi.(np).(2) *. t ) *. t ) *. a2r ) and di = (3600.0 *. dinc.(np).(0) +. (dinc.(np).(1) +. dinc.(np).(2) *. t ) *. t ) *. a2r and doh = anpm ((3600.0 *. omega.(np).(0) +. (omega.(np).(1) +. omega.(np).(2) *. t ) *. t ) *. a2r ) (* apply the trigonometric terms. *) and dmu = 0.35953620 *. t in (* loop invariant *) let kp = kp.(np) and kq = kq.(np) and ca = ca.(np) and sa = sa.(np) and cl = cl.(np) and sl = sl.(np) in for k = 0 to 7 do let arga = kp.(k) *. dmu and argl = kq.(k) *. dmu in da := !da +. (ca.(k) *. cos arga +. sa.(k) *. sin arga) *. 0.0000001; dl := !dl +. (cl.(k) *. cos argl +. sl.(k) *. sin argl) *. 0.0000001 done; begin let arga = kp.(8) *. dmu in da := !da +. t *. (ca.(8) *. cos arga +. sa.(8) *. sin arga ) *. 0.0000001; for k = 8 to 9 do let argl = kq.(k) *. dmu in dl := !dl +. t *. ( cl.(k) *. cos argl +. sl.(k) *. sin argl ) *. 0.0000001 done; end; dl := mod_float !dl twopi; (* iterative solution of kepler's equation to get eccentric anomaly. *) let am = !dl -. dp in let ae = ref (am +. de *. sin am) and k = ref 0 in let dae = ref ((am -. !ae +. de *. sin !ae) /. (1.0 -. de *. cos !ae)) in ae := !ae +. !dae; incr k; while !k < 10 or abs_float !dae >= 1e-12 do dae := (am -. !ae +. de *. sin !ae) /. (1.0 -. de *. cos !ae); ae := !ae +. !dae; incr k done; (* true anomaly. *) let ae2 = !ae /. 2.0 in let at = 2.0 *. atan2 (sqrt ((1.0 +. de) /. (1.0 -. de)) *. sin ae2) (cos ae2) (* distance (au) and speed (radians per day). *) and r = !da *. (1.0 -. de *. cos !ae) and v = gaussk *. sqrt ((1.0 +. 1.0 /. amas.(np) ) /. (!da *. !da *. !da)) and si2 = sin (di /. 2.0) in let xq = si2 *. cos doh and xp = si2 *. sin doh and tl = at +. dp in let xsw = sin tl and xcw = cos tl in let xm2 = 2.0 *. (xp *. xcw -. xq *. xsw ) and xf = !da /. sqrt (1.0 -. de *. de) and ci2 = cos (di /. 2.0) in let xms = (de *. sin dp +. xsw) *. xf and xmc = (de *. cos dp +. xcw) *. xf and xpxq2 = 2.0 *. xp *. xq in (* position (j2000 ecliptic x,y,z in au). *) let x = r *. (xcw -. xm2 *. xp) and y = r *. (xsw +. xm2 *. xq) and z = r *. (-.xm2 *. ci2) in (* 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). *) let x = v *. ((-1.0 +. 2.0 *. xp *. xp) *. xms +. xpxq2 *. xmc) and y = v *. (( 1.0 -. 2.0 *. xq *. xq ) *. xmc -. xpxq2 *. xms) and z = v *. (2.0 *. ci2 *. (xp *. xms +. xq *. xmc)) in (* 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. *) let radecdist state rdd = (* 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 then rdd.(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 *) let _ = let jd = [| 0.0; 0.0 |] and pv = [| [| 0.0; 0.0; 0.0 |]; [| 0.0; 0.0; 0.0 |] |] and position = [| 0.0; 0.0; 0.0 |] in (* Test *) jd.(0) <- j2000; jd.(1) <- 1.0; for p = 0 to 7 do planetpv jd p pv; radecdist pv position; Printf.printf "%d %.2f %.2f\n%!" p position.(0) position.(1) done; (* Benchmark *) for i = 0 to test_loops - 1 do jd.(0) <- j2000; jd.(1) <- 0.0; for n = 0 to test_length - 1 do jd.(0) <- jd.(0) +. 1.0; for p = 0 to 7 do planetpv jd p pv; radecdist pv position; done done done