using Plots, ComplexPhasePortrait, Interact, ApproxFun, Statistics, SpecialFunctions, LinearAlgebra
gr();

periodic_rule(n) = 2π/n*(0:(n-1)), 2π/n*ones(n)

function circle_rule(n, r) 
    θ = periodic_rule(n)[1]
    r*exp.(im*θ), 2π*im*r/n*exp.(im*θ)
end

function ellipse_rule(n, a, b) 
    θ = periodic_rule(n)[1]
    a*cos.(θ) + b*im*sin.(θ), 2π/n*(-a*sin.(θ) + im*b*cos.(θ))
end

f = θ ->  1/(2 + cos(θ))

errs = [((x, w) = periodic_rule(n); abs(sum(w.*f.(x)) - sum(Fun(f, 0 .. 2π)))) for n = 1:45];
plot(errs.+eps(); yscale=:log10, title="exponential convergence of n-point trapezium rule", legend=false, xlabel="n")

n=20
(x, w) = periodic_rule(n)
plot(Fun(f, 0 .. 2π); label = "integrand")
plot!(x, f.(x); label = "trapezium approximation")

ζ, w = circle_rule(20, 1.0)
scatter(ζ; title="quadrature points", legend=false, ratio=1.0)

ζ, w = circle_rule(20, 1.0)

f = z -> cos(z)

z = 0.1+0.2im

sum(f.(ζ)./(ζ .- z).*w)/(2π*im) - f(z)

f = z -> gamma(z)
fp = z -> gamma(z)polygamma(0,z) # exact derivative 
x = 1.2
fp_fd = [(h=2.0^(-n);  (f(x+h)-f(x))/h) for n = 1:50]
plot(abs.(fp_fd .- fp(x)); yscale=:log10, legend=false, title = "error of finite-difference with h=2^(-n)", xlabel="n")

trap_fp = [((ζ, w) = circle_rule(n, 0.5); 
        ζ .+= x; # circle around x
       sum(f.(ζ)./(ζ .- x).^2 .*w)/(2π*im)) for n=1:50]

plot(abs.(trap_fp .- fp.(x)); yscale=:log10, 
       title="Error of trapezium rule applied to Cauchy integral formula", xlabel="n", legend=false)

k=100
r = 1.0k
g = Fun( ζ -> exp(ζ)/(ζ - 0.1)^(k+1), Circle(0.1,r))
factorial(1.0k)/(2π*im) * sum(g) - exp(0.1)

θ = periodic_rule(100)[1]
plot(θ, real(f.(0.6exp.(im*θ) .+ x)./(0.5exp.(im*θ))))

scatter([1/5im,-1/5im]; label="singularities")
θ = range(0; stop=2π, length=2000)
a = 2; b= 0.1
plot!(a * cos.(θ) + im*b * sin.(θ); label="ellipse")

x = 0.1
f = z -> 1/(25z^2 + 1)

f_ellipse = [((z, w) = ellipse_rule(n, a, b); sum(f.(z)./(z.-x).*w)/(2π*im)) for n=1:1000]
plot(abs.(f_ellipse .- f(x)); yscale=:log10, title="convergence of n-point ellipse approximation", legend=false, xlabel="n")

f = z -> sqrt(z)

function sqrtⁿ(n,z,z₀) 
    ret = sqrt(z₀)
    c = 0.5/ret*(z-z₀)
    for k=1:n
        ret += c
        c *= -(2k-1)/(2*(k+1)*z₀)*(z-z₀)
    end
    ret
end


z₀ = 0.3
n = 40
phaseplot(-2..2, -2..2, z -> sqrtⁿ.(n,z,z₀))

(ζ, w) = ellipse_rule(20, 4.0, 1.0);
ζ  .= ζ .+ 4.5;
f_c = z ->  sum(f.(ζ)./(ζ.-z).*w)/(2π*im)
phaseplot(-2..10, -2 .. 2, f_c)

(ζ, w) = ellipse_rule(100, 4.0, 1.0);
ζ  .= ζ .+ 4.5;
f_c = z ->  sum(f.(ζ)./(ζ.-z).*w)/(2π*im)
f_c(5.3) - sqrt(5.3)

A = randn(5,5)
A = A + A'
λ = eigvals(A)

z,w = circle_rule(100,maximum(abs.(λ))+1)

plot(z)
scatter!(λ,zeros(5); label = "eigenvalues of A")

function circle_exp(A, n, z₀, r)
    z,w = circle_rule(n,r)
    z .+= z₀

    ret = zero(A)
    for j=1:n
        ret += w[j]*exp(z[j])*inv(z[j]*I - A)
    end

    ret/(2π*im)
end
   
circle_exp(A, 100, 0, 8.0) -exp(A) |>norm

function ellipse_exp(A, n, z₀, a, b)
    z,w = ellipse_rule(n,a,b)
    z .+= z₀

    ret = zero(A)
    for j=1:n
        ret += w[j]*exp(z[j])*inv(z[j]*I - A)
    end
    ret/(2π*im)
end


ellipse_exp(A, 50, 0, 8.0, 5.0) -exp(A) |>norm

function taylor_exp(A,n)
    ret = Matrix(I, size(A))
    for k=1:n
        ret += A^k/factorial(1.0k)
    end
    ret
end

B = A - 20I

taylor_exp(B, 200) -exp(B) |>norm

scatter(complex.(eigvals(B)))
plot!(ellipse_rule(50,8,5)[1] .- 20)

norm(ellipse_exp(B, 50, -20.0, 8.0, 5.0) - exp(B))

A = [1 2 3; 1 5 2; -4 1 6]

R = sum(abs.(A - Diagonal(diag(A))),dims=2)

drawcircle!(z0, R) = plot!(θ-> real(z0) + R[1]*cos(θ), θ-> imag(z0) + R[1]*sin(θ), 0, 2π, fill=(0,:red), α = 0.2, legend=false)

λ = eigvals(A)
p = plot()
for k = 1:size(A,1)
    drawcircle!(A[k,k], R[k])
end
scatter!(complex.(λ); label="eigenvalues")
scatter!(complex.(diag(A)); label="diagonals")
p

z₀ = (maximum(diag(A) .+ R) + minimum(diag(A) .- R)) /2 # average edges of circle
r = max(abs.(diag(A) .- R .- z₀)..., abs.(diag(A) .+ R .- z₀)...)

plot!(Circle(z₀, r))

g =Fun(θ -> 1/(2-cos(θ)), Laurent(-π .. π))
g₊ = g.coefficients[1:2:end]
scatter(abs.(g₊); yscale=:log10, label="|g_k|", legend=:bottomleft)
R = 1.1
scatter!(2/(4-R-inv(R))*R.^(-(0:length(g₊))), label = "R = $R")
R = 3.5
scatter!(2/(4-R-inv(R))*R.^(-(0:length(g₊))), label = "R = $R")
R = 2+sqrt(3)-0.1
scatter!(2/(4-R-inv(R))*R.^(-(0:length(g₊))), label = "R = $R")