(* Movie *)
(* See "wait_timed_write" in view-source:http://caml.inria.fr/pub/docs/manual-ocaml/libref/Thread.html#VALwait_timed_read *)
type pt = {mutable x  : float; mutable y  : float;
           mutable xd : float; mutable yd : float};;
open Graphics;;
open_graph "";
set_window_title "Thud e";
ignore (wait_next_event [Key_pressed]);
auto_synchronize (1<0);
let m = 20 and n = 20 and ai = Array.init
and pf = Printf.printf and qut = ref 1 and pls = ref [] in
let i = let sc =
   0.9 *. (min ((float (size_x ())) /. ((float m) +. 0.5 *. (float n)))
   ((float (size_y ())) /. (0.866 *. (float n)))) in
   (fun r -> 50 + (int_of_float (sc *. r))) in
let dl x y xe ye = moveto (i x) (i y); lineto (i xe) (i ye) in
let a = ai m (fun i -> ai n (fun j ->
  {x = (float i) +. (float j) *. 0.5; y = ((sqrt 3.) /. 2.) *. (float j);
     xd = 0.; yd = 0.})) in
let (+:=) x v = x := !x +. v
and ct v p q r =
  if !qut = 3 then pls := v::!pls;
  set_color (rgb (int_of_float (280. *. v)) 0 0);
  if 1<2 then fill_poly [| i p.x, i p.y; i q.x, i q.y; i r.x, i r.y |];
  set_color black
and tarry _ = let c = ref 0 and d = ref 1 in let rec t cnt =
  let cd = Char.code and c = read_key () in
    if '0' <= c && c <= '9' then t (10*cnt + (cd c) - (cd '0'))
    else if c = '\r' then cnt
    else if c = 'q' then (qut := 2; cnt)
    else (Printf.printf "BadChar %d >%s<\n" (cd c) (Char.escaped c); cnt)
    in  (fun _ -> if !c>0 then (c := pred !c; !c) else
    (c := (let q = t 0 in if q > 0 then d := q; !d) -1;
      flush stdout; !c)) in
let tr = tarry ()
and t = ref 0. and dt = 0.005 in
let prn _ = (pf "t = %13.9f\n" !t;
   for i=0 to m-1 do for j=0 to n-1 do let p = a.(i).(j) in
      pf "%2d %2d %16.12f %16.12f %16.12f %16.12f\n"
         i j  p.x p.y p.xd p.yd done done; flush stdout)
and tp = ref (-. dt /. 2.) in while 0<1 do
 clear_graph ();
 if !t > !tp then (tp +:= 1.; prn ());
 let vl = max 0.0001
 and dx = ai (n-1) (fun j -> a.(0).(j+1).x -. a.(0).(j).x)
 and dy = ai (n-1) (fun j -> a.(0).(j+1).y -. a.(0).(j).y) in
 for i=0 to m-2 do
   let dex = ref (a.(i+1).(0).x -. a.(i).(0).x)
   and dey = ref (a.(i+1).(0).y -. a.(i).(0).y) in
   for j=0 to n-2 do
     let ddx = a.(i).(j+1).x -. a.(i+1).(j).x
     and ddy = a.(i).(j+1).y -. a.(i+1).(j).y in
    (let p0 = (let vol = (dy.(j) *. !dex -. dx.(j) *. !dey) in
       ct vol a.(i).(j) a.(i).(j+1) a.(i+1).(j);
       dt *. (1./. (vl vol) -. 2.)) in
     a.(i).(j).xd <- a.(i).(j).xd -. p0 *. ddy;
     a.(i).(j).yd <- a.(i).(j).yd +. p0 *. ddx;
     a.(i).(j+1).xd <- a.(i).(j+1).xd -. p0 *. !dey;
     a.(i).(j+1).yd <- a.(i).(j+1).yd +. p0 *. !dex;
     a.(i+1).(j).xd <- a.(i+1).(j).xd +. p0 *. dy.(j);
     a.(i+1).(j).yd <- a.(i+1).(j).yd -. p0 *. dx.(j);
     dx.(j) <- a.(i+1).(j+1).x -. a.(i+1).(j).x;
     dy.(j) <- a.(i+1).(j+1).y -. a.(i+1).(j).y;
     dex := a.(i+1).(j+1).x -. a.(i).(j+1).x;
     dey := a.(i+1).(j+1).y -. a.(i).(j+1).y);
     let p1 = let vol = (dy.(j) *. !dex -. dx.(j) *. !dey) in
       ct vol a.(i+1).(j+1) a.(i).(j+1) a.(i+1).(j);
       dt *. (1./. (vl vol) -. 2.) in
     a.(i).(j+1).xd <- a.(i).(j+1).xd -. p1 *. dy.(j);
     a.(i).(j+1).yd <- a.(i).(j+1).yd +. p1 *. dx.(j);
     a.(i+1).(j).xd <- a.(i+1).(j).xd +. p1 *. !dey;
     a.(i+1).(j).yd <- a.(i+1).(j).yd -. p1 *. !dex;
     a.(i+1).(j+1).xd <- a.(i+1).(j+1).xd +. p1 *. ddy;
     a.(i+1).(j+1).yd <- a.(i+1).(j+1).yd -. p1 *. ddx;
     a.(i).(j).x <- a.(i).(j).x +. dt *. a.(i).(j).xd;
     a.(i).(j).y <- a.(i).(j).y +. dt *. a.(i).(j).yd;
     dl a.(i).(j).x a.(i).(j).y a.(i+1).(j).x a.(i+1).(j).y;
     dl a.(i).(j).x a.(i).(j).y a.(i).(j+1).x a.(i).(j+1).y;
     dl a.(i+1).(j).x a.(i+1).(j).y a.(i).(j+1).x a.(i).(j+1).y
   done;
   a.(i).(n-1).x <- a.(i).(n-1).x +. dt *. a.(i).(n-1).xd;
   a.(i).(n-1).y <- a.(i).(n-1).y +. dt *. a.(i).(n-1).yd;
   dl a.(i).(n-1).x a.(i).(n-1).y a.(i+1).(n-1).x a.(i+1).(n-1).y
 done;
 for j=0 to n-1 do
   a.(m-1).(j).x <- a.(m-1).(j).x +. dt *. a.(m-1).(j).xd;
   a.(m-1).(j).y <- a.(m-1).(j).y +. dt *. a.(m-1).(j).yd done;
 for j=0 to n-2 do
   dl a.(m-1).(j).x a.(m-1).(j).y a.(m-1).(j+1).x a.(m-1).(j+1).y done;
 if 1<0 then for i=0 to m-1 do for j=0 to n-1 do
   a.(i).(j).xd <- 0.9998 *. a.(i).(j).xd;
   a.(i).(j).yd <- 0.9998 *. a.(i).(j).yd done done;
 t +:= dt;
 moveto 2 2; draw_string (Printf.sprintf "%10.4f" !t);
 synchronize (); ignore (tr ());
 if !qut = 3 then (
   let z = Array.of_list (List.sort (fun x y -> if x<y then -1 else 1) !pls) in
   pls := []; let n = (Array.length z) / 20 in
   for j=0 to n-1 do pf "%3d %9.6f\n" j z.(20*j) done;
pf "Zimm %d\n" n; flush stdout; qut := 1)
   else if !qut = 2 then qut := 3
 done;
prn (); ignore (tr ())