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