module Jig = (struct
open Zone open Linear
exception Escape of Zone.fray
type cray = float array * float array

let al = Array.length and aini = Array.init

let mix (top, embd, sbl) perm = let pll t = List.map
   (fun zn -> List.map (fun v -> perm.(v)) zn) t in pll top,
   (let n = (al embd) in let e = Array.make n embd.(0) in
     for j=0 to n-1 do e.(perm.(j)) <- embd.(j) done; e),
   List.map (sort (<)) (pll sbl)

let nz (embd, sg) vl =
  let v = Array.of_list vl in let np = (al v) in
  let cm (i, j) = (let n = al embd.(0) and ds2 = ref 0. in
  for k = 0 to n-1 do let e = embd.(v.(i)).(k) -. embd.(v.(j)).(k) in
  ds2 := !ds2 +. sg.(k) *. e *. e done; !ds2) in
  let neigh = (aini np (fun j ->
    {ray = (fun lam dzray -> raise (Escape dzray)); 
     curl = (fun (_, x) -> pzi "Leak!" j; [|[|42.|]|])})) in
  let ds2 = aini np (fun i -> aini np (fun j ->
       cm (i, j))) in
   ((fun fn ex -> neigh.(fn) <- ex),
     newz ds2 neigh (fun m -> pzl (m^" zone entry ") vl))

let prep embd fvl o =
   (* let l2 v = let s = ref 0. in for i = 0 to ((al v) - 1)
      do s := !s +. v.(i) *. v.(i) done; !s in
   let norm v = sm (1. /. (sqrt (l2 v))) v in *)
   let fva = Array.of_list fvl in
     let n = al fva and e j = embd.(fva.(j)) in
     let no = (vneg (e (n-1))) in let cto = trans (let z =
       (aini (n-1) (fun j -> vadd (e j) no)) in
       aini n (fun k -> if k<n-1 then z.(k) else
         (let n = cross z in if o then n else vneg n))) in
      (n, e, cto)

type tooll = {lray : cray -> cray; lcurl : curlp}
let jig embd bl inf outf =
let (io, (_, launch)) = List.assoc inf bl
and (oo, (intro, _)) = List.assoc outf bl
and p fvl o = (match prep embd fvl o with (n, e, cto) -> (let cti = inv cto in
    ((fun (pos, dir) ->
      ((mtv cti (vadd (vneg (zza "foy" (e (n-1)))) pos)), (mtv cti dir))),
     (fun (pos, dir) -> ((vadd (mtv cto pos) (e (n-1))),
       (mtv cto dir))), cti, cto))) in
   let (_, out, _, ctoX) = p outf oo in
   intro {ray = (fun lam (wh, dir) ->
     let n = (al dir) - 1 in dir.(n) <- -. dir.(n);
     out (wh, dir));
     curl = (fun (j, x) ->
        print_string "out"; print_newline (); print_int j; x)};
   let (inT, _, cti, _) = p inf io in
{lray = (fun ingress -> (launch.ray 0. (inT ingress)));
  lcurl = (fun (j, a) ->
   Array.map (mtv (let na = Array.copy ctoX and n = (al ctoX)-1 in
   for j = 0 to n do na.(j) <- Array.copy ctoX.(j);
   (* This code is grotesque!! *)
   na.(j).(n) <- -. na.(j).(n) done; na))
   (launch.curl (j, Array.map (mtv cti) a)))}

let pyt = [|1.; 1.; 1.; 1.|]
end : sig
type cray = float array * float array
open Zone
val mix :
  (int list list * 'a array * int list list) ->
  int array -> int list list * 'a array * int list list
val nz :
  float array array (* ds2 *) * float array (* signature *) ->
           int list (* global vertex list for new zone *) ->
           (int -> cray nAcc -> unit) *
           (int -> cray nAcc)
type tooll = {lray : cray -> cray; lcurl : curlp}
val jig :
  float array array ->
  (int list *
   (bool *
    ((cray nAcc -> unit) * cray nAcc)))
  list -> int list -> int list -> tooll
val pyt : float array
end);;