diff --git a/Basis/Basis.ml b/Basis/Basis.ml
index 91f8f64..75c7561 100644
--- a/Basis/Basis.ml
+++ b/Basis/Basis.ml
@@ -1,48 +1,67 @@
 type t = 
   {
     size : int;
-    contracted_shells : ContractedShell.t array;
+    contracted_shells : ContractedShell.t array ;
+    contracted_atomic_shells : ContractedAtomicShell.t array lazy_t;
   }
 
+module Ca = ContractedAtomicShell
 module Cs = ContractedShell
 module Gb = GeneralBasis
 module Ps = PrimitiveShell
 
 
 (** Returns an array of the basis set per atom *)
-let of_nuclei_and_general_basis n b =
-  let result = 
+let of_nuclei_and_general_basis nucl bas =
+  let index_ = ref 0 in
+  let contracted_shells = 
     Array.map (fun (e, center) ->
-        List.assoc e b
+        List.assoc e bas
         |> Array.map (fun (totAngMom, shell) ->
             let lc = 
               Array.map (fun Gb.{exponent ; coefficient} ->
                 coefficient, Ps.make totAngMom center exponent) shell
             in
-            Cs.make lc)
-      ) n
+            let result = Cs.make ~index:!index_ lc in
+            index_ := !index_ + Cs.size_of_shell result;
+            result
+          )
+    ) nucl
     |> Array.to_list
     |> Array.concat
   in
-  Array.iteri (fun i x ->
-      if (i > 0) then
-        result.(i) <- Cs.with_index x (
-            Cs.index result.(i-1)  +
-            Cs.size_of_shell result.(i-1) ) 
-    ) result ;
+  let contracted_atomic_shells = lazy(
+    let uniq_center_angmom =
+      Array.map (fun x -> Cs.center x, Cs.totAngMom x) contracted_shells
+      |> Array.to_list
+      |> List.sort_uniq compare
+    in
+    let csl = 
+      Array.to_list contracted_shells
+    in
+    List.map (fun (center, totAngMom) ->
+      let a = 
+        List.filter (fun x -> Cs.center x = center && Cs.totAngMom x = totAngMom) csl
+        |> Array.of_list
+      in
+      Ca.make ~index:(Cs.index a.(0)) a
+    ) uniq_center_angmom
+    |> List.sort (fun x y -> compare (Ca.index x) (Ca.index y))
+    |> Array.of_list
+  ) in
+  { contracted_shells ; contracted_atomic_shells ; size = !index_ }
 
-  let size =
-    let n = ref 0 in
-    for i=0 to (Array.length result) - 1 do
-      n := !n + Cs.size_of_shell result.(i)
-    done; !n
-  in
-  { contracted_shells = result ; size }
+
+let size x = x.size
+
+let contracted_atomic_shells x = Lazy.force x.contracted_atomic_shells
+
+let contracted_shells x = x.contracted_shells
 
 
 
 let to_string b =
-  let b = b.contracted_shells in
+  let b = contracted_atomic_shells b in
   let line ="
 -----------------------------------------------------------------------
 " in 
@@ -56,7 +75,7 @@ let to_string b =
 -----------------------------------------------------------------------
 "
   ^
