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

φ = z -> sqrt(z-1)sqrt(z+1)/(2im*(1+z^2)) + 
        sqrt(im-1)sqrt(im+1)/4*1/(z-im) - 
        sqrt(-im-1)sqrt(-im+1)/4*1/(z+im)

phaseplot(-3..3, -3..3, φ)

φ(200.0+200.0im)

φ(2.0+2.0im)

x = Fun()
cauchy(sqrt(1-x^2)/(1+x^2), 2.0+2.0im)

φ =  z -> (log(z-1)-log(z+1)) / ((2π*im)*(2+z)) - log(3)/(2π*im*(2+z))

phaseplot(-3..3, -3..3, φ)

φ(2.0+2.0im)

cauchy(1/(2+x), 2.0+2.0im)

f = x/sqrt(1-x^2)
π*hilbert(f, 0.1)

C = randn()
φ = z -> -z/(2*sqrt(z-1)*sqrt(z+1))+1/2 + C/(sqrt(z-1)*sqrt(z+1))

φ(0.1+0.0im)+φ(0.1-0.0im)

φ(1E8)

C = randn()
φ = z -> z/(sqrt(z-1)*sqrt(z+1)) + C/(sqrt(z-1)*sqrt(z+1)) 

φ(0.1+0.0im)+φ(0.1-0.0im)

φ(1E9)

z = 2.0+2.0im
cauchy(x^2, z)

(z^2*(log(z-1)-log(z+1))+2z)/(2π*im)

C = randn()
φ = z -> im/(sqrt(z-1)*sqrt(z+1)) * ((1-z^2)*(log(z-1)-log(z+1))-2z)/(2π*im) + C/(sqrt(z-1)sqrt(z+1))

φ(0.1+0.0im) + φ(0.1-0.0im) - sqrt(1-0.1^2)

φ(1E5)

φ = z -> 1/(2*(1+z^2)) + 
        sqrt(im-1)sqrt(im+1)*im/(4sqrt(z-1)sqrt(z+1))*1/(z-im) - 
        sqrt(-im-1)sqrt(-im+1)*im/(4sqrt(z-1)sqrt(z+1))*1/(z+im)

φ(1E5)*1E5^2

φ(0.1+0.0im) + φ(0.1-0.0im)

1/(1+0.1^2)

t = Fun()
z = 2.0+2.0im
cauchy(t, z), z*(log(z-1)-log(z+1))/(2π*im) + 1/(im*π)

x = 0.1
hilbert(t,x), x*(log(1-x)-log(1+x))/(π) + 2/π

t = Fun()
z = 2.0+2.0im
cauchy(sqrt(1-t^2)/(2+t), z), (sqrt(z-1)sqrt(z+1)-z)/(2im*(2+z)) + (sqrt(3)-2)/(2im*(z+2))

x = 0.1
hilbert(sqrt(1-t^2)/(2+t), x), -x/(2+x)+(sqrt(3)-2)/(2+x)

u = (t/(2+t)-(sqrt(3)-2)/(2+t)+(sqrt(3)-3)/3)/sqrt(1-t^2)
plot(u; ylims=(-5,5))

x = 0.1
hilbert(u,x) , 1/(2+x)

a = 1.0+2.0im
b = -1.0-2.0im

@show log(a*b)
@show log(a) + log(b)


scatter([real(1/b)], [imag(1/b)]; label="1/b")
scatter!([real(a)], [imag(a)]; label="a")
scatter!([0.0], [0.0]; label="0")
plot!(Segment(1, 1/b) ∪ Segment(1/b, a) ∪ Segment(a, 1); label="contour")

a = 1.0+2.0im
b = -2.0+2.0im

@show log(a*b)
@show log(a) + log(b) - 2π*im
scatter([real(1/b)], [imag(1/b)]; label="1/b")
scatter!([real(a)], [imag(a)]; label="a")
scatter!([0.0], [0.0]; label="0")
plot!(Segment(1, 1/b) ∪ Segment(1/b, a) ∪ Segment(a, 1); label="contour")

a = -2.0
b = -3.0

@show log(a*b)
@show log(a+0.0im) + log(b-0.0im)
@show log(a-0.0im) + log(b+0.0im)
@show log(a-0.0im) + log(b-0.0im) + 2π*im
@show log(a+0.0im) + log(b+0.0im) - 2π*im;

a = -2.0
b = 3.0

