From 550727d150e67f82563a984e1ffc87b951283579 Mon Sep 17 00:00:00 2001
From: Fergus Baker <fergusbkr@gmail.com>
Date: Fri, 29 Apr 2022 10:00:25 +0100
Subject: [PATCH 1/3] fixed impact parameter mapping in second order

---
 src/GeodesicTracer/utility.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/GeodesicTracer/utility.jl b/src/GeodesicTracer/utility.jl
index 0f95b745..dbdd5dde 100644
--- a/src/GeodesicTracer/utility.jl
+++ b/src/GeodesicTracer/utility.jl
@@ -32,6 +32,6 @@ function map_impact_parameters(
 end
 
 function alpha_beta_to_vel(m::AbstractMetricParams{T}, u, α, β) where {T}
-    reg = u[2]^2
-    -1.0, -β / reg, -α / reg
+    mcomp = metric_components(m, @view(u[2:3]))
+    -1.0, -β / mcomp[3], -α / √(mcomp[3] * mcomp[4])
 end

From 1ac4f2090e97fa19a5318bed782ff5b70334da4f Mon Sep 17 00:00:00 2001
From: Fergus Baker <fergusbkr@gmail.com>
Date: Fri, 29 Apr 2022 10:04:16 +0100
Subject: [PATCH 2/3] added metric with refractive index

---
 src/metrics/kerr-refractive-ad.jl | 56 +++++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)
 create mode 100644 src/metrics/kerr-refractive-ad.jl

diff --git a/src/metrics/kerr-refractive-ad.jl b/src/metrics/kerr-refractive-ad.jl
new file mode 100644
index 00000000..dbb3a439
--- /dev/null
+++ b/src/metrics/kerr-refractive-ad.jl
@@ -0,0 +1,56 @@
+module __KerrRefractiveAD
+using ..StaticArrays
+
+@inline Σ(r, a, θ) = r^2 + a^2 * cos(θ)^2
+@inline Δ(r, R, a) = r^2 - R * r + a^2
+
+@fastmath function metric_components(M, a, n, corona_radius, rθ)
+    (r, θ) = rθ
+    R = 2M
+    Σ₀ = Σ(r, a, θ)
+    sinθ2 = sin(θ)^2
+
+    Δ₀ = Δ(r, R, a)
+
+    tt = -(1 - (R * r) / Σ₀)
+    rr = Σ₀ / Δ₀
+    θθ = Σ₀
+    ϕϕ = sinθ2 * (r^2 + a^2 + (sinθ2 * R * r * a^2) / Σ₀)
+
+    tϕ = (-R * r * a * sinθ2) / Σ₀
+
+    # need to interpolate the refractve index boundary
+    # else we don't have a gradient to compute, and we miss it
+    # hemorraging energy in the process
+    δr = 2.0 
+    if r ≤ corona_radius
+        # ...
+    elseif corona_radius ≤ r ≤ corona_radius + δr
+        t = (r - corona_radius) / δr
+        # use an arbitrarily steep smooth interpolation
+        # this one isn't perfect, but does a good job
+        k = atan(1e5t) * 2 / π
+        n = k + n * (1-k) 
+    else 
+        n = 1.0
+    end
+
+    @SVector [tt/n^2, rr, θθ, ϕϕ, tϕ/n]
+end
+
+end # module
+@with_kw struct KerrRefractiveAD{T} <: AbstractAutoDiffStaticAxisSymmetricParams{T}
+    @deftype T
+    M = 1.0
+    a = 0.0
+    n = 1.0
+    corona_radius = 20.0
+end
+
+# implementation
+GeodesicTracer.metric_components(m::KerrRefractiveAD{T}, rθ) where {T} =
+    __KerrRefractiveAD.metric_components(m.M, m.a, m.n, m.corona_radius, rθ)
+GradusBase.inner_radius(m::KerrRefractiveAD{T}) where {T} = m.M + √(m.M^2 - m.a^2)
+
+# special radii
+isco(m::KerrRefractiveAD{T}) where {T} = __BoyerLindquistFO.isco(m.M, m.a)

From 96053ebf064bbff8f4af9e6552c534e8d039504d Mon Sep 17 00:00:00 2001
From: Fergus Baker <fergusbkr@gmail.com>
Date: Fri, 29 Apr 2022 10:04:48 +0100
Subject: [PATCH 3/3] added exports

---
 src/Gradus.jl          | 2 +-
 src/metrics/metrics.jl | 1 +
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/Gradus.jl b/src/Gradus.jl
index 8edaefd9..4a906583 100644
--- a/src/Gradus.jl
+++ b/src/Gradus.jl
@@ -70,7 +70,7 @@ export AbstractCoronaModel,
 # pre-defined metrics
 include("metrics/metrics.jl")
 
-export BoyerLindquistAD, BoyerLindquistFO, JohannsenAD, MorrisThorneAD
+export BoyerLindquistAD, BoyerLindquistFO, JohannsenAD, MorrisThorneAD, KerrRefractiveAD
 
 # downstream modules and work
 include("special-radii.jl")
diff --git a/src/metrics/metrics.jl b/src/metrics/metrics.jl
index 5db91625..7dc5a16f 100644
--- a/src/metrics/metrics.jl
+++ b/src/metrics/metrics.jl
@@ -2,3 +2,4 @@ include("boyer-lindquist-ad.jl")
 include("boyer-lindquist-fo.jl")
 include("johannsen-ad.jl")
 include("morris-thorne-ad.jl")
+include("kerr-refractive-ad.jl")
\ No newline at end of file