ReCo.jl/graphics/Luxor_graphics/periodic_boundary_conditions.jl

using Luxor
using Random: Random
using StaticArrays: SVector

using ReCo: ReCo

function gen_periodic_boundary_conditions_graphics()
    Random.seed!(23)

    drawing_box_length = 300

    graphics_export_dir = "exports/graphics"
    mkpath(graphics_export_dir)

    Drawing(
        drawing_box_length,
        drawing_box_length,
        "$graphics_export_dir/periodic_boundary_conditions.pdf",
    )
    origin()

    particle_radius = 10
    box_len = drawing_box_length / 3
    half_box_len = box_len / 2

    N_in_one_box = 5

    colors = ["red", "blue", "orange", "green", "cyan"]

    distance_limit = 0.99 * 2 * particle_radius

    particles = Vector{ReCo.Particle}(undef, 9 * N_in_one_box)
    translation_vecs = [
        SVector(0.0, 0.0),
        SVector(box_len, 0.0),
        SVector(box_len, box_len),
        SVector(0.0, box_len),
        SVector(-box_len, box_len),
        SVector(-box_len, 0.0),
        SVector(-box_len, -box_len),
        SVector(0.0, -box_len),
        SVector(box_len, -box_len),
    ]

    p1_c = SVector(0.0, 0.0)

    for p1_ind in 1:N_in_one_box
        while true
            x = (rand() - 0.5) * (box_len - 2 * particle_radius)
            y = (rand() - 0.5) * (box_len - 2 * particle_radius)

            p1_c = SVector(x, y)

            no_collision = true

            for p2_ind in 1:(p1_ind - 1)
                if ReCo.norm2d(p1_c - particles[p2_ind].c) < distance_limit
                    no_collision = false

                    break
                end
            end

            if no_collision
                break
            end
        end

        setcolor(colors[p1_ind])
        for box_ind in 1:9
            c = p1_c + translation_vecs[box_ind]
            particle = ReCo.Particle(0, c, 0.0)
            particles[p1_ind * box_ind] = particle
            circle(c[1], c[2], particle_radius, :fill)
        end
    end

    setcolor("black")
    setline(1.0)
    for line_factor in (-3, -1, 1, 3)
        line(
            Point(line_factor * drawing_box_length / 6, drawing_box_length / 2),
            Point(line_factor * drawing_box_length / 6, -drawing_box_length / 2),
            :stroke,
        )
        line(
            Point(drawing_box_length / 2, line_factor * drawing_box_length / 6),
            Point(-drawing_box_length / 2, line_factor * drawing_box_length / 6),
            :stroke,
        )
    end

    p1_ind = 3
    p2_ind = 5

    p1_c = particles[p1_ind].c
    p2_c = particles[p2_ind].c
    normal_vec_form_p1_to_p2 = p2_c - p1_c
    min_image_vec = ReCo.restrict_coordinates(normal_vec_form_p1_to_p2, half_box_len)
    min_image_vec_from_p1 = p1_c + min_image_vec

    setcolor("black")
    arrow(Point(p1_c[1], p1_c[2]), Point(p2_c[1], p2_c[2]); linewidth=2.0)

    setcolor("red")
    arrow(
        Point(p1_c[1], p1_c[2]),
        Point(min_image_vec_from_p1[1], min_image_vec_from_p1[2]);
        linewidth=2.0,
    )

    finish()

    return nothing
end