diff --git a/.gitignore b/.gitignore
index 895b4eba6fb6804a90cfdda1070a688807490245..ac8618609e96bdc04799b0b15f34275ddd4080e5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,4 @@
 **__pycache__
 JupyterNotebooks/
+**.o
+**.so
diff --git a/src/biscotto/gaussian.f90 b/src/biscotto/gaussian.f90
new file mode 100644
index 0000000000000000000000000000000000000000..6f8daed2cfa5f304e8193ea6c182dcf844275ad3
--- /dev/null
+++ b/src/biscotto/gaussian.f90
@@ -0,0 +1,189 @@
+! Fortran library containing the performance-sensitive part of the code.
+! The main function takes care of the big fat integration loop
+
+! Notes section
+! To speed up shell lookup I can allocate an array of size (n_k, max_nshell, 3) where max_nshell = max(shell_count(k)) and init to zero. All smaller shells will have some extra zero padding at the end.
+
+module shell_utils
+
+  implicit none
+
+contains
+
+  function in_shell(Delta_k, k, q) result(in_check)
+
+    implicit none
+    real*8, intent(in) :: Delta_k, k
+    real*8, dimension(3), intent(in) :: q
+    real*8 :: half_bin, norm
+    logical*1 :: in_check
+
+    half_bin = Delta_k / 2.d0
+    norm = sqrt(sum(q*q))
+    in_check = (norm < k + half_bin) .and. (norm > k - half_bin)
+
+  end function in_shell
+
+  function shell_count(k_f, Delta_k, k) result(count)
+
+    implicit none
+    real*8, intent(in) :: k_f, Delta_k, k
+    real*8 :: r, half_bin
+    integer*8 :: k_norm, loop_bound, x, y, z, count
+
+    count = 0.d0
+    k_norm = nint(k / k_f)
+    half_bin = Delta_k / k_f / 2.d0
+    loop_bound = int(k / k_f + half_bin) + 1
+    do x = -loop_bound, loop_bound
+      do y = -loop_bound, loop_bound
+        do z = -loop_bound, loop_bound
+          r = sqrt(real(x, 8)**2 + real(y, 8)**2 + real(z, 8)**2)
+          if ((r < (k_norm + half_bin)) .and. (r > (k_norm - half_bin))) then
+            count = count + 1
+          end if
+        end do
+      end do
+    end do
+
+  end function shell_count
+
+  function get_shell(k_f, Delta_k, k, shell_count) result(shell)
+
+    implicit none
+    real*8, intent(in) :: k_f, Delta_k, k
+    integer*8, intent(in) :: shell_count
+    real*8, dimension(shell_count, 3) :: shell
+    real*8 :: r, half_bin
+    integer*8 :: k_norm, loop_bound, x, y, z, shell_index=1
+
+    shell = 0.d0
+    k_norm = nint(k / k_f)
+    half_bin = Delta_k / k_f / 2.d0
+    loop_bound = int(k / k_f + half_bin) + 1
+    do x = -loop_bound, loop_bound
+      do y = -loop_bound, loop_bound
+        do z = -loop_bound, loop_bound
+          r = sqrt(real(x, 8)**2 + real(y, 8)**2 + real(z, 8)**2)
+          if ((r < (k_norm + half_bin)) .and. (r > (k_norm - half_bin))) then
+            shell(shell_index, :) = (/x, y, z/) * k_f
+            shell_index = shell_index + 1
+          end if
+        end do
+      end do
+    end do
+
+  end function get_shell
+
+end module shell_utils
+
+module gaussian
+
+  use shell_utils
+  implicit none
+
+contains
+
+  function integral_kernel(l, p_1, p_2, p_3, W_0, W_2, W_4) result(element)
+
+    implicit none
+    integer*8, dimension(5), intent(in) :: l
+    real*8, dimension(3), intent(in) :: p_1, p_2, p_3
+    real*8, dimension(3), intent(in) :: w_0, w_2, w_4
+    real*8 :: element
+
+    ! Here I'll write the expression of the kernel for each value of the l's
+    ! For now, let's just write the all-0 case.
+
+    if ((l(3) == 0) .and. (l(4) == 0) .and. (l(5) == 0)) then
+      element = W_0(1) * W_0(2) * W_0(3)
+    end if
+
+  end function integral_kernel
+
+  function gaussian_integral( &
+    l, &
+    k_f, &
+    Delta_k, &
+    triangles_1, &
+    triangles_2, &
+    W_0, &
+    W_2, &
+    W_4, &
+    k_fft, &
+    sampling &
+  ) result(mix)
+    
+    implicit none
+    integer*8, dimension(5), intent(in) :: l
+    real*8, intent(in) :: k_f, Delta_k
+    real*8, intent(in) :: triangles_1(:, :), triangles_2(:, :)
+    real*8, intent(in) :: W_0(:, :, :), W_2(:, :, :, :), W_4(:, :, :, :)
+    real*8, intent(in) :: k_fft(:)
+    integer*8, intent(in) :: sampling
+    real*8 :: max_k
+    real*8, dimension(3, 3) :: q, p
+    real*8, dimension(:, :, :), allocatable :: shells
+    integer*8 :: k_num, max_shell_size
+    integer*8, dimension(:), allocatable :: shell_sizes
+    integer*8, dimension(6, 3) :: permutations
+    integer*8 :: i, i_t_1, i_t_2, i_q_1, i_q_2, i_perm
+    real*8 :: mix
+
+    ! Sanity check on matrices
+    if ((size(triangles_1, 2) /= 3) .and. (size(triangles_2, 2) /=3)) then
+      stop "[f90] ERROR - triangle configurations must be of size 3"
+    end if
+    do i = 1, 3
+      if ((size(k_fft) /= size(w_0, i)) & 
+      .or. (size(k_fft) /= size(w_2, i+1)) & 
+      .or. (size(k_fft) /= size(w_4, i+1))) then
+        stop "[f90] ERROR - Sizes of FFT kernels not matching FFT frequencies"
+      end if
+    end do
+    ! First step: prepare shells
+    ! Get largest shell size to allocate shell array
+    max_k = maxval(triangles_1)
+    k_num = nint(max_k / Delta_k)
+    max_shell_size = shell_count(k_f, Delta_k, max_k)
+    allocate(shells(k_num, max_shell_size, 3))
+    allocate(shell_sizes(k_num))
+    shells = 0.d0
+    shell_sizes = 0
+    ! Store shells
+    do i = 1, k_num
+      shell_sizes(i) = shell_count(k_f, Delta_k, i*Delta_k)
+      shells(i, 1:shell_sizes(i), :) = get_shell(k_f, Delta_k, i*Delta_k, shell_sizes(i))
+    end do
+    ! Define permutations
+    permutations = transpose( &
+      reshape( &
+        (/1, 2, 3,  1, 3, 2,  2, 3, 1,  2, 1, 3,  3, 1, 2,  3, 2, 1/), &
+        shape(transpose(permutations)) &
+      ) &
+    )
+    ! Loop over the two triangles
+    do i_t_1 = 1, size(triangles_1, 1)
+      do i_t_2 = 1, size(triangles_2, 1)
+        ! Loop over q_1 and q_2 (using the sampling size)
+        do i_q_1 = 1, sampling
+          i_k_1 = nint(triangles_1(i_t_1, 1) / Delta_k)
+          q(1, :) = shells(i_k_1)
+          do i_q_2 = 1, sampling
+          end do
+        end do
+      end do
+    end do
+
+
+
+
+
+
+
+  end function gaussian_integral
+
+end module gaussian
+
+
+  
\ No newline at end of file