Float = Union{Float32,Float64}

"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 this dimension and coordinate?"
which_index(n::Int, dim::Int, coord::Int)::Int = (dim - 1) * (n - 1) + coord

"k value of the given degree of freedom at the corresponding index, with coord=0 always returning 0"
get_k(n::Int, N::Int, i::CartesianIndex, dim::Int, coord::Int)::Int =
    coord == 0 ? 0 : i[which_index(n, dim, coord)] - N ÷ 2 - 1

"k value of the DOF at the specified cubic image"
get_shifted_k(n::Int, N::Int, i::CartesianIndex, dim::Int, coord::Int, image::Vector{Int})::Int =
    get_k(n, N, i, dim, coord) + N * image[dim]

"Difference of k values between two particles"
get_Δk(n::Int, N::Int, i::CartesianIndex, dim::Int, coord1::Int, coord2::Int, image::Vector{Int})::Int =
    get_k(n, N, i, dim, coord1) - get_shifted_k(n, N, i, dim, coord2, image)

"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
    coords = n - 1
    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 * coords)...)
    for image in images
        Threads.@threads for i in CartesianIndices(Vs)
            for coord1 in 1:coords
                for coord2 in 0:coord1-1
                    Δk² = 0
                    for dim in 1:d
                        Δk² += get_Δk(n, N, i, dim, coord1, coord2, image)^2
                    end
                    Vs[i] += V_twobody(Δk² * L²_over_N²)
                end
            end
        end
    end
    return Vs
end