(* 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";
ignore (wait_next_event [Key_pressed]);
auto_synchronize (1<0);
let m = 20 and n = 15 and ai = Array.init
and pf = Printf.printf in
let sc = let xl, yl = float (size_x ()), float (size_y ()) in
   0.9 *. (min (xl /. ((float m) +. 0.5 *. (float n)))
   (yl /. (0.866 *. (float n)))) in
let dr x y xe ye = let i r = 50 + (int_of_float (sc *. r)) in
  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 in
let 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 (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; !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 %13.9f %13.9f %13.9f %13.9f\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 p0 = dt *. (1. /. (vl (dy.(j) *. !dex -. dx.(j) *. !dey)) -. 2.)
     and ddx = a.(i).(j+1).x -. a.(i+1).(j).x
     and ddy = a.(i).(j+1).y -. a.(i+1).(j).y 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 =  dt *. (1./. (vl (dy.(j) *. !dex -. dx.(j) *. !dey)) -. 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;
     dr a.(i).(j).x a.(i).(j).y a.(i+1).(j).x a.(i+1).(j).y;
     dr a.(i).(j).x a.(i).(j).y a.(i).(j+1).x a.(i).(j+1).y;
     dr 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;
   dr 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
   dr a.(m-1).(j).x a.(m-1).(j).y a.(m-1).(j+1).x a.(m-1).(j+1).y done;
 let f = 0.995 in if f<1. then for i=0 to m-1 do for j=0 to n-1 do
   if i=0 or i=m-1 or j=0 or j=n-1 then (
   a.(i).(j).xd <- f *. a.(i).(j).xd;
   a.(i).(j).yd <- f *. a.(i).(j).yd) done done;
 t +:= dt;
 moveto 2 2; draw_string (Printf.sprintf "%10.4f" !t);
 synchronize ();
 ignore (tr ()) done;
prn (); ignore (tr ())