using ApproxFun, Plots, ComplexPhasePortrait, ApproxFun,  SingularIntegralEquations,
        SpecialFunctions, DualNumbers, SO

using SingularIntegralEquations.HypergeometricFunctions
gr();

r = 0.9  # radius we are solving in, can be up to but not including 1
S = Taylor(Circle(r))
α = 0.3
a = z -> α/z + 1/(1-z)
ã = Fun( z -> a(z) - α/z, S)
za = Fun( z -> z*a(z), S)
z = Fun(S)
D = Derivative() : S → S
L = z*D + (α+1-za)
v = L \ ã

# v is analytic in disk of radius r
phaseplot(v, (-1,1), (-1,1))

# u has a branch 
u = z -> z^α + z^(α+1)*v(z)
phaseplot(u, (-1,1), (-1,1))

up = dualpart(u(dual(0.1,1.0))) # uses DualNumbers.jl to calculate derivative

up - a(0.1)u(0.1)

S = Taylor()
z = Fun(S)
D = Derivative() : S

a₋₁ = -0.6
b₋₂ = 0.3

a = z -> a₋₁/z + 1/(2-z)
ã = Fun( z -> a(z) -  a₋₁/z, S)
za = Fun( z -> z*a(z), S)

b = z -> b₋₂/z^2 + 0.4/z + 1/(3-z)
zb̃ = Fun(z -> z*b(z) - b₋₂/z, S)
z²b = Fun( z -> z^2*b(z), S)

α = ((1-a₋₁) +  sqrt((a₋₁-1)^2-4b₋₂))/2
@show α
α*(α-1) + α*a₋₁ + b₋₂

L = z^2*D^2 + z*(2(α+1) + za)*D + (α*(α+1) + za*(α+1) + z²b)
v = L \ (-α*ã - zb̃)

phaseplot(v, (-3,3), (-3,3))

u = z -> z^α + z^(α+1)*v(z)

phaseplot(u, (-3,3), (-3,3))

# dual numbers don't work to second degree so let's check numerically: 
uf = Fun(u, 0.1 .. 0.2)
norm(uf'' + Fun(a, space(uf))*uf' + Fun(b, space(uf))*uf)

α̃ = ((1-a₋₁) -  sqrt((a₋₁-1)^2-4b₋₂))/2
@show α̃
α̃*(α̃-1) + α̃*a₋₁ + b₋₂

L = z^2*D^2 + z*(2(α̃+1) + za)*D + (α̃*(α̃+1) + za*(α̃+1) + z²b)
ṽ = L \ (-α̃*ã - zb̃)

phaseplot(ṽ, (-3,3), (-3,3))

ũ = z -> z^α + z^(α+1)*ṽ(z)

phaseplot(ũ, (-3,3), (-3,3))