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

a = 2.0
γ = Circle(a*exp(im*π/4), 0.1)
f = Fun(z -> z^3*sin(z)/(z^4+a^4), γ)
sum(f)/(2π*im), sin(a*exp(im*π/4))/4

γ = Circle(1, 0.1)
f = Fun(z -> (z+1)/(z^2-1)^2, γ)
sum(f)/(2π*im) # almost equals -1/4

a = 2.0

γ = Circle(a, 0.1)
f = Fun(z -> z^2*exp(z)/(z^3-a^3), γ)
sum(f)/(2π*im), exp(a)/3

θ = Fun(0 .. 2π)
sum(1/(5-4cos(θ))) -2π/3

θ = Fun(0 .. 2π)
sum(cos(2θ)/(5-4cos(θ)))/π

phaseplot(-3..3, -3..4, z-> 1/((z^2+1)*(z^2+4)))

x = Fun( Line())
sum(1/((x^2+1)*(x^2+4)))

π/6

phaseplot(-4..4, -4..4, z-> (z^2 - z + 2) / (z^4 + 10z^2 +9))

x = Fun(-1000 .. 1000)
sum(1/(x+im))

x = Fun(-100 .. 100)
sum(sin(2x)/(x^2+x+1))

-2π/sqrt(3) * sin(1)/exp(sqrt(3))

M = 200
x = Fun(-M .. M)
sum(cos(x)/(x^2+4)) # converges if we make M even bigger

π/(2*exp(2))

M = 2000
x = Fun(-M .. M)
sum(x*sin(x)/(x^2+1)) # Converges if we make M even bigger

π/ℯ

M = 10
ε = 1.0

plot(Segment(-M, -ε) ∪ Segment(ε, M);label="[-M,e] U [e,M]")
plot!(Arc(0.,ε, (π,0.)); label="C_e")
plot!(Arc(0., M, (0,π)); label = "C_M")

ε =0.001
M = 600.0
x = Fun(Segment(-M , -ε) ∪ Segment(ε, M))
a = 2.3; b = 3.8
sum((cos(a*x) - cos(b*x))/x^2) # Converges if we make M bigger

π*(b-a)

θ = Fun(0 .. 2π)
n = 4;
sum(cos(θ)^n)

π*factorial(1.0n)/(2^(n-1)*2*factorial(n/2))

k = 2.0
f = z -> exp(im*k*z)*sech(z)

-exp(-k*π)f(2.0)

f(2.0+im*π)

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

A = [1 2 -1; -2 2 1; 0 1 4]

λ = eigvals(A)
p = plot()
drawdisk!(1,3)
drawdisk!(2,3)
drawdisk!(4,1)
scatter!(complex.(λ); label="eigenvalues")
scatter!(complex.(diag(A)); label="diagonals")
p

λ = eigvals(A)
p = plot()
drawdisk!(1,2)
drawdisk!(2,3)
drawdisk!(4,2)
scatter!(complex.(λ); label="eigenvalues")
scatter!(complex.(diag(A)); label="diagonals")
p

λ = eigvals(A)
p = plot()
drawdisk!(2,3)
scatter!(complex.(λ); label="eigenvalues")
scatter!(complex.(diag(A)); label="diagonals")
p

n = 5
A = randn(n,n) 
A = A+ A' + 10I

λ, Q = eigen(A)

norm(A - Q*Diagonal(λ)*Q')

u₀ = randn(n)
v₀ = randn(n)
uv = solve(ODEProblem((uv,_,t) -> [uv[n+1:end]; A*uv[1:n]], [u₀; v₀], (0.,2.)); reltol=1E-10);

t = 2.0
uv(t)[1:n]

Q*Diagonal(cosh.(sqrt.(λ) .* t))*Q'*u₀ + Q*Diagonal(sinh.(sqrt.(λ) .* t) ./ sqrt.(λ))*Q'*v₀

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

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

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

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


n = 50
ellipse_f(z -> cosh(sqrt(z)*t), A, n, 10.0, 7.0, 2.0)*u₀ +
    ellipse_f(z -> sinh(sqrt(z)*t)/sqrt(z), A, n, 10.0, 7.0, 2.0)*v₀

phaseplot(-6..6, -3..3, z -> cosh(sqrt(z)))

phaseplot(-3..3, -3..3, z -> sinh(sqrt(z))/sqrt(z))

n = 5
A = randn(n,n) 
λ, V = eigen(A)

norm(A - V*Diagonal(λ)*inv(V))

u₀ = randn(n)
v₀ = randn(n)
uv = solve(ODEProblem((uv,_,t) -> [uv[n+1:end]; A*uv[1:n]], [u₀; v₀], (0.,2.)); reltol=1E-10);

t = 2.0
uv(t)[1:n]

V*Diagonal(cosh.(sqrt.(λ) .* t))*inv(V)*u₀ + 
   V*Diagonal(sinh.(sqrt.(λ) .* t) ./ sqrt.(λ))*inv(V)*v₀

n = 100
ellipse_f(z -> cosh(sqrt(z)*t), A, n, 0.0, 3.0, 3.0)*u₀ +
    ellipse_f(z -> sinh(sqrt(z)*t)/sqrt(z), A, n, 0.0, 3.0, 3.0)*v₀ 

f = z -> 2z/(4z-z^2-1)
f(exp(0.1im)) - 1/(2-cos(0.1))

phaseplot(-5..5, -5..5, f)

2+sqrt(3), 2- sqrt(3), 1/(2+sqrt(3))

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

g = θ -> 1/(2 - cos(θ))

Q = sum(Fun(g, 0 .. 2π))


err = Float64[ (
            (θ, w) = periodic_rule(n);
            sum(w.*g.(θ)) - Q
        ) for n=1:30];

N = length(err)

scatter(abs.(err), yscale=:log10, label="true error", title="Trapezium error bounds", legend=:bottomleft)
r = 1.1
scatter!(4/(4-r-inv(r)) .* r.^(-(1:N)) ./ (1 .- r.^(-(1:N))), label = "r = $r")
r = 3.0
scatter!(4/(4-r-inv(r)) .* r.^(-(1:N)) ./ (1 .- r.^(-(1:N))), label = "r = $r")
r = 2+sqrt(3)  - 0.01
scatter!(4/(4-r-inv(r)) .* r.^(-(1:N)) ./ (1 .- r.^(-(1:N))), label = "r = $r")