@show log(a*b +0.0im)
@show log(a+0.0im) + log(b);

@show log(a*b -0.0im)
@show log(a-0.0im) + log(b);

a = -2.0
b = 3.0 + im

@show log(a*b)
@show log(a-0.0im) + log(b);

b = 3.0 + im;
@show log(a*b)
@show log(a+0.0im) + log(b);

a = 1.0 + 1.0im
b = 1.0 + 1.0im
@show log(a*b +  eps()im) 
@show log((a+eps()im)*b)

a = 1.0 + 1.0im
b = -1.0 + 1.0im
@show log(a*b +  eps()im) 
@show log((a-eps()im)*b)

a = 1.0 + 1.0im
b = -1.0 + 1.0im
@show log(a*b -  eps()im) 
@show log(a)+log(b)-2π*im;

log1 = z -> begin
    if imag(z) > 0 
        log(z)
    elseif imag(z) == 0 && real(z) < 0
        log(z + 0.0im)
    elseif imag(z) < 0
        log(z) + 2π*im
    else
            error("log1 not defined on real axis")
    end
end

phaseplot(-3..3, -3..3, log1)

log1(-2.0+0.0im)

log1(-2.0-0.0im)

log1(-2.0)

log1(2.0+eps()im)

log1(2.0-eps()im)

θ = 0.1
v = (x,y) -> x^2 + y^2 < 1 ? 0 : imag(exp(-im*θ) * (x+im*y) + exp(im*θ) * (x+im*y)^(-1))

xx = yy = -5:0.01:5

contour(xx, yy, v.(xx', yy); nlevels=100)
plot!(Circle(); color=:black, legend=false)

θ =2.3
α = -θ/(2π)
κ = z -> (z-1)^(-α)*(z+1)^(α-1)
κ(0.1+0.0im) - exp(im*θ)*κ(0.1-0.0im)

κ(100.0)

phaseplot(-3..3, -3..3, κ)

θ =2.3
α = -θ/(2π)
κ = z -> (z-1)^(-α)*(z+1)^(α-1)


x = Fun()
κ₊ = exp(im*θ/2)*(1-x)^(-α)*(x+1)^(α-1)

f = Fun(exp)

z = 2+im
φ = z -> κ(z)*cauchy(f/κ₊, z)

φ(0.1+0.0im)-exp(im*θ)*φ(0.1-0.0im) - f(0.1)

β =  2.3;
x = -2.0
(x+0.0im)^(im*β) - exp(-2π*β)*(x-0.0im)^(im*β)

xx = yy = -1:0.01:1
contourf(xx, yy, abs.((xx' .+ im.*yy).^(im*β)))

θ =2.3
α = -θ/(2π)
κ = z -> (z-1)^(-α)*(z+1)^(α-1)
κ(0.1+0.0im) - exp(im*θ)*κ(0.1-0.0im)

r = 2.4
θ = 2.1
c = r*exp(im*θ)

α = -θ/(2π) + im*log(r)/(2π)

κ = z -> (z-1)^(-α)*(z+1)^(α-1)


κ(0.1+0.0im)-c*κ(0.1-0.0im)

V = x -> x^4
Vp = x -> 4x^3

N = 100
λ⁰ = randn(N)  # initial location
prob = ODEProblem((λ,t,_) -> Float64[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
    scatter(λ(t)/N^(1/4) ,zeros(N); label="charges", xlims=(-2,2), title="t = $t")
end

histogram(λ(10.0)/N^(1/4); nbins=30, normalize=true, label="histogram of charges")

b = 5
x = Fun(-b .. b)
(2im)cauchy(x^3*sqrt(b^2-x^2),z)

u = 4/(π*sqrt(b^2-x^2))*(-x^4 + b^2*x^2/2 + b^4/2)
hilbert(u, 0.1)

-4*0.1^3/π

plot(u; label = "b = $b")

sum(u)

cauchy(u, z)

z =100.0im;
2im/π*z^3 + 2im/π*(-z^4 + b^2*z^2/2 +b^4/2)/(sqrt(z-b)sqrt(z+b))

b = (2/3)^(1/4)
x = Fun(-b .. b)
u = 4/(π*sqrt(b^2-x^2))*(-x^4 + b^2*x^2/2 + b^4/2)
sum(u)

histogram(λ(10.0)/N^(1/4); nbins=25, normalize=true, label="histogram of charges")
plot!(u; label="u",xlims=(-2,2))