Initial Refactor (#1)

* minimally working

* moved exports top level of each subdir

* export all at module level

* added common metrics

* most work to date integrated

* added missing deps

* added CITATION, missing doi and date released

* various bugfixes

* some smoke tests and formatting

* tolerance fix in tests

* updated readme

* typo corrections

* bugfixes for disc profile tracing

* added new version of disc profile code

* lamp post model + domain & generators tests

* typo fixes

* fixed test tolerances

* accretion formulae bugfixes

* redshift pointfunction test

* actually include tests in runner

* punctuation fix

CITATION.cff

@@ -0,0 +1,12 @@
+cff-version: 1.2.0
+message: "If you use this software, please cite it as below."
+authors:
+  - family-names: Baker
+    given-names: Fergus
+  - family-names: Young
+    given-names: Andrew
+title: "Gradus.jl"
+version: 0.1.0
+doi:
+date-released:
+url: "https://github.com/astro-group-bristol/Gradus.jl"

Project.toml

@@ -0,0 +1,19 @@
+name = "Gradus"
+uuid = "c5b7b928-ea15-451c-ad6f-a331a0f3e4b7"
+authors = ["fjebaker <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+GeometricalPredicates = "fd0ad045-b25c-564e-8f9c-8ef5c5f21267"
+GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"

README.md

@@ -1,2 +1,37 @@
 # Gradus.jl
-Spacetime generic relativistic ray tracing toolkit, leveraging AD and CAS, capable of calculating reverberations lags around accreting black holes, and other observational signatures.
+
+A spacetime generic, relativistic ray tracing toolkit, leveraging AD and CAS, capable of calculating reverberation lags around accreting black holes and observational signatures.
+
+<p align="center"> <i> This package is in development. Please see the documentation.</i> </p>
+
+## About
+
+A pure Julia geodesic integration system build on [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) using automatic differentiation (AD) and computer algebra systems (CAS) to efficiently compute the geodesic equation. This package requires only a specification of the non-zero metric components in order to solve the 2nd order geodesic system. Alternatively, an implementation of the four velocity components may be specified to integrate a regular 1st order system.
+
+The motivation behind this package began with an interest in studying reverberation lags around accreting black holes, however the scope has since expanded to facilitate the exploration of generic metrics through time-like, space-like, and null geodesics. 
+
+Our aim is to make testing modified Kerr metrics and alternative gravity theories fast.
+
+Gradus.jl allows for drastically different relativistic simulations to be computed with a composable and reusable API, permitting an end user to simply and expressively calculate physical formulae, create observational signatures, and interface with other popular astrophysics tools. Gradus.jl implements a number of high level abstractions, on the path towards a fully parallelized, high performance numerical relativity ecosystem, scalable from personal computers to super computers.
+
+## Usage
+
+We assume you already have Julia >1.6 installed, in which case you can just add the package from the GitHub URL:
+```julia
+import Pkg;
+Pkg.add("https://github.com/astro-group-bristol/Gradus.jl")
+
+using Gradus
+```
+
+## Credits
+
+This work would not be possible without the incredible work of:
+
+- [The Julia programming language](https://github.com/JuliaLang/Julia)
+- [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl)
+- [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)
+
+<hr>
+
+<p align="center"> Astrophysics Group Bristol </p>

src/AccretionFormulae/AccretionFormulae.jl

@@ -0,0 +1,15 @@
+module AccretionFormulae
+
+import ..Gradus
+import ..GradusBase: AbstractMetricParams, metric
+
+using ..FirstOrderMethods: FirstOrderGeodesicPoint
+using ..Rendering: PointFunction
+
+include("redshift.jl")
+
+const redshift = PointFunction(_redshift_guard)
+
+export redshift
+
+end # module

src/AccretionFormulae/redshift.jl

@@ -0,0 +1,235 @@
+module RedshiftFunctions
+import ..AccretionFormulae.Gradus
+import ..AccretionFormulae.Gradus: __BoyerLindquistFO
+import ..AccretionFormulae: AbstractMetricParams, metric
+using StaticArrays
+using Tullio: @tullio
+
+"""
+    eⱽ(M, r, a, θ)
+Modified from Cunningham et al. (1975) eq. (A2a):
+```math
+e^\\nu = \\sqrt{\\frac{\\Delta \\Sigma}{A}}.
+```
+"""
+eⱽ(M, r, a, θ) = √(
+    __BoyerLindquistFO.Σ(r, a, θ) * __BoyerLindquistFO.Δ(M, r, a) /
+    __BoyerLindquistFO.A(M, r, a, θ),
+)
+
+
+"""
+    eᶲ(M, r, a, θ) 
+Modified from Cunningham et al. (1975) eq. (A2b):
+```math
+e^\\Phi = \\sin \\theta \\sqrt{\\frac{A}{\\Sigma}}.
+```
+"""
+eᶲ(M, r, a, θ) =
+    sin(θ) * √(__BoyerLindquistFO.A(M, r, a, θ) / __BoyerLindquistFO.Σ(r, a, θ))
+
+
+"""
+    ω(M, r, a, θ)
+From Cunningham et al. (1975) eq. (A2c):
+```math
+\\omega = \\frac{2 a M r}{A}.
+```
+"""
+ω(M, r, a, θ) = 2 * a * M * r / __BoyerLindquistFO.A(M, r, a, θ)
+
+
+"""
+    Ωₑ(M, r, a)
+Coordinate angular velocity of an accreting gas. 
+Taken from Cunningham et al. (1975) eq. (A7b):
+```math
+\\Omega_e = \\frac{\\sqrt{M}}{a \\sqrt{M} + r_e^{3/2}}.
+```
+# Notes
+Fanton et al. (1997) use
+```math
+\\Omega_e = \\frac{\\sqrt{M}}{a \\sqrt{M} \\pm r_e^{3/2}},
+```
+where the sign is dependent on co- or contra-rotation. This function may be extended in the future to support this definition.
+"""
+Ωₑ(M, r, a) = √M / (r^1.5 + a * √M)
+
+
+"""
+    Vₑ(M, r, a, θ)  
+Velocity of an accreting gas in a locally non-rotating reference frame (see Bardeen et al. 1973).
+Taken from Cunningham et al. (1975) eq. (A7b):
+```math
+V_e = (\\Omega_e - \\omega) e^{\\Phi - \\nu}.
+```
+"""
+Vₑ(M, r, a, θ) = (Ωₑ(M, r, a) - ω(M, r, a, θ)) * eᶲ(M, r, a, θ) / eⱽ(M, r, a, θ)
+
+
+"""
+    Lₑ(M, rms, a) 
+Angular momentum of an accreting gas within ``r_ms``.
+Taken from Cunningham et al. (1975) eq. (A11b):
+```math
+L_e = \\sqrt{M} \\frac{
+        r_{\\text{ms}}^2 - 2 a \\sqrt{M r_{\\text{ms}}} + a^2
+    }{
+        r_{\\text{ms}}^{3/2} - 2 M \\sqrt{r_{\\text{ms}}} + a \\sqrt{M}
+    }.
+```
+"""
+Lₑ(M, rms, a) = √M * (rms^2 - 2 * a * √(M * rms) + a^2) / (rms^1.5 - 2 * M * √rms + a * √M)
+
+
+"""
+    H(M, rms, r, a)
+Taken from Cunningham et al. (1975) eq. (A12e):
+```math
+H = \\frac{2 M r_e - a \\lambda_e}{\\Delta},
+```
+where we distinguing ``r_e`` as the position of the accreting gas. 
+"""
+H(M, rms, r, a) = (2 * M * r - a * Lₑ(M, rms, a)) / __BoyerLindquistFO.Δ(M, r, a)
+
+
+"""
+    γₑ(M, rms) 
+Taken from Cunningham et al. (1975) eq. (A11c):
+```math
+\\gamma_e = \\sqrt{1 - \\frac{
+        2M
+    }{
+        3 r_{\\text{ms}} 
+    }}.
+```
+"""
+γₑ(M, rms) = √(1 - (2 * M) / (3 * rms))
+
+
+"""
+    uʳ(M, rms, r)
+Taken from Cunningham et al. (1975) eq. (A12b):
+```math
+u^r = - \\sqrt{\\frac{
+        2M
+    }{
+        3 r_{\\text{ms}} 
+    }} \\left(
+        \\frac{ r_{\\text{ms}} }{r_e} - 1
+    \\right)^{3/2}.
+```
+"""
+uʳ(M, rms, r) = -√((2 * M) / (3 * rms)) * (rms / r - 1)^1.5
+
+
+"""
+    uᶲ(M, rms, r, a)
+Taken from Cunningham et al. (1975) eq. (A12c):
+```math
+u^\\phi = \\frac{\\gamma_e}{r_e^2} \\left( 
+        L_e + aH 
+    \\right).
+```
+"""
+uᶲ(M, rms, r, a) = γₑ(M, rms) / r^2 * (Lₑ(M, rms, a) + a * H(M, rms, r, a))
+
+
+"""
+    uᵗ(M, rms, r, a) 
+Taken from Cunningham et al. (1975) eq. (A12b):
+```math
+u^t = \\gamma_e \\left(
+        1 + \\frac{2 M (1 + H)}{r_e}
+    \\right).
+```
+"""
+uᵗ(M, rms, r, a) = γₑ(M, rms) * (1 + 2 * M * (1 + H(M, rms, r, a)) / r)
+
+
+"""
+Experimental API:
+"""
+function regular_pdotu_inv(L, M, r, a, θ)
+    (eⱽ(M, r, a, θ) * √(1 - Vₑ(M, r, a, θ)^2)) / (1 - L * Ωₑ(M, r, a))
+end
+
+@inline function regular_pdotu_inv(m::AbstractMetricParams{T}, u, v) where {T}
+    metric_matrix = metric(m, u)
+    # reverse signs of the velocity vector
+    # since we're integrating backwards
+    p = @inbounds @SVector [-v[1], 0, 0, -v[4]]
+
+    # TODO: this only works for Kerr
+    disc_norm = (eⱽ(m.M, u[2], m.a, u[3]) * √(1 - Vₑ(m.M, u[2], m.a, u[3])^2))
+
+    u_disc = @SVector [1 / disc_norm, 0, 0, Ωₑ(m.M, u[2], m.a) / disc_norm]
+
+    # use Tullio to do the einsum
+    @tullio g := metric_matrix[i, j] * (u_disc[i]) * p[j]
+    1 / g
+end
+
+function plunging_p_dot_u(E, a, M, L, Q, rms, r, sign_r)
+    inv(
+        uᵗ(M, rms, r, a) - uᶲ(M, rms, r, a) * L -
+        sign_r * uʳ(M, rms, r) * __BoyerLindquistFO.Σδr_δλ(E, L, M, Q, r, a) /
+        __BoyerLindquistFO.Δ(M, r, a),
+    )
+end
+
+function plunging_p_dot_u(m::AbstractMetricParams{T}, u, v, rms) where {T}
+    metric_matrix = metric(m, u)
+
+    # reverse signs of the velocity vector
+    # since we're integrating backwards
+    p = @inbounds @SVector [-v[1], v[2], 0, -v[4]]
+
+    u_disc = @inbounds @SVector [
+        uᵗ(m.M, rms, u[2], m.a),
+        uʳ(m.M, rms, u[2]),
+        0,
+        uᶲ(m.M, rms, u[2], m.a),
+    ]
+
+    @tullio g := metric_matrix[i, j] * u_disc[i] * p[j]
+    1 / g
+end
+
+@inline function redshift_function(m::Gradus.BoyerLindquistFO{T}, u, p, λ) where {T}
+    isco = Gradus.isco(m)
+    if u[2] > isco
+        @inbounds regular_pdotu_inv(p.L, m.M, u[2], m.a, u[3])
+    else
+        # change sign if we're after the sign flip
+        # TODO: this isn't robust to multiple sign changes 
+        # TODO: i feel like there should be a better way than checking this with two conditions
+        #       used to have λ > p.changes[1] to make sure we're ahead of the time flip (but when wouldn't we be???)
+        #       now p.changes[1] > 0.0 to make sure there was a time flip at all
+        sign_r = (p.changes[1] > 0.0 ? 1 : -1) * p.r
+        @inbounds plunging_p_dot_u(m.E, m.a, m.M, p.L, p.Q, isco, u[2], sign_r)
+    end
+end
+
+@inline function redshift_function(m::AbstractMetricParams{T}, u, v) where {T}
+    isco = Gradus.isco(m)
+    if u[2] > isco
+        regular_pdotu_inv(m, u, v)
+    else
+        plunging_p_dot_u(m, u, v, isco)
+    end
+end
+
+end # module
+
+# point functions exports
+function _redshift_guard(
+    m::Gradus.BoyerLindquistFO{T},
+    gp::FirstOrderGeodesicPoint{T},
+    max_time,
+) where {T}
+    RedshiftFunctions.redshift_function(m, gp.u, gp.p, gp.t)
+end
+function _redshift_guard(m::AbstractMetricParams{T}, gp, max_time) where {T}
+    RedshiftFunctions.redshift_function(m, gp.u, gp.v)
+end