(* Descended most directly from b.ml *)
type tv = float * float;;
type zn = {mutable p : tv; mutable v : tv};;
open Graphics;;
open_graph "";
set_window_title "Thud e";
ignore (wait_next_event [Key_pressed]);
auto_synchronize (1<0);
let (+^) (a, b) (c, d) = a +. c, b +. d
and (-^) (a, b) (c, d) = a -. c, b -. d
and (!^) (a, b) = (-. b , a)
and ( *^) a (b, c) = (a *. b, a *. c)
and ( *^.) (a, b) (c, d) = a *. c +. b *. d
and neg (a, b) = (-. a), (-. b)
and m = 20 and n = 20 and ai = Array.init and pf = Printf.printf in
let ci = int_of_float 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 -> let u, v = (sc *^ r) in (50+(ci u), 50+(ci v))) in
let dl p pe = let x, y = i p
   and xe, ye = i pe in moveto x y; lineto xe ye in
let av r v = r.v <- r.v +^ v and ap r v = r.p <- r.p +^ v
and a = ai m (fun i -> ai n (fun j ->
  {p = (float i) +. (float j) *. 0.5, ((sqrt 3.) /. 2.) *. (float j);
   v = 0., 0.}))
and dt = 0.005 and vl = max 0.0001 in
let ct v p q r =
  set_color (rgb (ci (280. *. v)) 0 0);
  if 1<2 then fill_poly [| i p; i q; i r |];
  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 (pf "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 () in
let pp (a, b) = pf "%15.11f %15.11f " a b in
let prn t = (pf "t = %10.6f\n" t;
    for i=0 to m-1 do for j=0 to n-1 do
      pf "%2d %2d " i j; pp a.(i).(j).p; pp a.(i).(j).v;
      pf "\n"; done done) in
let rec cycle t tp =
  if t > tp then (prn t; cycle t (tp +. 1.))
  else (
  clear_graph ();
  let dp = ai (n-1) (fun j -> a.(0).(j+1).p -^ a.(0).(j).p) in
  for i=0 to m-2 do
    let de = ref (!^ (a.(i+1).(0).p -^ a.(i).(0).p)) in
    for j=0 to n-2 do
      let dd = !^ (a.(i).(j+1).p -^ a.(i+1).(j).p) in
      let p0 = (let vol = ((dp.(j)) *^. !de) in
        ct vol a.(i).(j).p a.(i).(j+1).p a.(i+1).(j).p;
        dt *. (1./. (vl vol) -. 2.)) in
      av a.(i).(j) (p0 *^ dd);
      let r, s = p0 *^ !de, p0 *^ !^ (dp.(j)) in
      dp.(j) <- a.(i+1).(j+1).p -^ a.(i+1).(j).p;
      de := !^ (a.(i+1).(j+1).p -^ a.(i).(j+1).p);
      let p1 = let vol = ((dp.(j)) *^. !de) in
        ct vol a.(i+1).(j+1).p a.(i).(j+1).p a.(i+1).(j).p;
        dt *. (1./. (vl vol) -. 2.) in
      av a.(i).(j+1) ((p1 *^ !^ (dp.(j))) +^ r);
      av a.(i+1).(j) (neg (p1 *^ !de +^ s));
      av a.(i+1).(j+1) (neg (p1 *^ dd));
      ap a.(i).(j) (dt *^ a.(i).(j).v);
      dl a.(i).(j).p a.(i+1).(j).p;
      dl a.(i).(j).p a.(i).(j+1).p;
      dl a.(i+1).(j).p a.(i).(j+1).p
    done;
    ap a.(i).(n-1) (dt *^ a.(i).(n-1).v);
    dl a.(i).(n-1).p a.(i+1).(n-1).p
  done;
  for j=0 to n-1 do ap a.(m-1).(j) (dt *^ a.(m-1).(j).v) done;
  for j=0 to n-2 do dl a.(m-1).(j).p a.(m-1).(j+1).p done;
  moveto 2 2; draw_string (Printf.sprintf "%10.4f" t);
  synchronize (); ignore (tr ());
  cycle (t +. dt) tp)
in cycle 0. (-. dt /. 2.)