Random walk

analysis/mean_squared_displacement.jl

@@ -23,6 +23,54 @@ ReCo.restrict_coordinates!(::ReCo.Particle, ::Float64) = nothing
 ReCo.minimum_image_coordinate(value::Float64, ::Float64) = value
 ReCo.minimum_image(v::SVector{2,Float64}, ::Float64) = v
 
+const δt = 1e-4
+const Dₜ = ReCo.DEFAULT_Dₜ
+
+function fill_with_bundle_property!(
+    v::Vector, property::Symbol, sim_dir::String, first_bundle::Int64=1
+)
+    bundle_paths = ReCo.sorted_bundle_paths(sim_dir)
+
+    for i in first_bundle:length(bundle_paths)
+        bundle::ReCo.Bundle = JLD2.load_object(bundle_paths[i])
+
+        append!(v, getproperty(bundle, property))
+    end
+
+    return nothing
+end
+
+function max_possible_displacement(T::Float64, v₀::Float64, δt::Float64=δt, Dₜ::Float64=Dₜ)
+    return T * v₀ + T / δt * sqrt(2 * Dₜ * δt)
+end
+
+function msd_simulation(
+    v₀::Float64,
+    half_box_len::Float64,
+    T::Float64,
+    snapshot_at::Float64,
+    parent_dir::String,
+    comment::String="",
+)
+    dir = ReCo.init_sim(;
+        n_particles=1,
+        v₀=v₀,
+        parent_dir=parent_dir,
+        comment=comment,
+        half_box_len=half_box_len,
+    )
+
+    ReCo.run_sim(
+        dir;
+        duration=T,
+        seed=rand(1:typemax(Int64)),
+        snapshot_at=snapshot_at,
+        n_bundle_snapshots=1000,
+    )
+
+    return dir
+end
+
 function mean_squared_displacement(;
     n_simulations::Int64, v₀s::AbstractVector{Float64}, T::Float64
 )
