Electrostatic repulsion

using Makie
 using LinearAlgebra


 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