function animate(directory::String, sol::Solution, variables)
    println("Generating animation...")

    set_theme!(theme_black())

    fig = Figure(; resolution=(1080, 1080))
    ax = Axis(
        fig[1, 1];
        limits=(-variables.l, variables.l, -variables.l, variables.l),
        aspect=AxisAspect(1),
        xlabel=L"x",
        ylabel=L"y",
    )

    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 = "$directory/animation.mkv"

    record(fig, animation_path; framerate=variables.framerate) do io
        circles = Observable(Vector{Circle}(undef, variables.N))
        colors = Observable(Vector{RGBAf}(undef, variables.N))

        if variables.debug
            segments_x = Observable(zeros(2 * variables.N))
            segments_y = Observable(zeros(2 * variables.N))

            interaction_circles = Observable(Vector{Circle}(undef, variables.N))
            skin_circles = Observable(Vector{Circle}(undef, variables.N))

            interaction_colors = Observable(Vector{RGBAf}(undef, variables.N))
            skin_colors = Observable(Vector{RGBAf}(undef, variables.N))
        end

        @showprogress 0.6 for frame in 1:(variables.n_frames)
            @simd for i in 1:(variables.N)
                circles[][i] = Circle(
                    Point2(sol.center[i, frame]), variables.particle_diameter / 2
                )

                color = get(color_scheme, rem2pi(sol.φ[i, frame] / (2 * π), RoundDown))
                colors[][i] = RGBAf(color)

                if variables.debug
                    interaction_circles[][i] = Circle(
                        Point2(sol.center[i, frame]), variables.interaction_r
                    )
                    skin_circles[][i] = Circle(
                        Point2(sol.center[i, frame]), variables.skin_r
                    )

                    interaction_colors[][i] = RGBAf(color, 0.12)
                    skin_colors[][i] = RGBAf(color, 0.06)
                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]
                        segments_x[][2 * i] = sol.center[i, frame][1]

                        segments_y[][2 * i - 1] = sol.center[i, frame - 1][2]
                        segments_y[][2 * i] = sol.center[i, frame][2]
                    end

                    if frame == 2
                        linesegments!(ax, segments_x, segments_y; color=colors)
                    end
                end
            else # First frame
                poly!(ax, circles; color=colors)

                if variables.debug
                    poly!(ax, interaction_circles; color=interaction_colors)
                    poly!(ax, skin_circles; color=skin_colors)
                end

                println("Recording started!")
            end

            notify(circles)
            notify(colors)

            if variables.debug
                notify(interaction_circles)
                notify(interaction_colors)

                notify(skin_circles)
                notify(skin_colors)

                notify(segments_x)
                notify(segments_y)
            end

            ax.title = "t = $(round(sol.t[frame], digits=3))"

            recordframe!(io)
        end
    end

    println("Animation done and saved to $animation_path.")

    return nothing
end