@@ -30,9 +78,6 @@ function mean_squared_displacement(;
 
     n_v₀s = length(v₀s)
 
-    δt = 1e-4
-    Dₜ = ReCo.DEFAULT_Dₜ
-
     main_parent_dir = "mean_squared_displacement_$(Dates.now())"
 
     sim_dirs = Matrix{String}(undef, (n_simulations, n_v₀s))
@@ -40,42 +85,21 @@ function mean_squared_displacement(;
     progress = ProgressMeter.Progress(n_v₀s * n_simulations; dt=3, desc="MSD: ")
 
     for (v₀_ind, v₀) in enumerate(v₀s)
-        max_possible_displacement = T * v₀ + T / δt * sqrt(2 * Dₜ * δt)
+        half_box_len = max_possible_displacement(T, v₀)
 
         parent_dir = main_parent_dir * "/$v₀"
 
         Threads.@threads for sim_ind in 1:n_simulations
-            dir = ReCo.init_sim(;
-                n_particles=1,
-                v₀=v₀,
-                parent_dir=parent_dir,
-                comment="$sim_ind",
-                half_box_len=max_possible_displacement,
-            )
+            dir = msd_simulation(v₀, half_box_len, T, 0.5, parent_dir, "$sim_ind")
 
             sim_dirs[sim_ind, v₀_ind] = dir
 
-            ReCo.run_sim(
-                dir;
-                duration=T,
-                seed=rand(1:typemax(Int64)),
-                snapshot_at=0.5,
-                n_bundle_snapshots=1000,
-            )
-
             ProgressMeter.next!(progress; showvalues=[(:v₀, v₀)])
         end
     end
 
     ts = Float64[]
-
-    bundle_paths = ReCo.sorted_bundle_paths(sim_dirs[1, 1])
-
-    for i in 2:length(bundle_paths) # Skip the first bundle to avoid t = 0
-        bundle::ReCo.Bundle = JLD2.load_object(bundle_paths[i])
-
-        append!(ts, bundle.t)
-    end
+    fill_with_bundle_property!(ts, :t, sim_dirs[1, 1], 2) # Skip the first bundle to avoid t = 0
 
     mean_sq_displacements = zeros((length(ts), n_v₀s))
 
@@ -86,7 +110,7 @@ function mean_squared_displacement(;
             bundle_paths = ReCo.sorted_bundle_paths(sim_dir)
             snapshot_ind = 1
 
-            for i in 2:length(bundle_paths)
+            for i in 2:length(bundle_paths) # Skip the first bundle to avoid t = 0
                 bundle::ReCo.Bundle = JLD2.load_object(bundle_paths[i])
 
                 for c in bundle.c
@@ -113,21 +137,45 @@ function expected_mean_squared_displacement(t::Float64, v₀::Float64)
     return (4 * Dₜ + 2 * v₀^2 / Dᵣ) * t + 2 * v₀^2 * (exp(-Dᵣ * t) - 1) / (Dᵣ^2)
 end
 
+function init_cairomakie!()
+    CairoMakie.activate!()
+    set_theme!()
+
+    return nothing
+end
+
+function gen_figure()
+    text_width_in_pt = 405
+
+    return Figure(;
+        resolution=(text_width_in_pt, 0.55 * text_width_in_pt),
+        fontsize=10,
+        figure_padding=1,
+    )
+end
+
+function set_gaps!(fig::Figure)
+    colgap!(fig.layout, 5)
+    rowgap!(fig.layout, 5)
+
+    return nothing
+end
+
+function save_fig(filename::String, fig::Figure)
+    parent_dir = "exports/graphics"
+    mkpath(parent_dir)
+
+    return save("$parent_dir/$filename", fig; pt_per_unit=1)
+end
+
 function plot_mean_sq_displacement_with_expectation(
     ts::Vector{Float64},
     mean_sq_displacements::Matrix{Float64},
     v₀s::AbstractVector{Float64},
 )
-    CairoMakie.activate!()
-    set_theme!()
+    init_cairomakie!()
 
-    text_width_in_pt = 405
-
-    fig = Figure(;
-        resolution=(text_width_in_pt, 0.55 * text_width_in_pt),
-        fontsize=10,
-        figure_padding=1,
-    )
+    fig = gen_figure()
 
     ax = Axis(fig[1, 1]; xlabel=L"t", ylabel=L"\mathbf{MSD}", xscale=log10, yscale=log10)
 
@@ -149,18 +197,14 @@ function plot_mean_sq_displacement_with_expectation(
 
     Legend(fig[1, 2], v₀_scatter_plots, [L"v_0 = %$v₀" for v₀ in v₀s])
 
-    colgap!(fig.layout, 5)
-    rowgap!(fig.layout, 5)
+    set_gaps!(fig)
 
-    parent_dir = "exports/graphics/"
-    mkpath(parent_dir)
-
-    save("$parent_dir/mean_squared_displacement.pdf", fig; pt_per_unit=1)
+    save_fig("mean_squared_displacement.pdf", fig)
 
     return nothing
 end
 
-function run_analysis()
+function run_msd_analysis()
     v₀s = SVector(0.0, 20.0, 40.0, 60.0, 80.0)
 
     ts, mean_sq_displacements = mean_squared_displacement(;
@@ -172,4 +216,72 @@ function run_analysis()
     return nothing
 end
 
-# run_analysis()
+function plot_random_walk(T::Float64, v₀::Float64, seed::Int64)
+    Random.seed!(seed)
+
+    half_box_len = max_possible_displacement(T, v₀)
+    dir = msd_simulation(v₀, half_box_len, T, 8.0, "random_walk_$(Dates.now())")
+
+    ts = Float64[]
+    fill_with_bundle_property!(ts, :t, dir)
+
+    cs = SVector{2,Float64}[]
+    fill_with_bundle_property!(cs, :c, dir)
+
+    init_cairomakie!()
+
+    fig = gen_figure()
+
+    max_displacement = maximum(ReCo.norm2d.(cs))
+    expected_mean_displacement = sqrt(expected_mean_squared_displacement(ts[end], v₀))
+
+    limit = max(max_displacement, expected_mean_displacement)
+    limit *= 1.04
+
+    ax = Axis(
+        fig[1, 1];
+        xlabel=L"x",
+        ylabel=L"y",
+        aspect=AxisAspect(1),
+        limits=(-limit, limit, -limit, limit),
+    )
+
+    path_plot = lines!(ax, cs; linewidth=0.5)
+
+    marker_size = 4
+
+    start_point_plot = scatter!(ax, cs[1]; markersize=marker_size, color=:lime)
+
+    end_point_plot = scatter!(ax, cs[end]; markersize=marker_size, color=:red)
+
+    expected_mean_squared_displacement_plot = poly!(
+        ax,
+        Circle(Point2(0.0, 0.0), expected_mean_displacement);
+        strokecolor=:orange,
+        strokewidth=1,
+        color=RGBAf(0.0, 0.0, 0.0, 0.0),
+    )
+
+    Legend(
+        fig[1, 2],
+        [
+            start_point_plot,
+            end_point_plot,
+            path_plot,
+            expected_mean_squared_displacement_plot,
+        ],
+        ["Start point", "End point", "Path", "Expected mean displacement"],
+    )
+
+    set_gaps!(fig)
+
+    save_fig("random_walk.pdf", fig)
+
+    return nothing
+end
+
+function run_random_walk()
+    plot_random_walk(100_000.0, 0.0, 12345)
+
+    return nothing
+end