Electrostatic repulsion

using LinearAlgebra

using AbstractPlotting


clip11(x) = max(-1.0, min(1.0, x))

function repel(particles_node, N)
    particles = particles_node[]
    @inbounds for i in 1:N
        ftot = Vec3f0(0)
        p1 = particles[i]
        for j in 1:N
            if i != j
                p2 = particles[j]
                Δσ = acos(clip11(dot(p1, p2))) # great circle distance
                ftot += (p1 - p2)/max(1e-3, Δσ^2)
            end
        end
        particles[i] = normalize(p1 + 0.001 * ftot)
    end
    particles_node[] = particles
end

function addparticle!(particles, colors, nparticles)
    nparticles[] = nparticles[] + 1
    particles[][nparticles[]] = normalize(randn(Point3f0))
    colors[][nparticles[]] = to_color(:green)
    particles[] = particles[]
    colors[] = colors[]
end

s = Scene(show_axis = false)
mesh!(s, Sphere(Point3f0(0), 1f0), color = :gray)

max_particles = 5000
# Sadly, you currently can't resize 3D mesh particles, so we need to
# implement resize on our own...
particles = Node(fill(Point3f0(NaN), max_particles))
colors = Node(fill(RGBAf0(0, 0, 0, 0), max_particles))
meshscatter!(s, particles, color = colors, markersize = 0.05)
nparticles = Node(0)
for i=1:10
    addparticle!(particles, colors, nparticles)
end
update_cam!(s, FRect3D(Vec3f0(0), Vec3f0(1)))
s.center = false # don't reset the camera by display
N = 1000 # N gets replaced by 100 for testing
record(s, "output.mp4", 1:N) do iter
    isodd(iter) && addparticle!(particles, colors, nparticles)
    repel(particles, nparticles[])
end