tensor_element Function

real(dp) pure function tensor_element(i, j, k, r)

    !> order of the derivative wrt \(x\)-component of \(\mathbf{r}_{a}\)
    integer, intent(in) :: i
    !> order of the derivative wrt \(y\)-component of \(\mathbf{r}_{a}\)
    integer, intent(in) :: j
    !> order of the derivative wrt \(z\)-component of \(\mathbf{r}_{a}\)
    integer, intent(in) :: k
    !> coordinate vector \(\mathbf{r} = \mathbf{r}_{a}\) - \(\mathbf{r}_{b}\)
    real(dp), dimension(:), intent(in) :: r

    integer :: l, m, n, o, p
    real(dp) :: cl, cm, cn
    real(dp) :: norm

    tensor_element = 0.0_dp

    norm = norm2(r)

    do l = 0, i
        cl = tensor_coefficient(l,i,1) * (r(1) / norm)**l
        o = l + i + 1
        do m = 0, j
            cm = cl * tensor_coefficient(m,j,o) * (r(2) / norm)**m
            p = o + j + m
            do n = 0, k
                cn = cm * tensor_coefficient(n,k,p) * (r(3) / norm)**n
                tensor_element = tensor_element + cn
            end do
        end do
    end do

    tensor_element = tensor_element / norm**(i + j + k + 1)

end function tensor_element