-  ( Array.map (fun p -> Format.(fprintf str_formatter "%a" Cs.pp p;
+  ( Array.map (fun p -> Format.(fprintf str_formatter "%a" Ca.pp p;
     flush_str_formatter ())) b
     |> Array.to_list
     |> String.concat line
diff --git a/Basis/Basis.mli b/Basis/Basis.mli
index 070a758..5228cc8 100644
--- a/Basis/Basis.mli
+++ b/Basis/Basis.mli
@@ -1,14 +1,6 @@
 (** The atomic basis set is represented as an array of {!ContractedShell.t}. *)
 
-type t = private
-  {
-    size              : int    ; (** Number of contracted Gaussians *)
-    contracted_shells :
-      ContractedShell.t array ; (** Contracted shells *)
-  }
-
-val to_string : t -> string
-(** Pretty prints the basis set in a string. *)
+type t 
 
 
 val of_nuclei_and_general_basis : Nuclei.t -> GeneralBasis.t -> t
@@ -21,9 +13,24 @@ val of_nuclei_and_general_basis : Nuclei.t -> GeneralBasis.t -> t
 *)
 
 
-
 val of_nuclei_and_basis_filename : nuclei:Nuclei.t -> string -> t
 (** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t}
   from a file. 
   *)
 
+
+val size : t -> int
+(** Number of contracted basis functions. *)
+
+val contracted_atomic_shells : t -> ContractedAtomicShell.t array
+(** Returns the contracted basis functions per atom. *)
+
+val contracted_shells : t -> ContractedShell.t array
+(** Returns all the contracted basis functions. *)
+
+
+(** {2 Printers} *)
+
+val to_string : t -> string
+(** Pretty prints the basis set in a string. TODO *)
+
diff --git a/Basis/ContractedAtomicShell.ml b/Basis/ContractedAtomicShell.ml
new file mode 100644
index 0000000..570cb27
--- /dev/null
+++ b/Basis/ContractedAtomicShell.ml
@@ -0,0 +1,105 @@
+open Util
+open Constants
+
+type t = {
+  expo      : float array array;
+  coef      : float array array;
+  center    : Coordinate.t;
+  totAngMom : AngularMomentum.t;
+  norm_coef : float array array;
+  norm_coef_scale : float array;
+  index     : int;
+  contr     : ContractedShell.t array;
+}
+
+module Am = AngularMomentum
+module Co = Coordinate
+module Cs = ContractedShell
+
+
+let make ?(index=0) contr = 
+  assert (Array.length contr > 0);
+
+  let coef = Array.map Cs.coefficients contr
+  and expo = Array.map Cs.exponents contr
+  in
+
+  let center = Cs.center contr.(0) in
+  let rec unique_center = function
+  | 0 -> true
+  | i -> if Cs.center contr.(i) = center then unique_center (i-1) else false
+  in
+  if not (unique_center (Array.length contr - 1)) then
+    invalid_arg "ContractedAtomicShell.make Coordinate.t differ";
+    
+  let totAngMom = Cs.totAngMom contr.(0) in
+  let rec unique_angmom = function
+  | 0 -> true
+  | i -> if Cs.totAngMom contr.(i) = totAngMom then unique_angmom (i-1) else false
+  in
+  if not (unique_angmom (Array.length contr - 1)) then
+    invalid_arg "ContractedShell.make: AngularMomentum.t differ";
+    
+  let norm_coef =
+    Array.map Cs.normalizations contr 
+  in
+  let norm_coef_scale = Cs.norm_scales contr.(0)
+  in
+  { index ; expo ; coef ; center ; totAngMom ; norm_coef ;
+    norm_coef_scale ; contr  }
+
+
+let with_index a i =
+  { a with index = i }
+
+
+let exponents x = x.expo
+
+let coefficients x = x.coef
+
+let center x = x.center
+
+let totAngMom x = x.totAngMom
+
+let size x = Array.length x.contr
+
+let normalizations x = x.norm_coef
+
+let norm_scales x = x.norm_coef_scale
+
+let index x = x.index
+
+let size_of_shell x = Array.length x.norm_coef_scale
+
+let contracted_shells x = x.contr
+
+
+(** {2 Printers} *)
+
+open Format
+
+let pp_debug ppf x =
+  fprintf ppf "@[<2>{@ ";
+  fprintf ppf "@[<2>expo =@ %a ;@]@ " pp_float_2darray_size x.expo;
+  fprintf ppf "@[<2>coef =@ %a ;@]@ " pp_float_2darray_size x.coef;
+  fprintf ppf "@[<2>center =@ %a ;@]@ " Co.pp_angstrom x.center;
+  fprintf ppf "@[<2>totAngMom =@ %a ;@]@ " Am.pp_string x.totAngMom;
+  fprintf ppf "@[<2>norm_coef =@ %a ;@]@ " pp_float_2darray_size x.norm_coef;
+  fprintf ppf "@[<2>norm_coef_scale =@ %a ;@]@ " pp_float_array_size x.norm_coef_scale;
+  fprintf ppf "@[<2>index =@ %d ;@]@ " x.index;
+  fprintf ppf "}@,@]"
+
+let pp ppf s =
+  (match s.totAngMom with
+   | Am.S -> fprintf ppf "@[%3d@]    " (s.index+1)
+   | _    -> fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+(Array.length s.norm_coef_scale))
+  );
+  fprintf ppf "@[%a@ %a@] @[" Am.pp_string s.totAngMom Co.pp s.center;
+  Array.iter2 (fun e_arr c_arr -> fprintf ppf "[@[<v>";
+    Array.iter2 (fun e c -> fprintf ppf "@[%16.8e  %16.8e@]@;" e c)
+    e_arr c_arr;
+    fprintf ppf "@]]@;") s.expo s.coef;
+  fprintf ppf "@]"
+
+
+
diff --git a/Basis/ContractedAtomicShell.mli b/Basis/ContractedAtomicShell.mli
new file mode 100644
index 0000000..5d8c856
--- /dev/null
+++ b/Basis/ContractedAtomicShell.mli
@@ -0,0 +1,78 @@
+(** Set of contracted Gaussians differing only by the powers of x, y and z, with a
+    constant {!AngularMomentum.t}, all centered on the same center.
+
+{%
+\begin{align*}
+\chi_{n_x,n_y,n_z}(r) & = f(n_x,n_y,n_z) \sum_{j=1}^{n} \sum_{i=1}^{m} \mathcal{N}_{ij}\, d_{ij}\, g_{ij\,n_x,n_y,n_z}(r) \\
+        & = (x-X_A)^{n_x} (y-Y_A)^{n_y} (z-Z_A)^{n_z} f(n_x,n_y,n_z) \sum_{j=1}^{n} \sum_{i=1}^{m} \mathcal{N}_{ij}\, d_{ij}\, \exp \left( -\alpha_{ij} |r-R_A|^2 \right)
+\end{align*}
+%}
+
+where:
+
+- {% $g_{ij\,n_x,n_y,n_z}(r)$ %} is the i-th {!PrimitiveShell.t} of the j-th {!ContractedShell.t} 
+
+- {% $n_x + n_y + n_z = l$ %}, the total angular momentum
+
+- {% $\alpha_{ij}$ %} are the exponents (tabulated) of the j-th {!ContractedShell.t} 
+
+- {% $d_{ij}$ %} are the contraction coefficients of the j-th {!ContractedShell.t} 
+
+- {% $\mathcal{N}_{ij}$ %} is the normalization coefficient of the i-th primitive shell
+  ({!PrimitiveShell.norm_coef}) of the j-th {!ContractedShell.t} 
+
+- {% $f(n_x,n_y,n_z)$ %} is a scaling factor adjusting the normalization coefficient for the
+  particular powers of {% $x,y,z$ %} ({!PrimitiveShell.norm_coef_scale})
+
+*)
+
+type t 
+
+val make : ?index:int -> ContractedShell.t array -> t 
+(** Creates a contracted shell from a list of coefficients and primitives.  *)
+
+val with_index  : t -> int -> t
+(** Returns a copy of the contracted shell with a modified index. *)
+
+val index : t -> int
+(** Index in the basis set, represented as an array of contracted shells. *)
+
+val center : t -> Coordinate.t
+(** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %}. *)
+
+val totAngMom : t -> AngularMomentum.t
+(** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *)
+
+val size : t -> int
+(** Number of contracted functions, {% $n$ %} in the definition. *)
+
+val contracted_shells: t -> ContractedShell.t array
+(** Array of contracted gaussians *)
+
+val exponents : t -> float array array
+(** Array of exponents {% $\alpha_{ij}$ %}. The first index is the index of
+    the contracted function, and the second index is the index of the primitive.
+*)
+
+val coefficients : t -> float array array
+(** Array of contraction coefficients {% $d_{ij}$ %}. The first index is the index of
+    the contracted function, and the second index is the index of the primitive.
+*)
+
+val normalizations : t -> float array array
+(** Normalization coefficients {% $\mathcal{N}_{ij}$ %}. The first index is the index of
+    the contracted function, and the second index is the index of the primitive.
+*)
+
+val norm_scales : t -> float array
+(** Scaling factors {% $f(n_x,n_y,n_z)$ %}, given in the same order as
+    [AngularMomentum.zkey_array totAngMom]. *)
+
+val size_of_shell : t -> int
+(** Number of contracted functions in the shell: length of {!norm_coef_scale}. *)
+
+
+(** {2 Printers} *)
+
+val pp : Format.formatter -> t -> unit
+
diff --git a/Basis/ContractedShell.ml b/Basis/ContractedShell.ml
index 61b2525..914515a 100644
--- a/Basis/ContractedShell.ml
+++ b/Basis/ContractedShell.ml
@@ -1,16 +1,15 @@
 open Util
 open Constants
