pair_correlation_funciton

src/PreVector.jl

@@ -17,3 +17,11 @@ function reset!(pv::PreVector{T}) where {T}
 
     return nothing
 end
+
+function iterate(pv::PreVector{T}, state=1) where {T}
+    if state > pv.last_ind
+        return nothing
+    else
+        return (pv.v[state], state + 1)
+    end
+end

src/ReCo.jl

@@ -8,7 +8,7 @@ using GLMakie, ColorSchemes, LaTeXStrings
 using StaticArrays
 
 import Dates: now, CompoundPeriod, canonicalize
-import Base: push!, run
+import Base: push!, run, iterate
 
 # Development deps
 using Revise

src/animation.jl

@@ -12,12 +12,9 @@ function animate(sol::Solution, args, name_part::String; framerate::Int64=10)
         ylabel=L"y",
     )
 
-    Colorbar(
-        fig[1, 2];
-        limits=(0, 2),
-        colormap=ColorSchemes.cyclic_mrybm_35_75_c68_n256_s25,
-        label=L"\frac{\varphi}{\pi}",
-    )
+    color_scheme = ColorSchemes.cyclic_mrybm_35_75_c68_n256_s25
+
+    Colorbar(fig[1, 2]; limits=(0, 2), colormap=color_scheme, label=L"\frac{\varphi}{\pi}")
 
     animation_path = "exports/$name_part.mkv"
 
@@ -42,9 +39,7 @@ function animate(sol::Solution, args, name_part::String; framerate::Int64=10)
                     Point2(sol.center[i, frame]), args.particle_diameter / 2
                 )
 
-                color = get(
-                    ColorSchemes.rainbow, rem2pi(sol.φ[i, frame] / (2 * π), RoundDown)
-                )
+                color = get(color_scheme, rem2pi(sol.φ[i, frame] / (2 * π), RoundDown))
                 colors[][i] = RGBAf(color)
 
                 if args.debug

src/pair_correlation_function.jl

@@ -0,0 +1,67 @@
+function pair_correlation(frame=0, dr=0.1)
+    sol, args, name_part = run(;
+        N=1000,
+        T=1.0,
+        v=20.0,
+        δt=1e-5,
+        save_at=0.1,
+        framerate=10,
+        n_steps_before_verlet_list_update=1000,
+    )
+
+    n_r = 100
+    rs = range(0, args.l; length=n_r)
+    N_g = zeros((args.N, n_r))
+
+    if frame == 0
+        frame = args.n_frames
+    end
+
+    for r_ind in 1:n_r
+        r = rs[r_ind]
+
+        for i in 1:(args.N)
+            for j in 1:(args.N)
+                if i != j
+                    c1 = sol.center[i, frame]
+                    c2 = sol.center[j, frame]
+
+                    r⃗₁₂ = SVector{2}(c1) - SVector{2}(c2)
+
+                    r⃗₁₂ = minimum_image(r⃗₁₂; l=args.l)
+
+                    distance = sqrt(r⃗₁₂[1]^2 + r⃗₁₂[2]^2)
+
+                    if (distance >= r) && (distance <= r + dr)
+                        N_g[i, r_ind] += 1
+                    end
+                end
+            end
+        end
+    end
+
+    g = zeros(n_r)
+
+    for r_ind in 1:n_r
+        r = rs[r_ind]
+        tmp_g = 0.0
+
+        for i in 1:(args.N)
+            tmp_g += N_g[i, r_ind]
+        end
+
+        tmp_g = tmp_g * (2 * args.l)^2 / (args.N^2 * 2 * π * r * dr)
+        g[r_ind] = tmp_g
+    end
+
+    return plot_g(rs, g)
+end
+
+function plot_g(rs, g)
+    fig = Figure(; resolution=(1080, 1080))
+    ax = Axis(fig[1, 1]; xlabel=L"r", ylabel=L"g(r)")
+
+    scatterlines!(ax, rs, g)
+
+    return (fig, rs, g)
+end

src/run.jl

@@ -8,8 +8,7 @@ function run(;
     save_data::Bool=false,
     n_steps_before_verlet_list_update::Int64=100,
     debug::Bool=false,
-    )
-
+)
     Random.seed!(42)
 
     μ = 1.0
