using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, DifferentialEquations
gr();

V = x -> x^2/2
Vp = x -> x

λ⁰ = 2.3  # initial location
prob = ODEProblem((λ,_,t) -> -Vp(λ), λ⁰, (0.0, 10.0))
λ = solve(prob; reltol=1E-6);

@gif for t=0.0:0.05:7.0
    plot(V, range(-5, stop=5, length=100); label="potential", title="t = $t")
    scatter!([λ(t)] ,[0.0]; label="charge")
end

λ⁰ = [1.2,2.3]  # initial location
prob = ODEProblem((λ,_,t) -> [1/(λ[1] - λ[2]), 1/(λ[2] - λ[1])], λ⁰, (0.0, 10.0))
λ = solve(prob; reltol=1E-6);

@gif for t=0.0:0.05:7.0
    scatter(λ(t) ,zeros(2); label="charges", xlims=(-5,5), title="t = $t")
end

λ₀ = [1.2,2.3]  # initial location
prob = ODEProblem((λ,_,t) -> [1/(λ[1] - λ[2]) - Vp(λ[1]), 1/(λ[2] - λ[1]) - Vp(λ[2])], λ₀, (0.0, 10.0))
λ = solve(prob; reltol=1E-6);

@gif for t=0.0:0.05:7.0
    plot(V, range(-5, stop=5, length=100); label="potential", title="t = $t")
    scatter!(λ(t) ,zeros(2); label="charges", xlims=(-5,5), title="t = $t")
end

N = 100
λ⁰ = randn(N)  # initial location
prob = ODEProblem((λ,_,t) -> [sum(1 ./ (λ[k] .- λ[[1:k-1;k+1:end]])) - Vp(λ[k]) for k=1:N], λ⁰, (0.0, 10.0))
λ = solve(prob; reltol=1E-6);

@gif for t=0.0:0.05:7.0
    plot(V, range(-5, stop=5, length=100); label="potential", title="t = $t")
    scatter!(λ(t) ,zeros(N); label="charges", xlims=(-5,5), title="t = $t")
end

@gif for t=0.0:0.05:7.0
    scatter(λ(t)/sqrt(N) ,zeros(N); label="charges", xlims=(-2,2), title="t = $t")
end

histogram(λ(10.0)/sqrt(N); nbins=20, normalize=true, label="histogram of charges")
plot!(x -> sqrt(2-x^2)/(π), range(eps()-sqrt(2.0); stop=sqrt(2)-eps(), length=100), label="semicircle")

x = Fun(-sqrt(2) .. sqrt(2))
w = sqrt(2-x^2)/π

norm(hilbert(w) + x/π)

a = -0.1; b= 0.3;
w = Fun(x -> sqrt(2-x^2)/π, a .. b)
sum(w)  # integral of w(x) between a and b

length(filter(λ -> a ≤ λ ≤ b, λ(10.0)/sqrt(N)))/N # integral of w_N(x) between a and b 

x = Fun(-sqrt(2) .. sqrt(2))
w = sqrt(2-x^2)/π
z = 0.1+0.2im
cauchy(w, z)

(sum(1 ./((λ(10.0)/sqrt(N)) .- z))/N)/(2π*im)