Branch ambiguity fixed
This commit is contained in:
parent
bb7b9cb198
commit
b01b6e5d6c
|
|
@ -20,7 +20,7 @@ c0 = find_zero(quick_extrapolate, 0.85)
|
||||||
# Calculation of training and extrapolating E
|
# Calculation of training and extrapolating E
|
||||||
training_c = range(1.2, 0.9, 9) # original: range(1.35, 0.9, 5)
|
training_c = range(1.2, 0.9, 9) # original: range(1.35, 0.9, 5)
|
||||||
training_E = [quick_pole_E(V_system(c)) for c in training_c]
|
training_E = [quick_pole_E(V_system(c)) for c in training_c]
|
||||||
training_k = new_sqrt.(2μ .* training_E)
|
training_k = alt_sqrt.(2μ .* training_E)
|
||||||
|
|
||||||
extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
|
extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
|
||||||
exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
|
exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
|
||||||
|
|
@ -28,7 +28,7 @@ exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
|
||||||
order::Int = ceil((length(training_c) - 1) / 2) # order of the Pade approximant
|
order::Int = ceil((length(training_c) - 1) / 2) # order of the Pade approximant
|
||||||
|
|
||||||
# Solve coefficients as a linear system
|
# Solve coefficients as a linear system
|
||||||
M_left_element(c, i) = complex(c - c0)^(i/2)
|
M_left_element(c, i) = alt_sqrt(c - c0)^i
|
||||||
M_left = M_left_element.(training_c, (0:order)')
|
M_left = M_left_element.(training_c, (0:order)')
|
||||||
M_right = -training_k .* M_left[:, 2:end] # remove the first column
|
M_right = -training_k .* M_left[:, 2:end] # remove the first column
|
||||||
M = hcat(M_left, M_right) # M = [M_left | M_right]
|
M = hcat(M_left, M_right) # M = [M_left | M_right]
|
||||||
|
|
@ -37,14 +37,11 @@ a = sol[1:order+1]
|
||||||
b = [1; sol[order+2:end]]
|
b = [1; sol[order+2:end]]
|
||||||
|
|
||||||
# Pade approximant
|
# Pade approximant
|
||||||
polynomial(a, c) = sum(i -> a[i+1] * complex(c - c0)^(i/2), 0:order)
|
polynomial(a, c) = sum(i -> a[i+1] * alt_sqrt(c - c0)^i, 0:order)
|
||||||
pade_approx(c) = polynomial(a, c) / polynomial(b, c)
|
pade_approx(c) = polynomial(a, c) / polynomial(b, c)
|
||||||
|
|
||||||
# Extrapolate
|
# Extrapolate
|
||||||
extrapolated_k = pade_approx.([training_c; extrapolating_c])
|
extrapolated_k = pade_approx.([training_c; extrapolating_c])
|
||||||
if real.(extrapolated_k[end]) < 0 # flip if following anti-resonance
|
|
||||||
extrapolated_k = -conj.(extrapolated_k)
|
|
||||||
end
|
|
||||||
extrapolated_E = (extrapolated_k .^ 2) / (2μ)
|
extrapolated_E = (extrapolated_k .^ 2) / (2μ)
|
||||||
|
|
||||||
# Plotting
|
# Plotting
|
||||||
|
|
|
||||||
|
|
@ -2,8 +2,8 @@ using LinearAlgebra, DelimitedFiles, SparseArrays
|
||||||
|
|
||||||
@enum coordinate_system jacobi src
|
@enum coordinate_system jacobi src
|
||||||
|
|
||||||
"Square root function with the branch cut along the postive real axis"
|
"Square root function with the branch cut along the postive imaginary axis"
|
||||||
new_sqrt(x::Number)::ComplexF64 = im * sqrt(complex(-x))
|
alt_sqrt(x::Number)::ComplexF64 = sqrt(im * x) / sqrt(im)
|
||||||
|
|
||||||
"Sum over array while minimizing catastrophic cancellation as much as possible"
|
"Sum over array while minimizing catastrophic cancellation as much as possible"
|
||||||
function better_sum(arr::Array{T}) where T<:Real
|
function better_sum(arr::Array{T}) where T<:Real
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue