(* 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.);