type pt = {x : float; y : float; xd : float; yd : float};;
let m = 10 and n = 10 and ai = Array.init in
let a = Array.init 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 (-:=) x v = x := !x -. v
and pf = Printf.printf in
let t = ref 0. and dt = 0.005 in
let tp = ref (-. dt /. 2.) in while !t < 10.55 do
  if !t > !tp then t +:= 1.; pf "t = %11.7f\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 %11.7f %11.7f %11.7f %11.7f\n"
         i j  p.x p.y p.xd p.yd done done;
 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 -:= p0 *. ddy; a.(i).(j).yd +:= p0 *. ddx;
     a.(i).(j+1).xd -:= p0 *. dey; a.(i).(j+1).yd +:= p0 *. !dex;
     a.(i+1).(j).xd -:= p0 *. dy.(j); 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 -:= p1 *. dy.(j); a.(i).(j+1).yd +:= p1 *. dx.(j);
     a.(i+1).(j).xd +:= p1 *. !dey; a.(i+1).(j).yd +:= p1 *. !dex;
     a.(i+1).(j+1).xd +:= p1 *. ddy; a.(i+1).(j+1).yd -:= p1 *. ddx;
     a.(i).(j).x +:= dt *. a.(i).(j).xd; a.(i).(j).y +:= dt *. a.(i).(j).yd
   done;
   a.(i).(n-1).x += dt *. a.(i).(n-1).xd; a.(i).(n-1).y +:= dt *. a.(i).(n-1).yd
 done;
 for j=0 to n-1 do
   a.(m-1).(j).x +:= dt *. a.(m-1).(j).xd;
   a.(m-1).(j).y +:= dt *. a.(m-1).(j).yd done;
 t +:= dt
done