-open Coordinate
 
 type t = {
-  expo      : float array;       (** Array of exponents {% $\alpha_i$ %} *)
-  coef      : float array;       (** Array of contraction coefficients {% $d_i$ %} *)
-  center    : Coordinate.t;      (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %} *)
-  totAngMom : AngularMomentum.t; (** Total angular momentum : {% $l = n_x + n_y + n_z$ %} *)
-  norm_coef : float array;       (** Normalization coefficients of primitive functions {% $1/\mathcal{N}_i$ %} *)
-  norm_coef_scale : float array; (** Scaling factors {% $f_i$ %}, given in the same order as [AngularMomentum.zkey_array totAngMom]. *)
-  index    : int;                (** Index in the basis set, represented as an array of contracted shells. *)
-  prim : PrimitiveShell.t array;
+  expo      : float array;
+  coef      : float array;
+  center    : Coordinate.t;
+  totAngMom : AngularMomentum.t;
+  norm_coef : float array;
+  norm_coef_scale : float array;
+  index     : int;
+  prim      : PrimitiveShell.t array;
 }
 
 module Am = AngularMomentum
@@ -56,25 +55,25 @@ let with_index a i =
   { a with index = i }
 
 
-let expo x = x.expo
+let exponents x = x.expo
 
-let coef x = x.coef
+let coefficients x = x.coef
 
 let center x = x.center
 
 let totAngMom x = x.totAngMom
 
