open Util

type t = {
  expo     : float;
  expo_inv : float;
  center_ab: Coordinate.t;
  center   : Coordinate.t;
  norm_sq  : float;
  norm     : float;
  coef     : float;
  norm_fun : int array -> int array -> float;
  i        : int;
  j        : int;
}

exception Null_contribution

let create_array ?(cutoff=0.) p_a p_b =
  let log_cutoff =
    if (cutoff = 0.) then infinity
    else -. (log cutoff)
  in

  let center_ab = Coordinate.(
    Contracted_shell.center p_a |- Contracted_shell.center p_b )
  in
  let norm_sq =
    Coordinate.dot center_ab center_ab
  in
  Array.init (Contracted_shell.size p_a) (fun i ->
      let p_a_expo_center = Coordinate.(
        Contracted_shell.expo p_a i |. Contracted_shell.center p_a )
      in
      let f1 = 
        Contracted_shell.norm_coef p_a i
      in

      Array.init (Contracted_shell.size p_b) (fun j ->
          try
            let f2 =
              Contracted_shell.norm_coef p_b j
            in
            let norm_fun a b =
              f1 a *. f2 b
            in
            let norm =
              norm_fun
                [| Angular_momentum.to_int @@ Contracted_shell.totAngMom p_a ; 0 ; 0 |]
                [| Angular_momentum.to_int @@ Contracted_shell.totAngMom p_b ; 0 ; 0 |]
            in
            if (norm < cutoff) then
              raise Null_contribution;
            let p_b_expo_center = Coordinate.(
              Contracted_shell.expo p_b j |. Contracted_shell.center p_b )
            in
            let expo = Contracted_shell.(expo p_a i +. expo p_b j) in
            let expo_inv = 1. /. expo in
            let center = Coordinate.(
              expo_inv |. (p_a_expo_center |+ p_b_expo_center ) )
            in
            let argexpo =
              Contracted_shell.(expo p_a i *. expo p_b j) *. norm_sq *. expo_inv
            in
            if (argexpo > log_cutoff) then
              raise Null_contribution;
            let g =
              (pi *. expo_inv)**(1.5) *. exp(-. argexpo)
            in
            let norm_inv = 1./.norm in
            let norm_fun a b =
              norm_inv *. norm_fun a b
            in
            let coef =
              norm *. Contracted_shell.(coef p_a i *. coef p_b j) *. g
            in
            if (abs_float coef  < cutoff) then
              raise Null_contribution;
            Some { i ; j ; norm_fun ; norm ; coef ; expo ; expo_inv ; center ; center_ab ; norm_sq }
          with
          | Null_contribution -> None
        )
    )
  |> Array.to_list
  |> Array.concat
  |> Array.to_list
  |> List.filter (function Some _ -> true | None -> false)
  |> List.map (function Some x -> x | None -> assert false)
  |> Array.of_list