(* Descended most directly from b.ml *)
type tv = float * float;;
type zn = {mutable p : tv; mutable v : tv};;
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 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 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 (
  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 = dt *. (1. /. (vl ((dp.(j)) *^. !de)) -. 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 = dt *. (1. /. (vl ((dp.(j)) *^. !de)) -. 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)
    done;
    ap a.(i).(n-1) (dt *^ a.(i).(n-1).v)
  done;
  for j=0 to n-1 do ap a.(m-1).(j) (dt *^ a.(m-1).(j).v) done;
  cycle (t +. dt) tp)
in cycle 0. (-. dt /. 2.);