Source

mutated_ocaml / testsuite / tests / misc-unsafe / almabench.ml

Full commit
(*
 * ALMABENCH 1.0.1
 *      OCaml     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