-let size x = Array.length x.coef
+let size x = Array.length x.prim
 
-let norm_coef x = x.norm_coef
+let normalizations x = x.norm_coef
 
-let norm_coef_scale x = x.norm_coef_scale
+let norm_scales x = x.norm_coef_scale
 
 let index x = x.index
 
 let size_of_shell x = Array.length x.norm_coef_scale
 
-let prim x = x.prim
+let primitives x = x.prim
 
 
 (** {2 Printers} *)
diff --git a/Basis/ContractedShell.mli b/Basis/ContractedShell.mli
index c09ff1e..af0c1d1 100644
--- a/Basis/ContractedShell.mli
+++ b/Basis/ContractedShell.mli
@@ -34,15 +34,8 @@ val make : ?index:int -> (float * PrimitiveShell.t) array -> t
 val with_index  : t -> int -> t
 (** Returns a copy of the contracted shell with a modified index. *)
 
-
-val expo : t -> float array
-(** Array of exponents {% $\alpha_i$ %}. *)
-
-val coef : t -> float array
-(** Array of contraction coefficients {% $d_i$ %}. *)
-
-val prim : t -> PrimitiveShell.t array
-(** Array of primitive gaussians *)
+val index : t -> int
+(** Index in the basis set, represented as an array of contracted shells. *)
 
 val center : t -> Coordinate.t
 (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %}. *)
@@ -51,20 +44,26 @@ val totAngMom : t -> AngularMomentum.t
 (** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *)
 
 val size : t -> int
-(** Number of contracted functions, {% $m$ %} in the definition. *)
+(** Number of primitive functions, {% $m$ %} in the definition. *)
 
-val norm_coef : t -> float array
+val primitives : t -> PrimitiveShell.t array
+(** Array of primitive gaussians *)
+
+val exponents : t -> float array
+(** Array of exponents {% $\alpha_i$ %}. *)
+
+val coefficients : t -> float array
+(** Array of contraction coefficients {% $d_i$ %}. *)
+
+val normalizations : t -> float array
 (** Normalization coefficients {% $\mathcal{N}_i$ %} of the primitive shells. *)
 
