is_naked_singularity: fails for certain metric configurations

[fjebaker]

function calc_exclusion(as, ϵs)
    regions = zeros(Float64, (length(as), length(ϵs)))
    Threads.@threads for i in eachindex(as)
        a = as[i]
        for (j, ϵ) in enumerate(ϵs)
            m = JohannsenPsaltisMetric(M = 1.0, a = a, ϵ3 = ϵ)
            regions[i, j] = if is_naked_singularity(m)
                NaN
            else
                Gradus.isco(m)
            end
        end
    end
    regions
end

as = range(0, 1.0, 120)
ϵs = range(-8, 10, 120)

img = @time calc_exclusion(as, ϵs)

Hits

    Stacktrace:
      [1] assert_bracket
        @ ~/.julia/packages/Roots/X7yil/src/Bracketing/bracketing.jl:52 [inlined]
      [2] init_state(::Roots.Bisection, F::Roots.Callable_Function{…}, x₀::Float64, x₁::Float64, fx₀::Float64, fx₁::Float64; m::Float64, fm::Float64)
        @ Roots ~/.julia/packages/Roots/X7yil/src/Bracketing/bisection.jl:50
      [3] init_state(::Roots.Bisection, F::Roots.Callable_Function{…}, x₀::Float64, x₁::Float64, fx₀::Float64, fx₁::Float64)
        @ Roots ~/.julia/packages/Roots/X7yil/src/Bracketing/bisection.jl:34
      [4] init_state(M::Roots.Bisection, F::Roots.Callable_Function{…}, x::Tuple{…})
        @ Roots ~/.julia/packages/Roots/X7yil/src/Bracketing/bracketing.jl:6
      [5] #init#42
        @ ~/.julia/packages/Roots/X7yil/src/find_zero.jl:299 [inlined]
      [6] init
        @ ~/.julia/packages/Roots/X7yil/src/find_zero.jl:289 [inlined]
      [7] solve(𝑭𝑿::Roots.ZeroProblem{…}, M::Roots.Bisection, p::Nothing; verbose::Bool, kwargs::@Kwargs{…})
        @ Roots ~/.julia/packages/Roots/X7yil/src/find_zero.jl:491
      [8] solve
        @ ~/.julia/packages/Roots/X7yil/src/find_zero.jl:484 [inlined]
      [9] #find_zero#39
        @ ~/.julia/packages/Roots/X7yil/src/find_zero.jl:220 [inlined]
     [10] find_zero (repeats 2 times)
        @ ~/.julia/packages/Roots/X7yil/src/find_zero.jl:210 [inlined]
     [11] find_zero
        @ ~/.julia/packages/Roots/X7yil/src/find_zero.jl:243 [inlined]
     [12] #isco#477
        @ ~/.julia/packages/Gradus/2k8ka/src/special-radii.jl:23 [inlined]
     [13] isco
        @ ~/.julia/packages/Gradus/2k8ka/src/special-radii.jl:14 [inlined]

[fjebaker]

It's really hard to find out why by catching the error because Roots.jl panics instead of throwing an error, but my few observations are that it's for negative deformation parameter mainly and high spins.

[fjebaker]

This turned out to be a read herring; the ISCO calculations are working fine, but it's the is_naked_singularity guard that is failing, returning false positives and false negatives. That algorithm needs to be revisited.