(* ReactiveML : interactive simulation with rmltop  *)
(* n-body simulation                                *)

type planet = 
    { id : int;
      mass : float;
      pos : float * float * float;
      speed : float * float * float; }
;;

(* Constants *)
let g = 6.67;;
let dt = 0.1;;

(* Global signal *)
signal env default [] gather (fun x y -> x :: y);;

(* --------------------------------------------------------------------- *)
(* Auxiliary functions *)
let random_speed () =
  ((Random.float 100.0) -. 50.0,
   (Random.float 100.0) -. 50.0,
   (Random.float 100.0) -. 50.0)
;;

let new_pos x y =
  let max_x_2 = (Graphics.size_x()) / 2 in
  let max_y_2 = (Graphics.size_y()) / 2 in
  (float_of_int (x - max_x_2),
   float_of_int (y - max_y_2), 
   (Random.float 200.0) -. 100.0)
;;


let random_pos () =
  let x = ((Random.int 200) - 100) + (Graphics.size_x()) / 2 in
  let y = ((Random.int 200) - 100) + (Graphics.size_y()) / 2 in
  new_pos x y
;;

let distance2 (x,y,z) (x',y',z') =
  (x' -. x)*.(x' -. x)
    +. (y' -. y)*.(y' -. y)
    +. (z' -. z)*.(z' -. z)
;;

let distance pos1 pos2 = sqrt (distance2 pos1 pos2)
;;

let new_planet =
  let cpt = ref 0 in
  fun pos ->
    incr cpt;
    { id = !cpt;
      mass = 1.0;
      pos = pos;
      speed = random_speed(); }
;;

let random_planet () = 
  new_planet (random_pos());;

let color_of_int = function
  | 0 -> Graphics.yellow
  | 1 -> Graphics.blue
  | 2 -> Graphics.green
  | 3 -> Graphics.red
  | 4 -> Graphics.cyan
  | 5 -> Graphics.black
  | 6 -> Graphics.magenta
  | _ -> Graphics.black
;;

(* --------------------------------------------------------------------- *)

(* Graphical window *)

let update_display all =
  Graphics.clear_graph();
  let max_x_2 = (Graphics.size_x()) / 2 in
  let max_y_2 = (Graphics.size_y()) / 2 in
  List.iter 
    (fun { id=id; pos=(x,y,z) } -> 
      Graphics.set_color (color_of_int (id mod 7));
      Graphics.fill_circle 
        (int_of_float x + max_x_2)
        (int_of_float y  + max_y_2)
        ( if (z < 5000.0) & (z > -250.0)
        then 5+(int_of_float z / 50) else 1))
    (List.sort (fun {pos = (_,_,z1)} {pos = (_,_,z2)} -> compare z1 z2) all)
;;

let process window =
  Graphics.open_graph "";
  Graphics.auto_synchronize false;
  loop
    await env (all) in 
    update_display all;
    Graphics.synchronize()
  end
;;

(* #run window;; *)


(* --------------------------------------------------------------------- *)
(* planet definition *)

let compute_pos =
  let force 
        { pos= (x1,y1,z1) as pos1; mass=m1 } 
        { pos= (x2,y2,z2) as pos2; mass=m2 } =
    let d2 = distance2 pos1 pos2 in
    let d = sqrt d2 in
    if (d <> 0.0) then 
      let  f12 = g *. (m1 *. m2) /. d2 in
      (f12 *. (x2 -. x1) /. d,
       f12 *. (y2 -. y1) /. d,
       f12 *. (z2 -. z1) /. d)
    else
      (0.0, 0.0, 0.0)
  in
  fun ({ pos=(x,y,z); speed=(x',y',z') } as me) all ->
    let fx, fy, fz = 
      (List.fold_left 
         (fun (fx,fy,fz) p -> 
           let x,y,z = force me p in
           (fx +. x), 
           (fy +. y), 
           (fz +. z)) 
         (0.0, 0.0, 0.0) 
         all)
    in
    let (sx, sy, sz) as speed = 
      (x' +. fx *. dt,
       y' +. fy *. dt,
       z' +. fz *. dt)
    in
    let pos = (x +. sx *. dt,
               y +. sy *. dt,
               z +. sz *. dt)
    in
    { id = me.id;
      mass = me.mass;
      pos = pos;
      speed = speed; }
;;

let process planet =
  let me = ref (random_planet()) in
  loop 
    emit env !me;
    await env (all) in
    me := compute_pos !me all
  end
;;

(* #run planet;; *)

let process sun =
  let me = 
    { id = 0; 
      mass =  30000.0; 
      pos = (0.0, 0.0, 0.0); 
      speed = (0.0, 0.0, 0.0) } 
  in
  loop 
    emit env me;
    pause
  end
;;

(* #run sun;; *)

(* #exec (for i = 1 to 50 dopar run planet done);; *)


let eclipse { pos = (x, y, z) } = 
  abs_float x < 5. && abs_float y < 5. && z > 0.;;

(* List.filter eclipse (pre ?env);; *)

let process eclipse_observer =
  loop
    await env (all) in
    if List.exists eclipse all then #suspend
  end
;;

(* #run eclipse_observer;; *)