-val norm_coef_scale : t -> float array
+val norm_scales : t -> float array
 (** Scaling factors {% $f(n_x,n_y,n_z)$ %}, given in the same order as
     [AngularMomentum.zkey_array totAngMom]. *)
 
-val index : t -> int
-(** Index in the basis set, represented as an array of contracted shells. *)
-
 val size_of_shell : t -> int
-(** Number of contracted functions in the shell. *)
+(** Number of contracted functions in the shell: length of {!norm_coef_scale}. *)
 
 
 (** {2 Printers} *)
diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml
index 1005821..ab3d1ce 100644
--- a/Basis/ContractedShellPair.ml
+++ b/Basis/ContractedShellPair.ml
@@ -34,18 +34,18 @@ Format.printf "@[<2>shell_a :@ %a@]@;" Cs.pp s_a;
 Format.printf "@[<2>shell_b :@ %a@]@;" Cs.pp s_b;
 *)
 
-  let make = Psp.create_make_of (Cs.prim s_a).(0) (Cs.prim s_b).(0) in
+  let make = Psp.create_make_of (Cs.primitives s_a).(0) (Cs.primitives s_b).(0) in
 
   let shell_pairs = 
     Array.mapi (fun i p_a ->
-        let c_a = (Cs.coef s_a).(i) in
+        let c_a = (Cs.coefficients s_a).(i) in
         let make = make p_a in
         Array.mapi (fun j p_b ->
-            let c_b = (Cs.coef s_b).(j) in
+            let c_b = (Cs.coefficients s_b).(j) in
             let coef = c_a *. c_b in
             assert (coef <> 0.);
             let cutoff = cutoff /. abs_float coef in
-            coef, make p_b cutoff) (Cs.prim s_b)) (Cs.prim s_a)
+            coef, make p_b cutoff) (Cs.primitives s_b)) (Cs.primitives s_a)
     |> Array.to_list
     |> Array.concat
     |> Array.to_list
@@ -112,7 +112,7 @@ let cmp a b =
 (** The array of all shell pairs with their correspondance in the list
     of contracted shells.
 *)
-let of_basis basis =
+let of_contracted_shell_array basis =
   Array.mapi (fun i shell_a ->
     Array.mapi (fun j shell_b -> 
       create shell_a shell_b)
diff --git a/Basis/ContractedShellPair.mli b/Basis/ContractedShellPair.mli
index a457789..6edc98b 100644
--- a/Basis/ContractedShellPair.mli
+++ b/Basis/ContractedShellPair.mli
@@ -26,7 +26,7 @@ val create : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option
     [None].
 *)
 
-val of_basis : ContractedShell.t array -> t option array array
+val of_contracted_shell_array : ContractedShell.t array -> t option array array
 (** Creates all possible contracted shell pairs from the basis set.
     If the shell pair is not significant, sets the value to [None].
 *)
diff --git a/Basis/ERI.ml b/Basis/ERI.ml
index 4a7fca9..ad05f22 100644
--- a/Basis/ERI.ml
+++ b/Basis/ERI.ml
@@ -84,14 +84,14 @@ let of_basis basis =
     | _ -> assert false
   in
 
-  let n = basis.Bs.size
-  and shell = basis.Bs.contracted_shells 
+  let n = Bs.size basis
+  and shell = Bs.contracted_shells basis
   in
 
 
   (* Pre-compute all shell pairs *)
   let shell_pairs =
-    Csp.of_basis shell
+    Csp.of_contracted_shell_array shell
   in
 
   (* Pre-compute diagonal integrals for Schwartz *)
diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml
index 1b6d6cc..4271aca 100644
--- a/Basis/KinInt.ml
+++ b/Basis/KinInt.ml
@@ -127,8 +127,8 @@ let of_basis basis =
     | _ -> assert false
   in
 
