open Util
open Constants

exception Null_contribution

type t =
{
  atomic_shell_a         : AtomicShell.t;
  atomic_shell_b         : AtomicShell.t;
  contracted_shell_pairs : ContractedShellPair.t list;
}


module Am  = AngularMomentum
module As  = AtomicShell
module Co  = Coordinate
module Cs  = ContractedShell
module Ps  = PrimitiveShell
module Psp = PrimitiveShellPair

(** Creates an atomic shell pair : an array of pairs of contracted shells.
*)
let make ?(cutoff=1.e-32) atomic_shell_a atomic_shell_b =

  let l_a = Array.to_list (As.contracted_shells atomic_shell_a)
  and l_b = Array.to_list (As.contracted_shells atomic_shell_b)
  in

  let contracted_shell_pairs =
    List.map (fun s_a -> 
      List.map (fun s_b ->
        Csp.make ~cutoff s_a s_b
      ) l_b
    ) l_a
  in


let atomic_shell_a  x = x.atomic_shell_a
let atomic_shell_b  x = x.atomic_shell_b
let contracted_shell_pairs = x.contracted_shell_pairs

let monocentric x = Csp.monocentric @@ List.hd x.contracted_shell_pairs

let center_ab x = Csp.center_ab @@ List.hd x.contracted_shell_pairs

let totAngMon x = Csp.totAngMon @@ List.hd x.contracted_shell_pairs

let norm_scales x = Csp.norm_scales @@ List.hd x.contracted_shell_pairs

let norm_sq x = Csp.norm_sq @@ List.hd x.contracted_shell_pairs

(** The array of all shell pairs with their correspondance in the list
    of contracted shells.
*)
let of_atomic_shell_array basis =
  Array.mapi (fun i shell_a ->
    Array.mapi (fun j shell_b -> 
      make ~cutoff shell_a shell_b)
    (Array.sub basis 0 (i+1)) 
  ) basis