@@ -33,27 +32,29 @@ function run(;
     n_steps_before_save = round(Int64, save_at / δt)
 
     args = (
-        v = v,
-        c₁ = 4 * ϵ * 6 * σ^6 * δt * μ,
-        c₂ = 2 * σ^6,
-        c₃ = sqrt(2 * D₀ * δt),
-        c₄ = sqrt(2 * Dᵣ * δt),
-        vδt = v * δt,
-        μ = μ,
-        interaction_r = interaction_r,
-        interaction_r² = interaction_r^2,
-        N = N,
-        l = l,
-        particle_diameter = particle_diameter,
-        particles = generate_particles(grid_n, grid_box_width, l),
-        skin_r = skin_r,
-        skin_r² = skin_r^2,
-        verlet_list = [PreVector(Int64, N - 1) for i in 1:(N - 1)],
-        n_frames = floor(Int64, integration_steps / n_steps_before_save) + 1,
-        debug = debug,
+        v=v,
+        c₁=4 * ϵ * 6 * σ^6 * δt * μ,
+        c₂=2 * σ^6,
+        c₃=sqrt(2 * D₀ * δt),
+        c₄=sqrt(2 * Dᵣ * δt),
+        vδt=v * δt,
+        μ=μ,
+        interaction_r=interaction_r,
+        interaction_r²=interaction_r^2,
+        N=N,
+        l=l,
+        particle_diameter=particle_diameter,
+        particles=generate_particles(grid_n, grid_box_width, l),
+        skin_r=skin_r,
+        skin_r²=skin_r^2,
+        verlet_list=[PreVector(Int64, N - 1) for i in 1:(N - 1)],
+        n_frames=floor(Int64, integration_steps / n_steps_before_save) + 1,
+        debug=debug,
     )
 
-    sol, end_time = simulate(args, δt, T, n_steps_before_verlet_list_update, n_steps_before_save)
+    sol, end_time = simulate(
+        args, δt, T, n_steps_before_verlet_list_update, n_steps_before_save
+    )
 
     name_part = "$(end_time)_T=$(T)_N=$(N)_v=$(v)_dt=$(δt)"
 
@@ -65,5 +66,5 @@ function run(;
         animate(sol, args, name_part; framerate=framerate)
     end
 
-    return (sol = sol, args = args, name_part = name_part)
+    return (sol=sol, args=args, name_part=name_part)
 end

src/simulation.jl

@@ -4,9 +4,9 @@ function update_verlet_list!(args)
     @simd for pv in args.verlet_list
         reset!(pv)
     end
-    
+
     for i in 1:(args.N - 1)
-        for j in (i + 1):args.N
+        for j in (i + 1):(args.N)
             p1 = args.particles[i]
             p2 = args.particles[j]
 
@@ -17,6 +17,8 @@ function update_verlet_list!(args)
             end
         end
     end
+
+    return nothing
 end
 
 function euler!(args)
@@ -24,17 +26,19 @@ function euler!(args)
         p = args.particles[i]
         verlet_list = args.verlet_list[p.id]
 
-        for j in 1:verlet_list.last_ind
+        for j in 1:(verlet_list.last_ind)
             p2 = args.particles[verlet_list.v[j]]
 
-            overlapping, r⃗₁₂, distance² = are_overlapping(p, p2, args.interaction_r², args.l)
+            overlapping, r⃗₁₂, distance² = are_overlapping(
+                p, p2, args.interaction_r², args.l
+            )
 
             if overlapping
                 c = args.c₁ / (distance²^4) * (args.c₂ / (distance²^3) - 1)
-                
+
                 @simd for k in 1:2
                     dck = c * r⃗₁₂[k]
-                    
+
                     p.tmp_c[k] -= dck
                     p2.tmp_c[k] += dck
                 end
@@ -59,19 +63,19 @@ function euler!(args)
     return nothing
 end
 
-function simulate(args,
+function simulate(
+    args,
     δt::Float64,
     T::Float64,
     n_steps_before_verlet_list_update::Int64,
     n_steps_before_save::Int64,
-    )
-
+)
     sol = Solution(args.N, args.n_frames)
 
     frame = 0
-    
+
     frame = save_frame!(sol, frame, 0.0, args.particles)
-    
+
     start_time = now()
     println("Started simulation at $start_time.")