Source

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

Full commit
camlspotter d2116cf 

camlspotter 35cc436 
camlspotter d2116cf 














camlspotter 4eace9a 
camlspotter d2116cf 


















































camlspotter 4eace9a 
camlspotter d2116cf 

















































































camlspotter 4eace9a 
camlspotter d2116cf 




























camlspotter 4eace9a 
camlspotter d2116cf 




camlspotter 4eace9a 
camlspotter d2116cf 
camlspotter 4eace9a 
camlspotter d2116cf 













camlspotter 4eace9a 

camlspotter d2116cf 
camlspotter 4eace9a 
camlspotter d2116cf 






















camlspotter 4eace9a 
camlspotter d2116cf 



camlspotter 4eace9a 


camlspotter d2116cf 
camlspotter 4eace9a 
camlspotter d2116cf 

camlspotter 4eace9a 

camlspotter d2116cf 







camlspotter 4eace9a 
camlspotter d2116cf 










camlspotter 4eace9a 



camlspotter d2116cf 
camlspotter 4eace9a 



camlspotter d2116cf 
camlspotter 4eace9a 



camlspotter d2116cf 

camlspotter 4eace9a 
camlspotter d2116cf 



camlspotter 4eace9a 

camlspotter d2116cf 


camlspotter 4eace9a 
camlspotter d2116cf 

camlspotter 4eace9a 
camlspotter d2116cf 




camlspotter 4eace9a 
camlspotter d2116cf 















camlspotter 4eace9a 

camlspotter d2116cf 

(*
 * 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