diff --git a/Project.toml b/Project.toml
index b8e36b9..9fd25c5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
diff --git a/calculations/ACCC.jl b/calculations/ACCC.jl
new file mode 100644
index 0000000..94597f0
--- /dev/null
+++ b/calculations/ACCC.jl
@@ -0,0 +1,53 @@
+using Roots, LinearAlgebra, Plots
+
+include("../EC.jl")
+include("../common.jl")
+include("../p_space.jl")
+
+μ = 0.5
+V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
+
+# determining c0 with EC
+temp_c = range(1.1, 0.9, 3)
+p, w = get_mesh([0, 8], [256])
+H0 = get_T_matrix(p, μ)
+V = get_V_matrix(V_system(1), p, w)
+EC = affine_EC(H0, V, w)
+train!(EC, temp_c; ref_eval=-0.2, CAEC=false)
+quick_extrapolate(c) = minimum(abs2, get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, 0))
+c0 = find_zero(quick_extrapolate, 0.85)
+
+# Calculation of training and extrapolating E
+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_k = new_sqrt.(2μ .* training_E)
+
+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]
+
+order::Int = ceil((length(training_c) - 1) / 2) # order of the Pade approximant
+
+# Solve coefficients as a linear system
+M_left_element(c, i) = complex(c - c0)^(i/2)
+M_left = M_left_element.(training_c, (0:order)')
+M_right = -training_k .* M_left[:, 2:end] # remove the first column
+M = hcat(M_left, M_right) # M = [M_left | M_right]
+sol = M \ training_k
+a = sol[1:order+1]
+b = [1; sol[order+2:end]]
+
+# Pade approximant
+polynomial(a, c) = sum(i -> a[i+1] * complex(c - c0)^(i/2), 0:order)
+pade_approx(c) = polynomial(a, c) / polynomial(b, c)
+
+# Extrapolate
+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μ)
+
+# Plotting
+scatter(real.(training_E), imag.(training_E), label="training")
+scatter!(real.(exact_E), imag.(exact_E), label="exact")
+scatter!(real.(extrapolated_E), imag.(extrapolated_E), label="extrapolated", m=:star5)
diff --git a/common.jl b/common.jl
index 1201230..c5e537d 100644
--- a/common.jl
+++ b/common.jl
@@ -2,6 +2,9 @@ using LinearAlgebra, DelimitedFiles, SparseArrays
 
 @enum coordinate_system jacobi src
 
+"Square root function with the branch cut along the postive real axis"
+new_sqrt(x::Number)::ComplexF64 = im * sqrt(complex(-x))
+
 "Sum over array while minimizing catastrophic cancellation as much as possible"
 function better_sum(arr::Array{T}) where T<:Real
     pos_arr = arr[arr .> 0]