module Zone = struct exception Hit of float array exception NoBone type fray = float array * float array type 'a rayp = float -> fray -> 'a type curlp = int * float array array -> float array array type 'a nAcc = {ray : 'a rayp; curl : curlp} let aini = Array.init and len = Array.length and ex f z = if f>z then z else z+1 let newz ds2 (neigh : ('a nAcc) array) id = let n = (len ds2) - 1 in let gd = zzm "gd" (aini n (fun i -> aini n (fun j -> 0.5 *. (ds2.(i).(n) +. ds2.(j).(n) -. ds2.(i).(j))))) in let aug a = (let n = len a in let lr = Array.make (n+1) 0. in aini (n+1) (fun i -> if i < n then aini (n+1) (fun j -> let e = (if j < n then a.(i).(j) else Array.fold_left (-.) 0. a.(i)) in lr.(j) <- lr.(j) -. e; e) else lr)) and dp u v = let s = ref 0. in (for i = 0 to n-1 do for j=0 to n-1 do s := !s +. u.(i) *. v.(j) *. gd.(i).(j) done done; !s) and pi a k = let lv = Array.fold_left (-.) 1. a in (aini (n+1) (fun j -> if j < k then (if j < n-1 then a.(j) else lv) else if j=k then 0. else if j < n then a.(j-1) else lv)) and po z k = (if (abs_float z.(k) > 1e-11) then raise (Hit z); aini n (fun j -> if j=n-1 then 0. else z.(ex k j))) and vol = let g = Linear.determ gd in assert (g>0.); sqrt g (* n! times real vol *) in let (gu, nfn) = let agui = aug (zzm "gu" (Linear.inv gd)) in (agui, zzm "nfn" (aini (n+1) (fun j -> Linear.sm vol agui.(j)))) in let cvf2z v fn = pzi "fn" fn; let a = zzf "va" v.(n-1) and w = zza "nn" nfn.(fn) and e j = if j=n-1 then (let q = ref 0. in for j = 0 to n-2 do q := !q -. v.(j) done; !q) else v.(j) in aini (n+1) (fun j -> a *. w.(j) +. if j if j = n-1 then a else a *. w.(ex k j))) in aini n (fun j -> (if j < n-1 then v.(ex k j) else 0.) -. v.(k) *. rv.(j))) in (* let tst v = for i = 0 to n do (let v2 = cvf2z v i in let v3 = cvz2f v2 i in pza "v" v; pza "v2" v2; pza "v3" v3) done in for i = 0 to n-1 do (let r = Array.make n 0. in r.(i) <- 1.; tst r) done; *) if false then (* special for jo44 testing and such. *) (let va=nfn.(2) and vb= nfn.(3) in id "pek"; pza "nfna" va; pza "nfnb" vb; pzf "dpab" (dp va vb); pzf "dpaa" (dp va va); pzf "dpbb" (dp vb vb)); fun zn -> {ray = (fun lam ((wh, vel) : fray) -> id "ray"; let whz = pi wh zn and vz = cvf2z vel zn in let lm = ref infinity and efn = ref zn in for j=0 to n do if j <> zn then (if ((zzf "vzj" vz.(j)) < 0.) && ((zzf "offs" (whz.(j) +. !lm *. vz.(j))) < 0.) then (lm := -. whz.(j) /. vz.(j); efn := j)) done; if !lm = infinity then print_string "leaky zone\n"; (neigh.(!efn).ray (lam +. !lm) ((po (aini (n+1) (fun j -> whz.(j) +. !lm *. vz.(j))) !efn), cvz2f vz !efn) : fray)); curl = (fun (fb, x) -> id "curl"; let zb = ex zn fb in if zb=zn then raise NoBone; neigh.(zb).curl ((if zn < zb then zn else zn-1), aini (len x) (fun k -> (cvz2f (cvf2z x.(k) zn) zb))))} end;;