open Jog open Zone;;
let coo = [|[|0.; 0.; -0.; 1.|]; [|0.; 0.; 1.; -0.|];
  [| 0.; 1.; 0.; 0.|]; [|1.; 0.; 0.; 0.|]; [|0.; -1.; 0.; 0.|]; [|-1.; 0.; 0.; 0.|];
  [|0.; 0.; 0.; 0.|]|]
and sq w = w *. w in
let pyt i j = let s = ref 0. in for k=0 to (len coo.(0)) - 1 do s := !s +.
     (sq (coo.(i).(k) -. coo.(j).(k))) done; !s
and noise = mem (len coo) in
let top = [[6; 0; 1; 2; 3]; [6; 0; 1; 3; 4]; [6; 0; 1; 4; 5]; [6; 0; 1; 5; 2]]
and ini = ref (fun _ -> raise (OffTrack 3))
and ds2 i j = (pyt i j) *. (0.9 +. (0.5 *. (noise i j)))
and fac = [0; 1; 2; 6] and bid = 2 in
let fa = Array.of_list fac in
let n = len fa in
let gdp = aini (n-1) (fun i -> aini (n-1) (fun j -> 0.5 *. (
  (ds2 fa.(n-1) fa.(i)) +. (ds2 fa.(n-1) fa.(j)) -. (ds2 fa.(i) fa.(j))))) in
let gd = aini n (fun i -> aini n (fun j ->
  if i<n-1 && j<n-1 then gdp.(i).(j)
  else if i=n-1 && j=n-1 then Linear.determ gdp
  else 0.)) in
let ip a b = let s = ref 0. in for i=0 to n-1 do for j=0 to n-1 do
   s := !s +. a.(i) *. b.(j) *. gd.(i).(j) done done; !s in
let pON = Linear.gs ip (let fON = Linear.iden n in
     (let t = fON.(bid) in fON.(bid) <- fON.(n-2); fON.(n-2) <- t); fON) in
let _ = Glue.morphgen (nz ds2) top
    [(fac, ((fun  x -> x),
      (fun t -> ini := t.curl; {t with curl = (fun (k, ball) -> ball)})))] in
    (pON, Linear.mul  (zzm "u" (!ini (bid, pON))) (Linear.inv pON))

(* should yield (_, ((1, 0, 0, 0), (0, 1, 0, 0), (0, 0, cos r, sin r), (0, 0, -sin r, cos r))) *);;