diff --git a/Project.toml b/Project.toml
index 604a316..bd73673 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
diff --git a/src/ReCo.jl b/src/ReCo.jl
index deed03f..3514546 100644
--- a/src/ReCo.jl
+++ b/src/ReCo.jl
@@ -4,8 +4,7 @@ Pkg.activate(".")
 using Random, Distributions
 using JLD2
 using ProgressMeter
-using GLMakie
-using LaTeXStrings
+using GLMakie, ColorSchemes, LaTeXStrings
 using StaticArrays
 
 import Dates: now, CompoundPeriod, canonicalize
diff --git a/src/animation.jl b/src/animation.jl
index 0ba1125..6b4afc8 100644
--- a/src/animation.jl
+++ b/src/animation.jl
@@ -1,6 +1,8 @@
 function animate(sol::Solution, args, name_part::String; framerate::Int64=10)
     println("Generating animation...")
 
+    set_theme!(theme_black())
+
     fig = Figure()
     ax = Axis(
         fig[1, 1];
@@ -10,57 +12,56 @@ function animate(sol::Solution, args, name_part::String; framerate::Int64=10)
         ylabel=L"y",
     )
 
+    Colorbar(fig[1, 2]; limits=(0, 2), colormap=ColorSchemes.rainbow, label=L"\frac{\varphi}{\pi}")
+
     animation_path = "exports/$name_part.mkv"
 
     record(fig, animation_path; framerate=framerate) do io
-        old_cx = zeros(args.N)
-        old_cy = zeros(args.N)
-
-        new_cx = zeros(args.N)
-        new_cy = zeros(args.N)
-
         segments_x = zeros(2 * args.N)
         segments_y = zeros(2 * args.N)
 
-        done_first_iter = false
+        circles = Vector{Circle}(undef, args.N)
+        interaction_circles = Vector{Circle}(undef, args.N)
+        skin_circles = Vector{Circle}(undef, args.N)
+
+        colors = Vector{RGBAf}(undef, args.N)
+        interaction_colors = Vector{RGBAf}(undef, args.N)
+        skin_colors = Vector{RGBAf}(undef, args.N)
 
         @showprogress 0.5 for frame in 1:args.n_frames
-            new_cx .= get_c_x.(sol.position[frame])
-            new_cy .= get_c_y.(sol.position[frame])
-
-            if done_first_iter
+            if frame > 1
                 empty!(ax)
 
-                for i in 1:(args.N)
-                    segments_x[2 * i - 1] = old_cx[i]
-                    segments_x[2 * i] = new_cx[i]
+                @simd for i in 1:(args.N)
+                    segments_x[2 * i - 1] = sol.center[i, frame - 1][1]
+                    segments_x[2 * i] = sol.center[i, frame][1]
 
-                    segments_y[2 * i - 1] = old_cy[i]
-                    segments_y[2 * i] = new_cy[i]
+                    segments_y[2 * i - 1] = sol.center[i, frame - 1][2]
+                    segments_y[2 * i] = sol.center[i, frame][2]
                 end
 
                 linesegments!(ax, segments_x, segments_y; color=1:(args.N))
             else
-                done_first_iter = true
                 println("Started recording!")
             end
 
-            old_cx .= new_cx
-            old_cy .= new_cy
+            @simd for i in 1:args.N
+                circles[i] = Circle(Point2(sol.center[i, frame]), args.particle_diameter / 2)
+                interaction_circles[i] = Circle(Point2(sol.center[i, frame]), args.interaction_r)
+                skin_circles[i] = Circle(Point2(sol.center[i, frame]), args.skin_r)
 
-            poly!(
-                ax,
-                [Circle(Point2(x, y), args.particle_diameter / 2) for (x, y) in zip(new_cx, new_cy)];
-                color=1:(args.N),
-            )
+                color = get(ColorSchemes.rainbow, rem2pi(sol.φ[i, frame] / (2 * π), RoundDown))
 
-            # DEBUG BEGIN
-            for (x, y) in zip(new_cx, new_cy)
-                arc!(ax, Point2(x, y), args.interaction_r, 0.0, 2 * π)
+                colors[i] = RGBAf(color)
+                interaction_colors[i] = RGBAf(color, 0.2)
+                skin_colors[i] = RGBAf(color, 0.05)
             end
