yapanuwan
/
DVR-jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Almost correct

This commit is contained in:
ysyapa 2023-08-30 21:23:26 +00:00
parent fd262cd8c9
commit fa7fa0c09e
1 changed files with 15 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
common.jl
View File

@ -9,23 +9,25 @@ struct system{T}
n::Int n::Int
N::Int N::Int
L::T L::T
μs::Vector{Int} μs::Vector{T}
invU::Matrix{Int} invU::Matrix{T}
function system{T}(d::Int, n::Int, N::Int, L::Real) where {T<:Float} function system{T}(d::Int, n::Int, N::Int, L::Real) where {T<:Float}
μs = [Int((coord + 1)^2 * reducedMass(coord, 1)) for coord in 1:(n - 1)] μs = [1/((coord + 1)^2 * reducedMass(coord, 1)) for coord in 1:(n - 1)]
# TODO: Optimize # TODO: Optimize
invU = Matrix{Int}(undef, n, n - 1) U = Matrix{T}(undef, n, n)
for i in CartesianIndices(invU) for i in CartesianIndices(U)
if i[1] - 1 == i[2] if i[1] + 1 == i[2]
invU[i] = -i[2] U[i] = -i[1]
elseif i[1] > i[2] elseif i[1] >= i[2]
invU[i] = 0 U[i] = 1
else else
invU[i] = 1 U[i] = 0
end end
end end
U[n, :] .= 1/n
invU = inv(U)[:, 1:(n - 1)]
return new{T}(d, n, N, convert(T, L), μs, invU) return new{T}(d, n, N, convert(T, L), μs, invU)
end end
@ -44,7 +46,7 @@ end
which_index(s::system, dim::Int, coord::Int)::Int = (dim - 1) * (s.n - 1) + coord 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" "Get the distance to the nearest image of the particle"
function nearest(s::system, Δk::Int)::Int function nearest(s::system, Δk)
# TODO: Optimize # TODO: Optimize
while true while true
if Δk >= s.N ÷ 2 if Δk >= s.N ÷ 2
@ -67,9 +69,9 @@ function calculate_Vs(s::system{T}, V_twobody::Function, ϕ::T, n_image::Int)::A
rs = s.invU * xs rs = s.invU * xs
for p1 in 1:s.n for p1 in 1:s.n
for p2 in 1:(p1 - 1) for p2 in 1:(p1 - 1)
Δk = Array{Int}(undef, s.d) Δk = Array{T}(undef, s.d)
for dim in 1:s.d for dim in 1:s.d
Δk_temp = Int(rs[p1, dim] - rs[p2, dim]) Δk_temp = rs[p1, dim] - rs[p2, dim]
Δk[dim] = nearest(s, Δk_temp) Δk[dim] = nearest(s, Δk_temp)
end end
for image in images for image in images