(*
 * 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 if a < 0.0 then w +. twopi else w -. twopi 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;
  (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);
  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 || 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 _ = 0 to test_loops - 1 do
    jd.(0) <- j2000;
    jd.(1) <- 0.0;
    for _ = 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