Float = Union{Float32,Float64}

norm_square(x::Array{Int})::Int = sum(x .* x)

"Eq (46): Partial derivative matrix element for 1 degree of freedom"
function ∂_1DOF(L::T, N::Int, k::Int, l::Int)::Complex{T} where {T<:Float}
    if k == l
        return -im * (π / L)
    else
        return (π / L) * (-1)^(k - l) * exp(-im * π * (k - l) / N) / sin(π * (k - l) / N)
    end
end

"Which index (dimension of the multidimensional array) corresponds to spatial dimension 'dim' and particle 'p'?"
which_index(n::Int, dim::Int, p::Int)::Int = (dim - 1) * (n - 1) + p

"Δk (distance in terms of lattice paramter) between two particles along the given dimension"
function get_Δk(n::Int, N::Int, i::CartesianIndex, dim::Int, p1::Int, p2::Int)::Int
    if p1 == p2
        return 0
    elseif p1 == n
        return -(i[which_index(n, dim, p2)] - N ÷ 2 - 1)
    elseif p2 == n
        return i[which_index(n, dim, p1)] - N ÷ 2 - 1
    else
        return i[which_index(n, dim, p1)] - i[which_index(n, dim, p2)]
    end
end

"Calculate diagonal elements of the V matrix"
function calculate_Vs(V_twobody::Function, d::Int, n::Int, N::Int, L::T, n_image::Int)::Array{T} where {T<:Float}
    L²_over_N² = (L / N)^2
    images = collect.(Iterators.product(fill(-n_image:n_image, d)...)) # TODO: Learn how to use tuples instead of vectors
    Vs = zeros(T, fill(N, d * (n - 1))...)
    Threads.@threads for i in CartesianIndices(Vs)
        for p1 in 1:n
            for p2 in (p1 + 1):n
                min_Δk = Array{Int}(undef, d)
                for dim in 1:d
                    Δk = get_Δk(n, N, i, dim, p1, p2)
                    if Δk > N ÷ 2
                        min_Δk[dim] = Δk - N
                    elseif Δk < -N ÷ 2
                        min_Δk[dim] = Δk + N
                    else
                        min_Δk[dim] = Δk
                    end
                end
                for image in images
                    Δk² = norm_square(min_Δk .- (N .* image))
                    Vs[i] += V_twobody(Δk² * L²_over_N²)
                end
            end
        end
    end
    return Vs
end