type pt = {mutable x  : float; mutable y  : float;
           mutable xd : float; mutable yd : float};;
let m = 20 and n = 20 and ai = Array.init 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 pf = Printf.printf in
let t = ref 0. and dt = 0.005 in
let prn _ = (pf "t = %20.16f\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 %15.11f %15.11f %15.11f %15.11f\n"
         i j  p.x p.y p.xd p.yd done done)
and tp = ref (-. dt /. 2.) in while !t < 10.55 do
  if !t > !tp then (tp +:= 1.; prn ());
 let 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. /. (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./. (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
   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
 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;
 t +:= dt
done; prn ()