-            # DEBUG END
 
-            ax.title = "t = $(round(sol.t[frame], digits=2))"
+            poly!(ax, circles; color=colors)
+            poly!(ax, interaction_circles; color=interaction_colors)
+            poly!(ax, skin_circles; color=skin_colors)
+
+            ax.title = "t = $(round(sol.t[frame], digits=3))"
 
             recordframe!(io)
         end
diff --git a/src/data.jl b/src/data.jl
index 15469a4..b052a5a 100644
--- a/src/data.jl
+++ b/src/data.jl
@@ -1,6 +1,6 @@
 struct Solution
     t::Vector{Float64}
-    position::Matrix{Vector{Float64}}
+    center::Matrix{Vector{Float64}}
     φ::Matrix{Float64}
 
     function Solution(N::Int64, n_frames::Int64)
@@ -8,6 +8,26 @@ struct Solution
     end
 end
 
+function save_frame!(sol, frame, t, particles)
+    frame += 1
+
+    sol.t[frame] = t
+
+    @simd for p in particles
+        p_id = p.id
+
+        pre_pos = sol.center[p_id, frame]
+        
+        @simd for i in 1:2
+            pre_pos[i] = p.c[i]
+        end
+
+        sol.φ[p_id, frame] = p.φ
+    end
+
+    return frame
+end
+
 function save_data_jld(sol, args, name_part)
     println("Saving data...")
 
diff --git a/src/run.jl b/src/run.jl
index 590e835..dc189c4 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -28,6 +28,9 @@ function run(;
 
     skin_r = 1.1 * (2 * v * n_steps_before_verlet_list_update * δt + 2 * interaction_r)
 
+    integration_steps = floor(Int64, T / δt) + 1
+    n_steps_before_save = round(Int64, save_at / δt)
+
     args = (
         v = v,
         c₁ = 4 * ϵ * 6 * σ^6 * δt * μ,
@@ -42,12 +45,13 @@ function run(;
         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],
-        n_frames = floor(Int64, T / save_at) + 1,
+        n_frames = floor(Int64, integration_steps / n_steps_before_save) + 1,
     )
 
-    sol, end_time = simulate(args, save_at, δt, T, n_steps_before_verlet_list_update)
+    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)"
 
@@ -56,7 +60,7 @@ function run(;
     end
 
     if framerate > 0
-        generate_animation(sol, args, name_part; framerate=framerate)
+        animate(sol, args, name_part; framerate=framerate)
     end
 
     return (sol = sol, args = args, name_part = name_part)
diff --git a/src/simulation.jl b/src/simulation.jl
index 4a8d85a..4b0681e 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -55,46 +55,29 @@ function euler!(args)
     return nothing
 end
 
-function simulate(args, save_at::Float64, δt::Float64, T::Float64, n_steps_before_verlet_list_update::Int64)
-    sol = Solution(args.N, args.n_frames)
-    save_timer = save_at
-    frame = 0
+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.")
 
-    update_verlet_list_at = n_steps_before_verlet_list_update * δt
-    update_verlet_list_timer = update_verlet_list_at
-
-    @showprogress 0.2 for t in 0:δt:T
-        if save_timer >= save_at
-            frame += 1
-
-            sol.t[frame] = t
-
-            @simd for p in args.particles
-                p_id = p.id
-
-                pre_pos = sol.position[p_id, frame]
-                
-                @simd for i in 1:2
-                    pre_pos[i] = p.c[i]
-                end
-
-                sol.φ[p_id, frame] = p.φ
-            end
-
-            save_timer = 0.0
-        else
-            save_timer += δt
+    @showprogress 0.2 for (integration_step, t) in enumerate(0:δt:T)
+        if integration_step % n_steps_before_save == 0
+            frame = save_frame!(sol, frame, t, args.particles)
         end
 
-        if update_verlet_list_timer >= update_verlet_list_at
+        if integration_step % n_steps_before_verlet_list_update == 0
             update_verlet_list!(args)
-            
-            update_verlet_list_timer = 0.0
-        else
-            update_verlet_list_timer += δt
         end
 
         euler!(args)