pair_correlation

.gitignore

@@ -1,2 +1,4 @@
 Manifest.toml
-exports
+exports
+.ipynb_checkpoints
+*.ipynb

Project.toml

@@ -1,5 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"

analysis/pair_correlation.jl

@@ -0,0 +1,26 @@
+if splitdir(pwd())[2] == "analysis"
+    cd("..")
+end
+
+if splitdir(pwd())[2] != "ReCo"
+    error("You have to be in the main directeory ReCo!")
+else
+    include("src/ReCo.jl")
+    includet("src/analysis/pair_correlation_function.jl")
+end
+
+##
+
+using CairoMakie
+CairoMakie.activate!()
+
+##
+
+dirs = readdir("exports/2021_11_19"; join=true)
+
+##
+
+for dir in dirs
+    data = JLD2.load("$dir/data.jld2")
+    display(pair_correlation(data["sol"], data["variables"]).fig)
+end

src/ReCo.jl

@@ -16,6 +16,8 @@ import Base: push!, run, iterate
 using Revise
 using BenchmarkTools
 
+set_theme!(theme_black())
+
 includet("PreVector.jl")
 includet("Particle.jl")
 includet("setup.jl")
@@ -23,4 +25,4 @@ includet("simulation.jl")
 includet("data.jl")
 includet("animation.jl")
 includet("postprocess.jl")
-includet("run.jl")
+includet("run.jl")

src/analysis/pair_correlation_function.jl

@@ -0,0 +1,72 @@
+function plot_g(radius, g, variables)
+    fig = Figure()
+    ax = Axis(
+        fig[1, 1];
+        xticks=0:ceil(Int64, maximum(radius)),
+        yticks=0:ceil(Int64, maximum(g)),
+        xlabel=L"r",
+        ylabel=L"g(r)",
+        title="v = $(variables.v)",
+    )
+
+    scatterlines!(ax, radius, g; color=:white, markercolor=:red)
+
+    return fig
+end
+
+function pair_correlation(sol, variables)
+    n_r = 100
+    n_last_frames = 200
+    min_radius = variables.particle_diameter / 2
+    max_radius = 10.0 * variables.particle_diameter
+    dr = (max_radius - min_radius) / n_r
+
+    radius = range(min_radius, max_radius; length=n_r)
+    N_g = zeros(variables.N, n_r)
+
+    Threads.@threads for r_ind in 1:n_r
+        r = radius[r_ind]
+
+        @simd for i in 1:(variables.N)
+            for j in 1:(variables.N)
+                if i != j
+                    for k in 1:n_last_frames
+                        frame = variables.n_frames - k + 1
+                        c1 = sol.center[i, frame]
+                        c2 = sol.center[j, frame]
+
+                        r⃗₁₂ = SVector(c1[1] - c2[1], c1[2] - c2[2])
+
+                        r⃗₁₂ = minimum_image(r⃗₁₂; l=variables.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
+    end
+
+    g = zeros(n_r)
+
+    @simd for r_ind in 1:n_r
+        r = radius[r_ind]
+        tmp_g = 0.0
+
+        @tturbo for i in 1:(variables.N)
+            tmp_g += N_g[i, r_ind]
+        end
+
+        tmp_g *=
+            (2 * variables.l)^2 /
+            ((variables.N * n_last_frames) * variables.N * 2 * π * r * dr)
+        g[r_ind] = tmp_g
+    end
+
+    fig = plot_g(radius, g, variables)
+
+    return (; fig, radius, g)
+end

src/animation.jl

@@ -1,8 +1,6 @@
-function animate(directory::String, sol::Solution, variables)
+function animate(directory::String, sol::Solution, variables; framerate=0)
     println("Generating animation...")
 
-    set_theme!(theme_black())
-
     fig = Figure(; resolution=(1080, 1080))
     ax = Axis(
         fig[1, 1];
@@ -18,7 +16,14 @@ function animate(directory::String, sol::Solution, variables)
 
     animation_path = "$directory/animation.mkv"
 
-    record(fig, animation_path; framerate=variables.framerate) do io
+    record_framerate = 0
+    if framerate > 0
+        record_framerate = framerate
+    else
+        record_framerate = variables.framerate
+    end
+
+    record(fig, animation_path; framerate=record_framerate) do io
         circles = Observable(Vector{Circle}(undef, variables.N))
         colors = Observable(Vector{RGBAf}(undef, variables.N))
 
@@ -33,13 +38,14 @@ function animate(directory::String, sol::Solution, variables)
             skin_colors = Observable(Vector{RGBAf}(undef, variables.N))
         end
 
+        π2 = 2 * π
+        particle_radius = variables.particle_diameter / 2
+
         @showprogress 0.6 for frame in 1:(variables.n_frames)
             @simd for i in 1:(variables.N)
-                circles[][i] = Circle(
+                circles[][i] = Circle(Point2(sol.center[i, frame]), particle_radius)
-                    Point2(sol.center[i, frame]), variables.particle_diameter / 2
-                )
 
-                color = get(color_scheme, rem2pi(sol.φ[i, frame] / (2 * π), RoundDown))
+                color = get(color_scheme, rem2pi(sol.φ[i, frame], RoundDown) / π2)
                 colors[][i] = RGBAf(color)
 
                 if variables.debug
@@ -50,19 +56,23 @@ function animate(directory::String, sol::Solution, variables)
                         Point2(sol.center[i, frame]), variables.skin_r
                     )
 
-                    interaction_colors[][i] = RGBAf(color, 0.12)
+                    interaction_colors[][i] = RGBAf(color, 0.08)
-                    skin_colors[][i] = RGBAf(color, 0.06)
+                    skin_colors[][i] = RGBAf(color, 0.04)
                 end
             end
 
             if frame > 1
                 if variables.debug
                     @simd for i in 1:(variables.N)
-                        segments_x[][2 * i - 1] = sol.center[i, frame - 1][1]
+                        first_ind = 2 * i - 1
-                        segments_x[][2 * i] = sol.center[i, frame][1]
+                        second_ind = 2 * i
+                        frame_min_1 = frame - 1
 
-                        segments_y[][2 * i - 1] = sol.center[i, frame - 1][2]
+                        segments_x[][first_ind] = sol.center[i, frame_min_1][1]
-                        segments_y[][2 * i] = sol.center[i, frame][2]
+                        segments_x[][second_ind] = sol.center[i, frame][1]
+
+                        segments_y[][first_ind] = sol.center[i, frame_min_1][2]
+                        segments_y[][second_ind] = sol.center[i, frame][2]
                     end
 
                     if frame == 2

src/pair_correlation_function.jl

@@ -1,67 +0,0 @@
-function pair_correlation(frame=0, dr=0.1)
-    sol, args, filename = 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
-    radius = 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 = radius[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 = radius[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(radius, g)
-end
-
-function plot_g(radius, g)
-    fig = Figure(; resolution=(1080, 1080))
-    ax = Axis(fig[1, 1]; xlabel=L"r", ylabel=L"g(r)")
-
-    scatterlines!(ax, radius, g)
-
-    return (; fig, radius, g)
-end