Modified Jacobi (incorrect results)

common.jl

@@ -1,7 +1,7 @@
 Float = Union{Float32,Float64}
 
 norm_square(x) = sum(x .* x)
-reducedMass(m1::Float, m2::Float) = 1 / (1/m1 + 1/m2)
+reducedMass(m1, m2) = 1 / (1/m1 + 1/m2)
 
 "A few-body system defined by its physical parameters"
 struct system{T}
@@ -9,20 +9,23 @@ struct system{T}
     n::Int
     N::Int
     L::T
-    μs::Vector{T}
-    invU::Matrix{T}
+    μs::Vector{Int}
+    invU::Matrix{Int}
 
     function system{T}(d::Int, n::Int, N::Int, L::Real) where {T<:Float}
-        μs = collect(1:n) ./ collect(2:(n + 1))
+        μs = [Int((coord + 1)^2 * reducedMass(coord, 1)) for coord in 1:(n - 1)]
 
-        U = zeros(T, n, n)
-        for coord in 1:n
-            U[coord, 1:coord] .= 1 / coord
-            if coord != n
-                U[coord, coord + 1] = -1
+        # TODO: Optimize
+        invU = Matrix{Int}(undef, n, n - 1)
+        for i in CartesianIndices(invU)
+            if i[1] - 1 == i[2]
+                invU[i] = -i[2]
+            elseif i[1] > i[2]
+                invU[i] = 0
+            else
+                invU[i] = 1
             end
         end
-        invU = inv(U)[:, 1:(n - 1)]
         
         return new{T}(d, n, N, convert(T, L), μs, invU)
     end
@@ -41,13 +44,16 @@ end
 which_index(s::system, dim::Int, coord::Int)::Int = (dim - 1) * (s.n - 1) + coord
 
 "Get the distance to the nearest image of the particle"
-function nearest(s, Δk)
-    if Δk > s.N ÷ 2
-        return Δk - s.N
-    elseif Δk < -s.N ÷ 2
-        return Δk + s.N
-    else
-        return Δk
+function nearest(s::system, Δk::Int)::Int
+    # TODO: Optimize
+    while true
+        if Δk > s.N ÷ 2
+            Δk -= s.N
+        elseif Δk < -s.N ÷ 2
+            Δk += s.N
+        else
+            return Δk
+        end            
     end
 end
 
@@ -61,9 +67,10 @@ function calculate_Vs(s::system{T}, V_twobody::Function, ϕ::T, n_image::Int)::A
         rs = s.invU * xs
         for p1 in 1:s.n
             for p2 in 1:(p1 - 1)
-                Δk = Array{T}(undef, s.d)
+                Δk = Array{Int}(undef, s.d)
                 for dim in 1:s.d
-                    Δk[dim] = nearest(s, rs[p1, dim] - rs[p2, dim])
+                    Δk_temp = Int(rs[p1, dim] - rs[p2, dim])
+                    Δk[dim] = nearest(s, Δk_temp)
                 end
                 for image in images
                     Δk² = norm_square(Δk .- (s.N .* image))