-  let n     = basis.Bs.size 
-  and shell = basis.Bs.contracted_shells
+  let n     = Bs.size basis
+  and shell = Bs.contracted_shells basis
   in
 
   let result = Mat.create n n in
diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml
index b845948..e6d4f71 100644
--- a/Basis/NucInt.ml
+++ b/Basis/NucInt.ml
@@ -49,8 +49,8 @@ let of_basis_nuclei basis nuclei =
     | _ -> assert false
   in
 
-  let n     = basis.Bs.size 
-  and shell = basis.Bs.contracted_shells 
+  let n     = Bs.size basis
+  and shell = Bs.contracted_shells basis
   in
 
   let eni_array = Mat.create n n in
diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml
index 08a1ecb..afc1753 100644
--- a/Basis/Overlap.ml
+++ b/Basis/Overlap.ml
@@ -103,8 +103,8 @@ let of_basis basis =
     | _ -> assert false
   in
 
-  let n     = basis.Bs.size 
-  and shell = basis.Bs.contracted_shells 
+  let n     = Bs.size basis
+  and shell = Bs.contracted_shells basis
   in
 
   let result = Mat.create n n in
diff --git a/INSTALL.md b/INSTALL.md
index 4e2d32d..014a5f8 100644
--- a/INSTALL.md
+++ b/INSTALL.md
@@ -14,6 +14,14 @@ LaCAML is the OCaml binding to the LAPACK library.
 opam install lacaml
 ```
 
+To use MKL with LaCaml:
+
+```bash
+export  LACAML_LIBS="-L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_rt -lpthread -lm -ldl"
+opam install lacaml
+```
+
+
 # odoc-ltxhtml
 
 This plugin allows to embed equations in the documentation generated by Ocamldoc.
diff --git a/README.rst b/README.rst
index 90cbef1..c4c6ea0 100644
--- a/README.rst
+++ b/README.rst
@@ -7,3 +7,4 @@ Requirements
 * gmp : GNU Multiple Precision arithmetic library
 * odoc-ltxhtml : https://github.com/akabe/odoc-ltxhtml
 
+
diff --git a/Utils/Util.ml b/Utils/Util.ml
index eaac21a..1ac05d9 100644
--- a/Utils/Util.ml
+++ b/Utils/Util.ml
@@ -178,12 +178,22 @@ let xt_o_x ~o ~x =
 (** {2 Printers} *)
 
 let pp_float_array_size ppf a = 
-  Format.fprintf ppf "@[<2>[@ %d: " (Array.length a);
+  Format.fprintf ppf "@[<2>@[ %d:@[<2>" (Array.length a);
   Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a;
-  Format.fprintf ppf "]@]"
+  Format.fprintf ppf "]@]@]"
 
 let pp_float_array ppf a = 
   Format.fprintf ppf "@[<2>[@ ";
   Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a;
   Format.fprintf ppf "]@]"
 
+let pp_float_2darray ppf a = 
+  Format.fprintf ppf "@[<2>[@ ";
+  Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array f) a;
+  Format.fprintf ppf "]@]"
+
+let pp_float_2darray_size ppf a = 
+  Format.fprintf ppf "@[<2>@[ %d:@[" (Array.length a);
+  Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array_size f) a;
+  Format.fprintf ppf "]@]@]"
+
diff --git a/Utils/Util.mli b/Utils/Util.mli
index 4a911a6..27c42a4 100644
--- a/Utils/Util.mli
+++ b/Utils/Util.mli
@@ -83,3 +83,20 @@ val pp_float_array : Format.formatter -> float array -> unit
 ]}
 *)
 
+val pp_float_2darray_size : Format.formatter -> float array array -> unit
+(** Example:
+{[
+[
+  2:[ 6:  1.000000   1.732051   1.732051   1.000000   1.732051   1.000000 ]
+    [ 4:  1.000000   2.000000   3.000000   4.000000 ] ]
+]}
+*)
+
+val pp_float_2darray : Format.formatter -> float array array -> unit
+(** Example:
+{[
+[ [   1.000000   1.732051   1.732051   1.000000   1.732051   1.000000 ]
+  [   1.000000   2.000000   3.000000   4.000000 ] ]
+]}
+*)
+