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<fn then e j else if j=fn then 0. else e (j-1))
      and cvz2f v k = (let w = nfn.(k) in let a = 1. /. w.(k) in
        (* array rv below expresses zone basis element k in terms of the bcc of the facet.
            Note that this does not depend on v and might be computed in morphogenesis. 
            Alternatively its creation might be made unnecessary at its point of use. *)
        let rv = (aini n (fun 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;;