diff --git a/docs/make.jl b/docs/make.jl
index 94d6f42..95e5f6a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -9,13 +9,14 @@ GENERATED = joinpath(@__DIR__, "..", "examples")
 OUTDIR = joinpath(@__DIR__, "src", "examples")
 
 blacklist = ["gpu", "JuKeBOX", "Krang_logo", "enzyme"]
-SOURCE_FILES = filter(x-> all(i->!occursin(i, x), blacklist), Glob.glob("*.jl", GENERATED))
-foreach(fn -> Literate.markdown(fn, OUTDIR, documenter=true), SOURCE_FILES)
+SOURCE_FILES =
+    filter(x -> all(i -> !occursin(i, x), blacklist), Glob.glob("*.jl", GENERATED))
+foreach(fn -> Literate.markdown(fn, OUTDIR, documenter = true), SOURCE_FILES)
 MD_FILES = [joinpath("examples", file) for file in readdir(OUTDIR)]
 
-Documenter.DocMeta.setdocmeta!(Krang, :DocTestSetup, :(using Krang); recursive=true)
+Documenter.DocMeta.setdocmeta!(Krang, :DocTestSetup, :(using Krang); recursive = true)
 
-format=DocumenterVitepress.MarkdownVitepress(
+format = DocumenterVitepress.MarkdownVitepress(
     repo = "https://github.com/dominic-chang/Krang.jl", # this must be the full URL!
     devbranch = "main",
     devurl = "dev",
@@ -24,35 +25,29 @@ format=DocumenterVitepress.MarkdownVitepress(
 )
 
 makedocs(;
-    sitename="Krang.jl",
-    authors="Dominic <dochang@g.harvard.edu> and contributors",
-    modules=[Krang],
-    repo="https://github.com/dominic-chang/Krang.jl/blob/{commit}{path}#{line}",
-    format=format,
+    sitename = "Krang.jl",
+    authors = "Dominic <dochang@g.harvard.edu> and contributors",
+    modules = [Krang],
+    repo = "https://github.com/dominic-chang/Krang.jl/blob/{commit}{path}#{line}",
+    format = format,
     draft = false,
     source = "src",
     build = "build",
-    pages=[
+    pages = [
         "Home" => "index.md",
         "Getting Started" => "getting_started.md",
         "Examples" => MD_FILES,
         "Theory" => [
-            "Raytracing" => [
-                "kerr_geodesic_summary.md",
-                "time_regularization.md",
-            ],
-            "Polarization" => [
-                "newmann_penrose.md",
-                "polarization.md",
-            ]
+            "Raytracing" => ["kerr_geodesic_summary.md", "time_regularization.md"],
+            "Polarization" => ["newmann_penrose.md", "polarization.md"],
         ],
         "api.md",
     ],
 )
 
 deploydocs(;
-    repo="github.com/dominic-chang/Krang.jl",
+    repo = "github.com/dominic-chang/Krang.jl",
     target = "build", # this is where Vitepress stores its output
-    devbranch="main",
-    push_preview=true,
+    devbranch = "main",
+    push_preview = true,
 )
diff --git a/examples/coordinate-example.jl b/examples/coordinate-example.jl
index 287f34b..46477d1 100644
--- a/examples/coordinate-example.jl
+++ b/examples/coordinate-example.jl
@@ -11,15 +11,13 @@ using CairoMakie
 
 curr_theme = Theme(
     Axis = (
-        xticksvisible = false, 
+        xticksvisible = false,
         xticklabelsvisible = false,
         yticksvisible = false,
         yticklabelsvisible = false,
-        aspect=1
-        ),
-    Heatmap = (
-        rasterize=true,
-    )
+        aspect = 1,
+    ),
+    Heatmap = (rasterize = true,),
 )
 set_theme!(merge!(curr_theme, theme_latexfonts()))
 
@@ -35,7 +33,7 @@ rmax = 10.0; # maximum radius to be ray traced
 ρmax = 15.0; # horizontal and vertical limits of the screen
 
 # Initialize Camera and pre-allocate memory for data to be plotted
-coordinates = (zeros(sze, sze) for _ in 1:3)
+coordinates = (zeros(sze, sze) for _ = 1:3)
 camera = Krang.SlowLightIntensityCamera(metric, θo, -ρmax, ρmax, -ρmax, ρmax, sze);
 colormaps = (:afmhot, :afmhot, :hsv)
 colorrange = ((-20, 20), (0, rmax), (0, 2π))
@@ -46,32 +44,35 @@ colorrange = ((-20, 20), (0, rmax), (0, 2π))
 # coordiante information for a sub-image `n` associated with the ray at the pixel `pixel` which intersects a cone with opening angle `θs`.
 # We will define a function which introduces some basic occlusion effects by checking if the ray has intersected with the cone on 
 # the 'far-side' or the 'near-side' from the obsever.
-function coordinate_point(pix::Krang.AbstractPixel, geometry::Krang.ConeGeometry{T,A}) where {T, A}
+function coordinate_point(
+    pix::Krang.AbstractPixel,
+    geometry::Krang.ConeGeometry{T,A},
+) where {T,A}
     n, rmin, rmax = geometry.attributes
     θs = geometry.opening_angle
 
     coords = ntuple(j -> zero(T), Val(4))
 
-    isindir = false 
-    for _ in 1:2 # Looping over isindir this way is needed to get Metal.jl to work
+    isindir = false
+    for _ = 1:2 # Looping over isindir this way is needed to get Metal.jl to work
         isindir ⊻= true
-        ts, rs, ϕs =  emission_coordinates(pix, geometry.opening_angle, isindir, n)
+        ts, rs, ϕs = emission_coordinates(pix, geometry.opening_angle, isindir, n)
         if rs ≤ rmin || rs ≥ rmax
             continue
         end
-        coords = isnan(rs) ? observation :  (ts, rs, θs, ϕs)
+        coords = isnan(rs) ? observation : (ts, rs, θs, ϕs)
     end
-    return coords 
+    return coords
 end
 
 # ## Drawing coordinate points
 # This function plots the coordinates associated with the n=0,1,2 subimages of a cone with opening angle θs.
 function draw!(axes_list, camera, coordinates, rmin, rmax, θs)
-    times, radii, azimuths = coordinates 
+    times, radii, azimuths = coordinates
     map(axes -> empty!.(axes), axes_list)
 
-    geometries = (Krang.ConeGeometry(θs, (i, rmin, rmax)) for i in 0:2)
-    
+    geometries = (Krang.ConeGeometry(θs, (i, rmin, rmax)) for i = 0:2)
+
     for (i, geometry) in enumerate(geometries)
         rendered_scene = coordinate_point.(camera.screen.pixels, Ref(geometry))
         for I in CartesianIndices(rendered_scene)
@@ -79,25 +80,41 @@ function draw!(axes_list, camera, coordinates, rmin, rmax, θs)
             radii[I] = rendered_scene[I][2]
             azimuths[I] = rendered_scene[I][4]
         end
-        coordinates = (times, radii, mod2pi.(azimuths ))
-        for j in 1:3
-            heatmap!(axes_list[i][j], coordinates[j], colormap = colormaps[j], colorrange=colorrange[j])
+        coordinates = (times, radii, mod2pi.(azimuths))
+        for j = 1:3
+            heatmap!(
+                axes_list[i][j],
+                coordinates[j],
+                colormap = colormaps[j],
+                colorrange = colorrange[j],
+            )
         end
     end
 end
 
 # Create Figure
-fig = Figure(resolution=(700, 700));
+fig = Figure(resolution = (700, 700));
 axes_list = [
     [
-        Axis(fig[i, 1], title=(i==1 ? "Regularized Time" : ""), titlesize=20, ylabel=(i==1 ? L"n=0" : i==2 ? L"n=1" : L"n=2"), ylabelsize=20),
-        Axis(fig[i, 2], title=(i==1 ? "Radius" : ""), titlesize=20),
-        Axis(fig[i, 3], title=(i==1 ? "Azimuth" : ""), titlesize=20),
-    ] for i in 1:3
+        Axis(
+            fig[i, 1],
+            title = (i == 1 ? "Regularized Time" : ""),
+            titlesize = 20,
+            ylabel = (i == 1 ? L"n=0" : i == 2 ? L"n=1" : L"n=2"),
+            ylabelsize = 20,
+        ),
+        Axis(fig[i, 2], title = (i == 1 ? "Radius" : ""), titlesize = 20),
+        Axis(fig[i, 3], title = (i == 1 ? "Azimuth" : ""), titlesize = 20),
+    ] for i = 1:3
 ]
 
 # Create the animation of Cone of Emission Coordinates
-recording = CairoMakie.record(fig, "coordinate.gif", range(0.0, π, length=180), framerate=12) do θs
+recording = CairoMakie.record(
+    fig,
+    "coordinate.gif",
+    range(0.0, π, length = 180),
+    framerate = 12,
+) do θs
     draw!(axes_list, camera, coordinates, rmin, rmax, θs)
 end
 
diff --git a/examples/custom-material-example.jl b/examples/custom-material-example.jl
index 75dde79..13ff177 100644
--- a/examples/custom-material-example.jl
+++ b/examples/custom-material-example.jl
@@ -18,9 +18,12 @@ end
 # This functor will take in a pixel and a geometry and return the redshift associated with a given sub image.
 # You must include the relevant physics in the functor definition. 
 # Here we will include redshift effects associated with a zero angular momentum observer (ZAMO).
-function (m::ZAMORedshifts)(pix::Krang.AbstractPixel{T}, intersection::Krang.Intersection) where {T}
+function (m::ZAMORedshifts)(
+    pix::Krang.AbstractPixel{T},
+    intersection::Krang.Intersection,
+) where {T}
     (; rs, θs, νr, νθ) = intersection
-    α,β=Krang.screen_coordinate(pix)
+    α, β = Krang.screen_coordinate(pix)
 
     ηtemp = η(metric, α, β, θo)
     λtemp = λ(metric, α, θo)
@@ -34,7 +37,7 @@ function (m::ZAMORedshifts)(pix::Krang.AbstractPixel{T}, intersection::Krang.Int
 end
 
 # We need to define the return type or the material for ray tracing. This can be things like a `Float` of a `StokesParams`. Let's use a float for this case.
-function Krang.yield(::ZAMORedshifts{T}) where T
+function Krang.yield(::ZAMORedshifts{T}) where {T}
     return zero(T)
 end
 
@@ -61,23 +64,23 @@ import CairoMakie as CMk
 
 theme = CMk.Theme(
     Axis = (
-        xticksvisible = false, 
+        xticksvisible = false,
         xticklabelsvisible = false,
         yticksvisible = false,
         yticklabelsvisible = false,
-        ),
+    ),
 )
 
 CMk.set_theme!(CMk.merge!(theme, CMk.theme_latexfonts()))
 
-fig = CMk.Figure(resolution=(700, 700));
-ax = CMk.Axis(fig[1, 1], title="Redshifts", titlesize=20, aspect=1)
-hm = CMk.heatmap!(ax, redshifts, colormap=:afmhot)
-CMk.Colorbar(fig[1, 2], hm, label="Redshifts", labelsize=20)
+fig = CMk.Figure(resolution = (700, 700));
+ax = CMk.Axis(fig[1, 1], title = "Redshifts", titlesize = 20, aspect = 1)
+hm = CMk.heatmap!(ax, redshifts, colormap = :afmhot)
+CMk.Colorbar(fig[1, 2], hm, label = "Redshifts", labelsize = 20)
 CMk.save("redshifts.png", fig)
 
 # ![redshifts](redshifts.png)
 
 # Saving the redshifts to a .npy file
 using NPZ
-npzwrite("redshifts_n1.npy", redshifts)
\ No newline at end of file
+npzwrite("redshifts_n1.npy", redshifts)
diff --git a/examples/mino-time-example.jl b/examples/mino-time-example.jl
index 5a56ec1..126010b 100644
--- a/examples/mino-time-example.jl
+++ b/examples/mino-time-example.jl
@@ -11,17 +11,17 @@ import GLMakie as GLMk
 GLMk.Makie.inline!(true)
 
 curr_theme = GLMk.Theme(# Makie theme
-    fontsize=20,
-    Axis=(
-        xticksvisible=false,
-        xticklabelsvisible=false,
-        yticksvisible=false,
-        yticklabelsvisible=false,
-        leftspinevisible=false,
-        rightspinevisible=false,
-        topspinevisible=false,
-        bottomspinevisible=false,
-        titlefontsize=30,
+    fontsize = 20,
+    Axis = (
+        xticksvisible = false,
+        xticklabelsvisible = false,
+        yticksvisible = false,
+        yticklabelsvisible = false,
+        leftspinevisible = false,
+        rightspinevisible = false,
+        topspinevisible = false,
+        bottomspinevisible = false,
+        titlefontsize = 30,
     ),
 )
 
@@ -45,41 +45,56 @@ camera = Krang.SlowLightIntensityCamera(metric, θo, -ρmax, ρmax, -ρmax, ρma
 # We will create a loop to plot the emission coordinates for each `τ` using the `emission_coordinates!` function.
 # Let us now create a figure to plot the emission coordinates on.
 
-fig = GLMk.Figure(resolution=(500, 600));
+fig = GLMk.Figure(resolution = (500, 600));
 
 
-recording = GLMk.record(fig, "raytrace.gif", range(0.1, 3, length=290), framerate=15) do τ
-    GLMk.empty!(fig)
+recording =
+    GLMk.record(fig, "raytrace.gif", range(0.1, 3, length = 290), framerate = 15) do τ
+        GLMk.empty!(fig)
 
-    coordinates = zeros(4, size(camera.screen.pixels)...) # Pre allocated array to store coordinates
-    emission_coordinates!(coordinates, camera, τ)
-    time = coordinates[1,:,:]
-    radius = coordinates[2,:,:]
-    inclination = coordinates[3,:,:]
-    azimuth = mod2pi.(coordinates[4,:,:])
+        coordinates = zeros(4, size(camera.screen.pixels)...) # Pre allocated array to store coordinates
+        emission_coordinates!(coordinates, camera, τ)
+        time = coordinates[1, :, :]
+        radius = coordinates[2, :, :]
+        inclination = coordinates[3, :, :]
+        azimuth = mod2pi.(coordinates[4, :, :])
 
-    data = (time, radius, inclination, azimuth)
-    titles = (GLMk.L"\text{Regularized Time }(t_s)", GLMk.L"\text{Radius }(r_s)", GLMk.L"\text{Inclination }(\theta_s)", GLMk.L"\text{Azimuth } (\phi_s)")
-    colormaps = (:afmhot, :afmhot, :afmhot, :hsv)
-    colorrange = ((-20, 20), (0, 10.0), (0,π), (0, 2π))
-    indices = ((1,1), ())
-
-    for i in 1:4
-        hm = GLMk.heatmap!(
-            GLMk.Axis(getindex(fig, (i > 2 ? 2 : 1), (iszero(i%2) ? 3 : 1)); aspect=1, title=titles[i]),
-            data[i],
-            colormap=colormaps[i],
-            colorrange=colorrange[i]
+        data = (time, radius, inclination, azimuth)
+        titles = (
+            GLMk.L"\text{Regularized Time }(t_s)",
+            GLMk.L"\text{Radius }(r_s)",
+            GLMk.L"\text{Inclination }(\theta_s)",
+            GLMk.L"\text{Azimuth } (\phi_s)",
         )
-        cb = GLMk.Colorbar(fig[(i > 2 ? 2 : 1), (iszero(i%2) ? 3 : 1)+1], hm; labelsize=30, ticklabelsize=20)
+        colormaps = (:afmhot, :afmhot, :afmhot, :hsv)
+        colorrange = ((-20, 20), (0, 10.0), (0, π), (0, 2π))
+        indices = ((1, 1), ())
+
+        for i = 1:4
+            hm = GLMk.heatmap!(
+                GLMk.Axis(
+                    getindex(fig, (i > 2 ? 2 : 1), (iszero(i % 2) ? 3 : 1));
+                    aspect = 1,
+                    title = titles[i],
+                ),
+                data[i],
+                colormap = colormaps[i],
+                colorrange = colorrange[i],
+            )
+            cb = GLMk.Colorbar(
+                fig[(i > 2 ? 2 : 1), (iszero(i % 2) ? 3 : 1)+1],
+                hm;
+                labelsize = 30,
+                ticklabelsize = 20,
+            )
+        end
+
+        ax = GLMk.Axis(fig[3, 1:3], height = 60)
+        GLMk.hidedecorations!(ax)
+        GLMk.text!(ax, 0, 100; text = GLMk.L"θ_o=%$(Int(floor(θo*180/π)))^\circ")
+        GLMk.rowgap!(fig.layout, 1, GLMk.Fixed(0))
     end
 
-    ax = GLMk.Axis(fig[3, 1:3], height=60)
-    GLMk.hidedecorations!(ax)
-    GLMk.text!(ax,0,100; text=GLMk.L"θ_o=%$(Int(floor(θo*180/π)))^\circ")
-    GLMk.rowgap!(fig.layout, 1, GLMk.Fixed(0))
-end
-
 # ![image](raytrace.gif)
 
 # ## Plotting rays
@@ -90,26 +105,26 @@ end
 camera = Krang.SlowLightIntensityCamera(metric, θo, -3, 3, -3, 3, 4);
 
 fig = GLMk.Figure()
-ax = GLMk.Axis3(fig[1,1], aspect=(1,1,1))
-GLMk.xlims!(ax, (-3, 3)) 
-GLMk.ylims!(ax, (-3, 3)) 
-GLMk.zlims!(ax, (-3, 3)) 
+ax = GLMk.Axis3(fig[1, 1], aspect = (1, 1, 1))
+GLMk.xlims!(ax, (-3, 3))
+GLMk.ylims!(ax, (-3, 3))
+GLMk.zlims!(ax, (-3, 3))
 lines_to_plot = []
 lines_to_plot = Krang.generate_ray.(camera.screen.pixels, 5_000)
 
-sphere = GLMk.Sphere(GLMk.Point(0.0,0.0,0.0), horizon(metric))
-GLMk.mesh!(ax, sphere, color=:black) # Sphere to represent black hole
+sphere = GLMk.Sphere(GLMk.Point(0.0, 0.0, 0.0), horizon(metric))
+GLMk.mesh!(ax, sphere, color = :black) # Sphere to represent black hole
 
-for i in lines_to_plot; 
-    ray = map(x-> begin
-    (;rs, θs, ϕs) = x
-    [rs*sin(θs)*cos(ϕs), rs*sin(θs)*sin(ϕs), rs*cos(θs)]
+for i in lines_to_plot
+    ray = map(x -> begin
+        (; rs, θs, ϕs) = x
+        [rs * sin(θs) * cos(ϕs), rs * sin(θs) * sin(ϕs), rs * cos(θs)]
     end, i)
     ray = hcat(ray...)
-    GLMk.lines!(ax, ray) 
+    GLMk.lines!(ax, ray)
 end
 fig
 
 GLMk.save("mino_time_rays.png", fig)
 
-# ![Photons trajectories around Kerr black hole in Boyer-Lindquist Coordinates](mino_time_rays.png)
\ No newline at end of file
+# ![Photons trajectories around Kerr black hole in Boyer-Lindquist Coordinates](mino_time_rays.png)
diff --git a/examples/neural-net-example.jl b/examples/neural-net-example.jl
index ec8922c..b609ec8 100644
--- a/examples/neural-net-example.jl
+++ b/examples/neural-net-example.jl
@@ -18,7 +18,7 @@ rng = Random.GLOBAL_RNG
 # Lets define an `ImageModel` which will be comprised of an emission layer that we will raytrace.
 # We will do this by first creating a struct to represent our image model that will store our emission model as a layer.
 
-struct ImageModel{T <: Chain}
+struct ImageModel{T<:Chain}
     emission_layer::T
 end
 
@@ -27,44 +27,45 @@ end
 # We will assume that the emission is coming from emission that originates in the equatorial plane.
 function (m::ImageModel)(x, ps, st)
     metric = Krang.Kerr(0.99e0)
-    θo = Float64(20/180*π)
-    pixels = Krang.IntensityPixel.(Ref(metric), x[1,:], x[2,:], θo)
+    θo = Float64(20 / 180 * π)
+    pixels = Krang.IntensityPixel.(Ref(metric), x[1, :], x[2, :], θo)
 
     sze = unsafe_trunc(Int, sqrt(size(x)[2]))
-    coords = zeros(Float64, 2,sze*sze)
-    emission_vals = zeros(Float64, 1, sze*sze)
-    for n in 0:1
-        for i in 1:sze
-            Threads.@threads for j in 1:sze
+    coords = zeros(Float64, 2, sze * sze)
+    emission_vals = zeros(Float64, 1, sze * sze)
+    for n = 0:1
+        for i = 1:sze
+            Threads.@threads for j = 1:sze
                 pix = pixels[i+(j-1)*sze]
                 α, β = Krang.screen_coordinate(pix)
                 T = typeof(α)
-                rs, ϕs = Krang.emission_coordinates_fast_light(pix, Float64(π/2), β > 0, n)[1:3]
+                rs, ϕs =
+                    Krang.emission_coordinates_fast_light(pix, Float64(π / 2), β > 0, n)[1:3]
                 xs = rs * cos(ϕs)
                 ys = rs * sin(ϕs)
                 if hypot(xs, ys) ≤ Krang.horizon(metric)
-                    coords[1,i+(j-1)*sze] = zero(T)
-                    coords[2,i+(j-1)*sze] = zero(T)
+                    coords[1, i+(j-1)*sze] = zero(T)
+                    coords[2, i+(j-1)*sze] = zero(T)
                 else
-                    coords[1,i+(j-1)*sze] = xs
-                    coords[2,i+(j-1)*sze] = ys
+                    coords[1, i+(j-1)*sze] = xs
+                    coords[2, i+(j-1)*sze] = ys
                 end
 
             end
         end
         emission_vals .+= m.emission_layer(coords, ps, st)[1]
     end
-    emission_vals,st 
+    emission_vals, st
 end
 
 # Lets define an emisison layer for our model as a simple fully connected neural network with 2 hidden layers.
 # The emission layer will take in 2D coordinates on an equatorial disk in the bulk spacetime and return a scalar intensity value.
 
 emission_model = Chain(
-    Dense(2 => 20, Lux.sigmoid),  
+    Dense(2 => 20, Lux.sigmoid),
     Dense(20 => 20, Lux.sigmoid),
-    Dense(20 => 1, Lux.sigmoid)
-    ) 
+    Dense(20 => 1, Lux.sigmoid),
+)
 
 ps, st = Lux.setup(rng, emission_model); # Get the emission model parameters and state
 
@@ -74,22 +75,27 @@ image_model = ImageModel(emission_model)
 ## Plotting the model
 
 # Lets create an 20x20 pixel image of the `image_model` with a field of view of $10 MG/c^2$.
- 
+
 sze = 20
 ρmax = 10e0
-pixels = zeros(Float64, 2, sze*sze)
+pixels = zeros(Float64, 2, sze * sze)
 for (iiter, i) in enumerate(range(-ρmax, ρmax, sze))
     for (jiter, j) in enumerate(range(-ρmax, ρmax, sze))
-        pixels[1,iiter+(jiter-1)*sze] = Float64(i)
-        pixels[2,iiter+(jiter-1)*sze] = Float64(j)
+        pixels[1, iiter+(jiter-1)*sze] = Float64(i)
+        pixels[2, iiter+(jiter-1)*sze] = Float64(j)
     end
 end
 
 # We can see the effects of raytracing on emission in the bulk spacetime by plotting an image of the emission model and the image model.
 using CairoMakie
 curr_theme = Theme(
-    Axis = (xticksvisible=false, xticklabelsvisible=false, yticksvisible=false, yticklabelsvisible=false,),
-    Heatmap =(colormap=:afmhot, ),
+    Axis = (
+        xticksvisible = false,
+        xticklabelsvisible = false,
+        yticksvisible = false,
+        yticklabelsvisible = false,
+    ),
+    Heatmap = (colormap = :afmhot,),
 )
 set_theme!(merge(curr_theme, theme_latexfonts()))
 
@@ -97,8 +103,11 @@ emitted_intensity = reshape(emission_model(pixels, ps, st)[1], sze, sze)
 received_intensity = reshape(image_model(pixels, ps, st)[1], sze, sze)
 
 fig = Figure();
-heatmap!(Axis(fig[1,1], aspect=1, title="Emission Model"), emitted_intensity)
-heatmap!(Axis(fig[1,2], aspect=1, title="Image Model (Lensed Emission Model)"), received_intensity)
+heatmap!(Axis(fig[1, 1], aspect = 1, title = "Emission Model"), emitted_intensity)
+heatmap!(
+    Axis(fig[1, 2], aspect = 1, title = "Image Model (Lensed Emission Model)"),
+    received_intensity,
+)
 save("emission_model_and_target_model.png", fig)
 
 # ![image](emission_model_and_target_model.png)
@@ -106,18 +115,18 @@ save("emission_model_and_target_model.png", fig)
 # ## Fitting the NeRF model
 # This will be a toy example showing the mechanics of fitting our ImageModel to a target image using the normalized cross correlation as a kernel for our loss function.
 # This will be the image we will try to fit our model to.
-target_img = reshape(received_intensity, 1, sze*sze);
+target_img = reshape(received_intensity, 1, sze * sze);
 
 # Lets fit our model using the normalized cross correlation as a kernel for our loss function.
 
 using Enzyme
-using Optimization 
+using Optimization
 using OptimizationOptimisers
 using StatsBase
 using ComponentArrays
 
-function mse(img1::Matrix{T}, img2::Matrix{T}) where T
-    mean(((img1 ./ sum(img1))  .- (img2 ./ sum(img2))) .^ 2)
+function mse(img1::Matrix{T}, img2::Matrix{T}) where {T}
+    mean(((img1 ./ sum(img1)) .- (img2 ./ sum(img2))) .^ 2)
 end
 
 function loss_function(pixels, y, ps, st)
@@ -135,21 +144,29 @@ received_intensity = reshape(image_model(pixels, ps, st)[1], sze, sze)
 loss_function(pixels, target_img, ps, st)
 
 fig = Figure();
-heatmap!(Axis(fig[1,1], aspect=1, title="Emission Model"), emitted_intensity, colormap=:afmhot)
-heatmap!(Axis(fig[1,2], aspect=1, title="Imgage Model (Lensed Emission Model)"), received_intensity, colormap=:afmhot)
+heatmap!(
+    Axis(fig[1, 1], aspect = 1, title = "Emission Model"),
+    emitted_intensity,
+    colormap = :afmhot,
+)
+heatmap!(
+    Axis(fig[1, 2], aspect = 1, title = "Imgage Model (Lensed Emission Model)"),
+    received_intensity,
+    colormap = :afmhot,
+)
 save("emission_model_and_image_model.png", fig)
 
 # ![image](emission_model_and_image_model.png)
 
 # Lets define callback function to print the loss as we optimize our model.
-mutable struct Callback 
+mutable struct Callback
     counter::Int
     stride::Int
     const f::Function
 end
 Callback(stride, f) = Callback(0, stride, f)
 function (c::Callback)(state, loss, others...)
-    c.counter += 1  
+    c.counter += 1
     if c.counter % c.stride == 0
         @info "On step $(c.counter) loss = $(loss)"
         return false
@@ -159,40 +176,58 @@ function (c::Callback)(state, loss, others...)
 end
 
 # We can now optimize our model using the ADAM optimizer.
-ps_trained, st_trained = let st=Ref(st), x=pixels, y=reshape(target_img, 1, sze*sze)
-    
+ps_trained, st_trained = let st = Ref(st), x = pixels, y = reshape(target_img, 1, sze * sze)
+
     optprob = Optimization.OptimizationProblem(
         Optimization.OptimizationFunction(
-            function(ps, constants)
+            function (ps, constants)
                 loss, st[] = loss_function(x, y, ps, st[])
                 loss
             end,
-            Optimization.AutoEnzyme()
+            Optimization.AutoEnzyme(),
         ),
-        ComponentArrays.ComponentVector{Float64}(ps)
+        ComponentArrays.ComponentVector{Float64}(ps),
     )
-    
+
     solution = Optimization.solve(
         optprob,
         OptimizationOptimisers.Adam(),
-        maxiters = 5_000, 
-    callback=Callback(100,()->nothing)
+        maxiters = 5_000,
+        callback = Callback(100, () -> nothing),
     )
-    
+
     solution.u, st[]
 end
 
 # Let's plot the results of our optimization. and compare it to the target image.
-received_intensity, st = ((x) -> (reshape(x[1], sze, sze), x[2]))(image_model(pixels, ps_trained, st_trained))
+received_intensity, st =
+    ((x) -> (reshape(x[1], sze, sze), x[2]))(image_model(pixels, ps_trained, st_trained))
 acc_intensity, st = ((x) -> (reshape(x[1], sze, sze), x[2]))(image_model(pixels, ps, st))
 loss_function(pixels, target_img, ps, st)
 loss_function(pixels, target_img, ps_trained, st_trained)
 using Printf
-begin 
-    fig = Figure(size=(700, 300));
-    heatmap!(Axis(fig[1,1], aspect=1, title="Target Image"), reshape(target_img, sze, sze))
-    heatmap!(Axis(fig[1,2], aspect=1, title="Starting State (loss=$(@sprintf("%0.2e", loss_function(pixels, target_img, ps, st)[1])))"), acc_intensity)  
-    heatmap!(Axis(fig[1,3], aspect=1, title="Fitted State (loss=$(@sprintf("%0.2e", loss_function(pixels, target_img, ps_trained, st_trained)[1])))"), received_intensity)
+begin
+    fig = Figure(size = (700, 300))
+    heatmap!(
+        Axis(fig[1, 1], aspect = 1, title = "Target Image"),
+        reshape(target_img, sze, sze),
+    )
+    heatmap!(
+        Axis(
+            fig[1, 2],
+            aspect = 1,
+            title = "Starting State (loss=$(@sprintf("%0.2e", loss_function(pixels, target_img, ps, st)[1])))",
+        ),
+        acc_intensity,
+    )
+    heatmap!(
+        Axis(
+            fig[1, 3],
+            aspect = 1,
+            title = "Fitted State (loss=$(@sprintf("%0.2e", loss_function(pixels, target_img, ps_trained, st_trained)[1])))",
+        ),
+        received_intensity,
+    )
     save("neural_net_results.png", fig)
 end
 
diff --git a/examples/polarization-example.jl b/examples/polarization-example.jl
index 5bfd39a..6f8c2de 100644
--- a/examples/polarization-example.jl
+++ b/examples/polarization-example.jl
@@ -7,15 +7,15 @@
 #
 # First, let's import Krang and CairoMakie for plotting.
 using Krang
-using CairoMakie 
+using CairoMakie
 
 curr_theme = Theme(
     Axis = (
-        xticksvisible = false, 
+        xticksvisible = false,
         xticklabelsvisible = false,
         yticksvisible = false,
         yticklabelsvisible = false,
-        ),
+    ),
 )
 set_theme!(merge!(curr_theme, theme_latexfonts()))
 
@@ -38,13 +38,13 @@ camera = Krang.IntensityCamera(metric, θo, -ρmax, ρmax, -ρmax, ρmax, 400);
 βv = 0.8776461626924748
 σ = 0.7335172899224874
 η1 = 2.6444786738735804
-η2 = π-η1
+η2 = π - η1
 
 # These will be used to define the magnetic field and fluid velocity.
 
-magfield1 = Krang.SVector(sin(ι)*cos(η1), sin(ι)*sin(η1), cos(ι));
-magfield2 = Krang.SVector(sin(ι)*cos(η2), sin(ι)*sin(η2), cos(ι));
-vel = Krang.SVector(βv, (π/2), χ);
+magfield1 = Krang.SVector(sin(ι) * cos(η1), sin(ι) * sin(η1), cos(ι));
+magfield2 = Krang.SVector(sin(ι) * cos(η2), sin(ι) * sin(η2), cos(ι));
+vel = Krang.SVector(βv, (π / 2), χ);
 R = 3.3266154761905455
 p1 = 4.05269835622511
 p2 = 4.411852974336667
@@ -54,10 +54,26 @@ p2 = 4.411852974336667
 # This includes the sub-images to ray trace, which in our case will be the n=0 and n=1 sub-images.
 
 θs = (75 * π / 180)
-material1 = Krang.ElectronSynchrotronPowerLawPolarization(magfield1..., vel..., σ, R, p1, p2, (0,1,));
+material1 = Krang.ElectronSynchrotronPowerLawPolarization(
+    magfield1...,
+    vel...,
+    σ,
+    R,
+    p1,
+    p2,
+    (0, 1),
+);
 geometry1 = Krang.ConeGeometry((θs))
-material2 = Krang.ElectronSynchrotronPowerLawPolarization(magfield2..., vel..., σ, R, p1, p2, (0,1,));
-geometry2 = Krang.ConeGeometry((π-θs))
+material2 = Krang.ElectronSynchrotronPowerLawPolarization(
+    magfield2...,
+    vel...,
+    σ,
+    R,
+    p1,
+    p2,
+    (0, 1),
+);
+geometry2 = Krang.ConeGeometry((π - θs))
 
 
 # We will create two meshes, one for each geometry anc create a scene with both meshes.
@@ -69,16 +85,25 @@ mesh2 = Krang.Mesh(geometry2, material2)
 scene = Krang.Scene((mesh1, mesh2))
 stokesvals = render(camera, scene)
 
-fig = Figure(resolution=(700, 700));
-ax1 = Axis(fig[1, 1], aspect=1, title="I")
-ax2 = Axis(fig[1, 2], aspect=1, title="Q")
-ax3 = Axis(fig[2, 1], aspect=1, title="U")
-ax4 = Axis(fig[2, 2], aspect=1, title="V")
+fig = Figure(resolution = (700, 700));
+ax1 = Axis(fig[1, 1], aspect = 1, title = "I")
+ax2 = Axis(fig[1, 2], aspect = 1, title = "Q")
+ax3 = Axis(fig[2, 1], aspect = 1, title = "U")
+ax4 = Axis(fig[2, 2], aspect = 1, title = "V")
 colormaps = [:afmhot, :redsblues, :redsblues, :redsblues]
 
-zip([ax1, ax2, ax3, ax4], [getproperty.(stokesvals,:I), getproperty.(stokesvals,:Q), getproperty.(stokesvals,:U), getproperty.(stokesvals,:V)], colormaps) .|> x->heatmap!(x[1], x[2], colormap=x[3])
+zip(
+    [ax1, ax2, ax3, ax4],
+    [
+        getproperty.(stokesvals, :I),
+        getproperty.(stokesvals, :Q),
+        getproperty.(stokesvals, :U),
+        getproperty.(stokesvals, :V),
+    ],
+    colormaps,
+) .|> x -> heatmap!(x[1], x[2], colormap = x[3])
 fig
 
 save("polarization_example.png", fig)
 
-# ![polarization of emission model](polarization_example.png)
\ No newline at end of file
+# ![polarization of emission model](polarization_example.png)
diff --git a/examples/raytracing-mesh-example.jl b/examples/raytracing-mesh-example.jl
index 9bea3e2..5df1451 100644
--- a/examples/raytracing-mesh-example.jl
+++ b/examples/raytracing-mesh-example.jl
@@ -5,7 +5,7 @@ using GeometryBasics
 using FileIO
 
 metric = Krang.Kerr(0.99) # Kerr metric with a spin of 0.99
-θo = 90/180*π # Inclination angle of the observer
+θo = 90 / 180 * π # Inclination angle of the observer
 ρmax = 12.0 # Horizontal and Vertical extent of the screen
 sze = 100 # Resolution of the screen is sze x sze
 
@@ -15,9 +15,22 @@ camera = Krang.SlowLightIntensityCamera(metric, θo, -ρmax, ρmax, -ρmax, ρma
 bunny_mesh = translate(
     rotate(
         scale(
-            load(download("https://graphics.stanford.edu/~mdfisher/Data/Meshes/bunny.obj", "bunny.obj")), 100
-        ),π/2, 1.0, 0.0, 0.0
-    ), 2.0,7.0,-10.0
+            load(
+                download(
+                    "https://graphics.stanford.edu/~mdfisher/Data/Meshes/bunny.obj",
+                    "bunny.obj",
+                ),
+            ),
+            100,
+        ),
+        π / 2,
+        1.0,
+        0.0,
+        0.0,
+    ),
+    2.0,
+    7.0,
+    -10.0,
 );
 
 # Ray trace the mesh embedded in the Kerr space time. The emission picked up by a ray will be the sum of all the intersections of the ray with the mesh
@@ -31,8 +44,8 @@ intersections = raytrace(camera, bunny_mesh);
 # And plot the image with GLMakie, 
 
 fig = GLMk.Figure()
-ax = GLMk.Axis(fig[1,1], aspect=1)  
-GLMk.heatmap!(ax, intersections, colormap=:afmhot);
+ax = GLMk.Axis(fig[1, 1], aspect = 1)
+GLMk.heatmap!(ax, intersections, colormap = :afmhot);
 
 save("mesh_geometry_example.png", fig);
 
@@ -41,26 +54,38 @@ save("mesh_geometry_example.png", fig);
 # Below is a short plotting routine that generates a video showing the scene from three different perspectives, and live renders the image.
 
 fig = GLMk.Figure()
-ax1 = GLMk.Axis3(fig[1,1], aspect=(1,1,1), title="Top view of scene", elevation = π/2, azimuth = π)
-ax2 = GLMk.Axis3(fig[1,2], aspect=(1,1,1), title="Side view of scene", elevation = 0.0, azimuth = 3π/2)
-ax3 = GLMk.Axis3(fig[2,1], aspect=(1,1,1), title="Isometric view of scene")
-
-ax = GLMk.Axis(fig[2,2], aspect=1, title="Heatmap of ray intersections with mesh")
+ax1 = GLMk.Axis3(
+    fig[1, 1],
+    aspect = (1, 1, 1),
+    title = "Top view of scene",
+    elevation = π / 2,
+    azimuth = π,
+)
+ax2 = GLMk.Axis3(
+    fig[1, 2],
+    aspect = (1, 1, 1),
+    title = "Side view of scene",
+    elevation = 0.0,
+    azimuth = 3π / 2,
+)
+ax3 = GLMk.Axis3(fig[2, 1], aspect = (1, 1, 1), title = "Isometric view of scene")
+
+ax = GLMk.Axis(fig[2, 2], aspect = 1, title = "Heatmap of ray intersections with mesh")
 
 for a in [ax1, ax2, ax3]
-    GLMk.xlims!(a, (-10, 10)) 
-    GLMk.ylims!(a, (-10, 10)) 
-    GLMk.zlims!(a, (-10, 10)) 
+    GLMk.xlims!(a, (-10, 10))
+    GLMk.ylims!(a, (-10, 10))
+    GLMk.zlims!(a, (-10, 10))
     GLMk.hidedecorations!(a)
 end
 
 GLMk.hidedecorations!(ax)
-sphere = GLMk.Sphere(GLMk.Point(0.0,0.0,0.0), horizon(metric)) # Sphere to represent black hole
+sphere = GLMk.Sphere(GLMk.Point(0.0, 0.0, 0.0), horizon(metric)) # Sphere to represent black hole
 lines_to_plot = Krang.generate_ray.(camera.screen.pixels, 100) # 100 is the number of steps to take along the ray
 
 img = zeros(sze, sze)
-recording = GLMk.record(fig, "mesh.mp4", 1:sze*sze, framerate=120) do i
-    line = lines_to_plot[i] 
+recording = GLMk.record(fig, "mesh.mp4", 1:sze*sze, framerate = 120) do i
+    line = lines_to_plot[i]
 
     img[i] = intersections[i]
 
@@ -70,15 +95,26 @@ recording = GLMk.record(fig, "mesh.mp4", 1:sze*sze, framerate=120) do i
     GLMk.empty!(ax3)
 
     for a in [ax1, ax2, ax3]
-        GLMk.mesh!(a, bunny_mesh, color=[parse(GLMk.Colorant, "rgba(0%, 50%, 70%,1.0)") for tri in bunny_mesh.position], colormap = :blues, transparency=true)
-        GLMk.mesh!(a, sphere, color=:black) 
+        GLMk.mesh!(
+            a,
+            bunny_mesh,
+            color = [
+                parse(GLMk.Colorant, "rgba(0%, 50%, 70%,1.0)") for
+                tri in bunny_mesh.position
+            ],
+            colormap = :blues,
+            transparency = true,
+        )
+        GLMk.mesh!(a, sphere, color = :black)
     end
 
-    cart_line = map(x->(x.rs*sin(x.θs)*cos(x.ϕs), x.rs*sin(x.θs)*sin(x.ϕs), x.rs*cos(x.θs)), line)
-    GLMk.lines!(ax3, cart_line, color=:red)
-    GLMk.heatmap!(ax, img, colormap=:afmhot, colorrange=(0, 8))
+    cart_line = map(
+        x -> (x.rs * sin(x.θs) * cos(x.ϕs), x.rs * sin(x.θs) * sin(x.ϕs), x.rs * cos(x.θs)),
+        line,
+    )
+    GLMk.lines!(ax3, cart_line, color = :red)
+    GLMk.heatmap!(ax, img, colormap = :afmhot, colorrange = (0, 8))
 end
 # ```@raw html 
 # <video loop muted playsinline controls src="./mesh.mp4" />
 # ```
-
diff --git a/src/Krang.jl b/src/Krang.jl
index 0f8a12d..c529a0c 100644
--- a/src/Krang.jl
+++ b/src/Krang.jl
@@ -7,11 +7,10 @@ import GeometryBasics
 using Rotations
 using Roots
 using KernelAbstractions
-@template (FUNCTIONS, METHODS, MACROS) =
-    """
-    $(TYPEDSIGNATURES)
-    $(DOCSTRING)
-    """
+@template (FUNCTIONS, METHODS, MACROS) = """
+                                         $(TYPEDSIGNATURES)
+                                         $(DOCSTRING)
+                                         """
 
 include("metrics/AbstractMetric.jl")
 include("metrics/Kerr/Kerr.jl")
diff --git a/src/cameras/IntensityCamera.jl b/src/cameras/IntensityCamera.jl
index f67dfa7..cb6142b 100644
--- a/src/cameras/IntensityCamera.jl
+++ b/src/cameras/IntensityCamera.jl
@@ -7,7 +7,7 @@ Each Pixel is associated with a single ray, and caches some information about th
 """
 struct IntensityPixel{T} <: AbstractPixel{T}
     metric::Kerr{T}
-    screen_coordinate::NTuple{2, T}
+    screen_coordinate::NTuple{2,T}
     "Radial roots"
     roots::NTuple{4,Complex{T}}
     "Radial antiderivative"
@@ -50,12 +50,14 @@ struct IntensityPixel{T} <: AbstractPixel{T}
         I0_inf = Krang.Ir_inf(met, roots)
         new{T}(
             met,
-            (α, β), 
+            (α, β),
             roots,
-            I0_inf, 
+            I0_inf,
             total_mino_time(met, roots),
-            Krang._absGθo_Gθhat(met, θo, tempη, tempλ), 
-            θo, tempη, tempλ
+            Krang._absGθo_Gθhat(met, θo, tempη, tempλ),
+            θo,
+            tempη,
+            tempλ,
         )
     end
 end
@@ -65,20 +67,29 @@ end
 
 Screen made of `IntensityPixel`s.
 """
-struct IntensityScreen{T, A <:AbstractMatrix} <: AbstractScreen
+struct IntensityScreen{T,A<:AbstractMatrix} <: AbstractScreen
     "Minimum and Maximum Bardeen α values"
-    αrange::NTuple{2, T}
+    αrange::NTuple{2,T}
 
     "Minimum and Maximum Bardeen β values"
-    βrange::NTuple{2, T}
+    βrange::NTuple{2,T}
 
     "Data type that stores screen pixel information"
     pixels::A
 
-    @kernel function _generate_screen!(screen, met::Kerr{T}, αmin, αmax, βmin, βmax, θo, res) where T
-        I,J = @index(Global, NTuple)
-        α = αmin + (αmax - αmin) * (T(I)-1) / (res-1)
-        β = βmin + (βmax - βmin) * (T(J)-1) / (res-1)
+    @kernel function _generate_screen!(
+        screen,
+        met::Kerr{T},
+        αmin,
+        αmax,
+        βmin,
+        βmax,
+        θo,
+        res,
+    ) where {T}
+        I, J = @index(Global, NTuple)
+        α = αmin + (αmax - αmin) * (T(I) - 1) / (res - 1)
+        β = βmin + (βmax - βmin) * (T(J) - 1) / (res - 1)
         screen[I, J] = IntensityPixel(met, α, β, θo)
     end
 
@@ -101,14 +112,33 @@ struct IntensityScreen{T, A <:AbstractMatrix} <: AbstractScreen
     # Returns
     - `IntensityScreen{T, A}`: An intensity screen object.
     """
-    function IntensityScreen(met::Kerr{T}, αmin::T, αmax::T, βmin::T, βmax::T, θo::T, res; A=Matrix) where {T}
+    function IntensityScreen(
+        met::Kerr{T},
+        αmin::T,
+        αmax::T,
+        βmin::T,
+        βmax::T,
+        θo::T,
+        res;
+        A = Matrix,
+    ) where {T}
         screen = A(Matrix{IntensityPixel{T}}(undef, res, res))
 
         backend = get_backend(screen)
 
-        _generate_screen!(backend)(screen, met, αmin, αmax, βmin, βmax, θo, res, ndrange = (res, res))
-        
-        new{T, A}((αmin, αmax), (βmin, βmax), screen)
+        _generate_screen!(backend)(
+            screen,
+            met,
+            αmin,
+            αmax,
+            βmin,
+            βmax,
+            θo,
+            res,
+            ndrange = (res, res),
+        )
+
+        new{T,A}((αmin, αmax), (βmin, βmax), screen)
     end
 end
 
@@ -118,12 +148,12 @@ end
 Camera that caches fast light raytracing information for an observer sitting at radial infinity.
 The frame of this observer is alligned with the Boyer-Lindquist frame.
 """
-struct IntensityCamera{T, A} <: AbstractCamera
+struct IntensityCamera{T,A} <: AbstractCamera
     metric::Kerr{T}
     "Data type that stores screen pixel information"
-    screen::IntensityScreen{T, A}
+    screen::IntensityScreen{T,A}
     "Observer screen_coordinate"
-    screen_coordinate::NTuple{2, T}
+    screen_coordinate::NTuple{2,T}
 
     @doc """
         IntensityCamera(met::Kerr{T}, θo, αmin, αmax, βmin, βmax, res::Int; A=Matrix) where {T}
@@ -143,17 +173,48 @@ struct IntensityCamera{T, A} <: AbstractCamera
     # Returns
     - `IntensityCamera{T, A}`: An intensity camera object.
     """
-    function IntensityCamera(met::Kerr{T}, θo, αmin, αmax, βmin, βmax, res::Int; A=Matrix) where {T}
-        new{T, A}(met, IntensityScreen(met, αmin, αmax, βmin, βmax, θo, res; A=A), (T(Inf), θo))
+    function IntensityCamera(
+        met::Kerr{T},
+        θo,
+        αmin,
+        αmax,
+        βmin,
+        βmax,
+        res::Int;
+        A = Matrix,
+    ) where {T}
+        new{T,A}(
+            met,
+            IntensityScreen(met, αmin, αmax, βmin, βmax, θo, res; A = A),
+            (T(Inf), θo),
+        )
     end
 end
 
-function η(pix::IntensityPixel) return pix.η end
-function λ(pix::IntensityPixel) return pix.λ end
-function roots(pix::IntensityPixel) return pix.roots end
-function screen_coordinate(pix::IntensityPixel) return pix.screen_coordinate end
-function inclination(pix::IntensityPixel) return pix.θo end
-function I0_inf(pix::IntensityPixel) return pix.I0_inf end
-function total_mino_time(pix::IntensityPixel) return pix.total_mino_time end
-function Ir_inf(pix::IntensityPixel) return pix.I0_inf end
-function absGθo_Gθhat(pix::IntensityPixel) return pix.absGθo_Gθhat end
\ No newline at end of file
+function η(pix::IntensityPixel)
+    return pix.η
+end
+function λ(pix::IntensityPixel)
+    return pix.λ
+end
+function roots(pix::IntensityPixel)
+    return pix.roots
+end
+function screen_coordinate(pix::IntensityPixel)
+    return pix.screen_coordinate
+end
+function inclination(pix::IntensityPixel)
+    return pix.θo
+end
+function I0_inf(pix::IntensityPixel)
+    return pix.I0_inf
+end
+function total_mino_time(pix::IntensityPixel)
+    return pix.total_mino_time
+end
+function Ir_inf(pix::IntensityPixel)
+    return pix.I0_inf
+end
+function absGθo_Gθhat(pix::IntensityPixel)
+    return pix.absGθo_Gθhat
+end
diff --git a/src/cameras/SlowLightIntensityCamera.jl b/src/cameras/SlowLightIntensityCamera.jl
index 1387edd..70e615b 100644
--- a/src/cameras/SlowLightIntensityCamera.jl
+++ b/src/cameras/SlowLightIntensityCamera.jl
@@ -8,7 +8,7 @@ Intensity Pixel Type.
 struct SlowLightIntensityPixel{T} <: AbstractPixel{T}
     metric::Kerr{T}
     "Pixel screen_coordinate"
-    screen_coordinate::NTuple{2, T}
+    screen_coordinate::NTuple{2,T}
     "Radial roots"
     roots::NTuple{4,Complex{T}}
     "Radial antiderivative"
@@ -66,17 +66,22 @@ struct SlowLightIntensityPixel{T} <: AbstractPixel{T}
         I0_inf = Krang.Ir_inf(met, roots)
         new{T}(
             met,
-            (α, β), 
-            roots, 
+            (α, β),
+            roots,
             I0_inf,
             total_mino_time(met, roots),
-            Krang.Iϕ_inf(met, roots, tempλ), 
-            Krang.It_inf(met, roots, tempλ), 
-            I1, I2, Ip, Im, 
-            Krang._absGθo_Gθhat(met, θo, tempη, tempλ), 
-            Krang._absGϕo_Gϕhat(met, θo, tempη, tempλ), 
+            Krang.Iϕ_inf(met, roots, tempλ),
+            Krang.It_inf(met, roots, tempλ),
+            I1,
+            I2,
+            Ip,
+            Im,
+            Krang._absGθo_Gθhat(met, θo, tempη, tempλ),
+            Krang._absGϕo_Gϕhat(met, θo, tempη, tempλ),
             Krang._absGto_Gthat(met, θo, tempη, tempλ),
-            θo, tempη, tempλ
+            θo,
+            tempη,
+            tempλ,
         )
     end
 end
@@ -86,20 +91,29 @@ end
 
 Screen made of `SlowLightIntensityPixel`s.
 """
-struct SlowLightIntensityScreen{T, A <:AbstractMatrix} <: AbstractScreen
+struct SlowLightIntensityScreen{T,A<:AbstractMatrix} <: AbstractScreen
     "Minimum and Maximum Bardeen α values"
-    αrange::NTuple{2, T}
+    αrange::NTuple{2,T}
 
     "Minimum and Maximum Bardeen β values"
-    βrange::NTuple{2, T}
+    βrange::NTuple{2,T}
 
     "Data type that stores screen pixel information"
     pixels::A
 
-    @kernel function _generate_screen!(screen, met::Kerr{T}, αmin, αmax, βmin, βmax, θo, res) where T
-        I,J = @index(Global, NTuple)
-        α = αmin + (αmax - αmin) * (T(I)-1) / (res-1)
-        β = βmin + (βmax - βmin) * (T(J)-1) / (res-1)
+    @kernel function _generate_screen!(
+        screen,
+        met::Kerr{T},
+        αmin,
+        αmax,
+        βmin,
+        βmax,
+        θo,
+        res,
+    ) where {T}
+        I, J = @index(Global, NTuple)
+        α = αmin + (αmax - αmin) * (T(I) - 1) / (res - 1)
+        β = βmin + (βmax - βmin) * (T(J) - 1) / (res - 1)
         screen[I, J] = SlowLightIntensityPixel(met, α, β, θo)
     end
     @doc """
@@ -120,14 +134,33 @@ struct SlowLightIntensityScreen{T, A <:AbstractMatrix} <: AbstractScreen
     # Returns
     A `SlowLightIntensityScreen` object.
     """
-    function SlowLightIntensityScreen(met::Kerr{T}, αmin, αmax, βmin, βmax, θo, res; A=Matrix) where {T}
+    function SlowLightIntensityScreen(
+        met::Kerr{T},
+        αmin,
+        αmax,
+        βmin,
+        βmax,
+        θo,
+        res;
+        A = Matrix,
+    ) where {T}
         screen = A(Matrix{SlowLightIntensityPixel{T}}(undef, res, res))
 
         backend = get_backend(screen)
 
-        _generate_screen!(backend)(screen, met, αmin, αmax, βmin, βmax, θo, res, ndrange = (res, res))
-        
-        new{T, typeof(screen)}((αmin, αmax), (βmin, βmax), screen)
+        _generate_screen!(backend)(
+            screen,
+            met,
+            αmin,
+            αmax,
+            βmin,
+            βmax,
+            θo,
+            res,
+            ndrange = (res, res),
+        )
+
+        new{T,typeof(screen)}((αmin, αmax), (βmin, βmax), screen)
     end
 end
 
@@ -137,12 +170,12 @@ end
 Camera that caches slow light raytracing information for an observer sitting at radial infinity.
 The frame of this observer is alligned with the Boyer-Lindquist frame.
 """
-struct SlowLightIntensityCamera{T, A} <: AbstractCamera
+struct SlowLightIntensityCamera{T,A} <: AbstractCamera
     metric::Kerr{T}
     "Data type that stores screen pixel information"
-    screen::SlowLightIntensityScreen{T, A}
+    screen::SlowLightIntensityScreen{T,A}
     "Observer screen_coordinate"
-    screen_coordinate::NTuple{2, T}
+    screen_coordinate::NTuple{2,T}
     @doc """
         SlowLightIntensityCamera(met::Kerr{T}, θo, αmin, αmax, βmin, βmax, res; A=Matrix) where {T}
 
@@ -161,27 +194,75 @@ struct SlowLightIntensityCamera{T, A} <: AbstractCamera
     # Returns
     - A `SlowLightIntensityCamera` object.
     """
-    function SlowLightIntensityCamera(met::Kerr{T}, θo, αmin, αmax, βmin, βmax, res; A=Matrix) where {T}
+    function SlowLightIntensityCamera(
+        met::Kerr{T},
+        θo,
+        αmin,
+        αmax,
+        βmin,
+        βmax,
+        res;
+        A = Matrix,
+    ) where {T}
         screen = SlowLightIntensityScreen(met, αmin, αmax, βmin, βmax, θo, res; A)
         new{T,typeof(screen.pixels)}(met, screen, (T(Inf), θo))
     end
 end
 
-function η(pix::SlowLightIntensityPixel) return pix.η end
-function λ(pix::SlowLightIntensityPixel) return pix.λ end
-function roots(pix::SlowLightIntensityPixel) return pix.roots end
-function screen_coordinate(pix::SlowLightIntensityPixel) return pix.screen_coordinate end
-function inclination(pix::SlowLightIntensityPixel) return pix.θo end
-function I0_inf(pix::SlowLightIntensityPixel) return pix.I0_inf end
-function total_mino_time(pix::SlowLightIntensityPixel) return pix.total_mino_time end
-function Ir_inf(pix::SlowLightIntensityPixel) return pix.I0_inf end
-function I1_inf_m_I0_terms(pix::SlowLightIntensityPixel) return pix.I1_inf_m_I0_terms end
-function I2_inf_m_I0_terms(pix::SlowLightIntensityPixel) return pix.I2_inf_m_I0_terms end
-function Ip_inf_m_I0_terms(pix::SlowLightIntensityPixel) return pix.Ip_inf_m_I0_terms end
-function Im_inf_m_I0_terms(pix::SlowLightIntensityPixel) return pix.Im_inf_m_I0_terms end
-function radial_inf_integrals_m_I0_terms(pix::SlowLightIntensityPixel) return I1_inf_m_I0_terms(pix), I2_inf_m_I0_terms(pix), Ip_inf_m_I0_terms(pix), Im_inf_m_I0_terms(pix) end
-function Iϕ_inf(pix::SlowLightIntensityPixel) return pix.Iϕ_inf end
-function It_inf(pix::SlowLightIntensityPixel) return pix.It_inf end
-function absGθo_Gθhat(pix::SlowLightIntensityPixel) return pix.absGθo_Gθhat end
-function absGϕo_Gϕhat(pix::SlowLightIntensityPixel) return pix.absGϕo_Gϕhat end
-function absGto_Gthat(pix::SlowLightIntensityPixel) return pix.absGto_Gthat end
+function η(pix::SlowLightIntensityPixel)
+    return pix.η
+end
+function λ(pix::SlowLightIntensityPixel)
+    return pix.λ
+end
+function roots(pix::SlowLightIntensityPixel)
+    return pix.roots
+end
+function screen_coordinate(pix::SlowLightIntensityPixel)
+    return pix.screen_coordinate
+end
+function inclination(pix::SlowLightIntensityPixel)
+    return pix.θo
+end
+function I0_inf(pix::SlowLightIntensityPixel)
+    return pix.I0_inf
+end
+function total_mino_time(pix::SlowLightIntensityPixel)
+    return pix.total_mino_time
+end
+function Ir_inf(pix::SlowLightIntensityPixel)
+    return pix.I0_inf
+end
+function I1_inf_m_I0_terms(pix::SlowLightIntensityPixel)
+    return pix.I1_inf_m_I0_terms
+end
+function I2_inf_m_I0_terms(pix::SlowLightIntensityPixel)
+    return pix.I2_inf_m_I0_terms
+end
+function Ip_inf_m_I0_terms(pix::SlowLightIntensityPixel)
+    return pix.Ip_inf_m_I0_terms
+end
+function Im_inf_m_I0_terms(pix::SlowLightIntensityPixel)
+    return pix.Im_inf_m_I0_terms
+end
+function radial_inf_integrals_m_I0_terms(pix::SlowLightIntensityPixel)
+    return I1_inf_m_I0_terms(pix),
+    I2_inf_m_I0_terms(pix),
+    Ip_inf_m_I0_terms(pix),
+    Im_inf_m_I0_terms(pix)
+end
+function Iϕ_inf(pix::SlowLightIntensityPixel)
+    return pix.Iϕ_inf
+end
+function It_inf(pix::SlowLightIntensityPixel)
+    return pix.It_inf
+end
+function absGθo_Gθhat(pix::SlowLightIntensityPixel)
+    return pix.absGθo_Gθhat
+end
+function absGϕo_Gϕhat(pix::SlowLightIntensityPixel)
+    return pix.absGϕo_Gϕhat
+end
+function absGto_Gthat(pix::SlowLightIntensityPixel)
+    return pix.absGto_Gthat
+end
diff --git a/src/cameras/camera_types.jl b/src/cameras/camera_types.jl
index c93edc2..d5a8677 100644
--- a/src/cameras/camera_types.jl
+++ b/src/cameras/camera_types.jl
@@ -30,8 +30,8 @@ Get the screen coordinates of a given pixel.
 # Returns
 - The screen coordinates of the pixel.
 """
-function screen_coordinate(pix::AbstractPixel) 
-    return pix.screen_coordinate 
+function screen_coordinate(pix::AbstractPixel)
+    return pix.screen_coordinate
 end
 
 """
@@ -45,8 +45,8 @@ Get the space time metric.
 # Returns
 - The space time metric
 """
-function metric(pix::AbstractPixel) 
-    return pix.metric 
+function metric(pix::AbstractPixel)
+    return pix.metric
 end
 
 """
@@ -60,8 +60,8 @@ Get the inclination angle (θo) of a given pixel.
 # Returns
 - The inclination angle (θo) of the pixel.
 """
-function inclination(pix::AbstractPixel) 
-    return pix.θo 
+function inclination(pix::AbstractPixel)
+    return pix.θo
 end
 
 """
@@ -306,4 +306,4 @@ Calculate the absolute value of Gto divided by Gthat for a given pixel.
 """
 function absGto_Gthat(pix::AbstractPixel)
     return _absGto_Gthat(metric(pix), inclination(pix), η(pix), λ(pix))
-end
\ No newline at end of file
+end
diff --git a/src/geometries/geometry_types.jl b/src/geometries/geometry_types.jl
index a3a123f..bafbc33 100644
--- a/src/geometries/geometry_types.jl
+++ b/src/geometries/geometry_types.jl
@@ -24,8 +24,8 @@ struct Mesh{G<:AbstractGeometry,M<:AbstractMaterial}
     material::M
 end
 
-const Scene = NTuple{N, Mesh} where {N}
-Scene() = NTuple{0, Mesh}()
+const Scene = NTuple{N,Mesh} where {N}
+Scene() = NTuple{0,Mesh}()
 
 function add(scene::Scene, mesh::Mesh)
     return (scene..., mesh)
@@ -36,21 +36,27 @@ end
 
 Cone Geometry with half opening angle `opening_angle`.
 """
-struct ConeGeometry{T, A} <: AbstractGeometry
+struct ConeGeometry{T,A} <: AbstractGeometry
     opening_angle::T
     attributes::A
-    ConeGeometry(opening_angle::T) where {T} = new{T, Nothing}(opening_angle)
-    ConeGeometry(opening_angle::T, attributes::A) where {T,A} = new{T, A}(opening_angle,attributes)
+    ConeGeometry(opening_angle::T) where {T} = new{T,Nothing}(opening_angle)
+    ConeGeometry(opening_angle::T, attributes::A) where {T,A} =
+        new{T,A}(opening_angle, attributes)
 end
 
-function _raytrace(observation, pix::AbstractPixel, mesh::Mesh{<:ConeGeometry{T,A}, <:AbstractMaterial}; res) where {T,A}
+function _raytrace(
+    observation,
+    pix::AbstractPixel,
+    mesh::Mesh{<:ConeGeometry{T,A},<:AbstractMaterial};
+    res,
+) where {T,A}
     geometry = mesh.geometry
     material = mesh.material
     θs = geometry.opening_angle
-    subimgs= material.subimgs
+    subimgs = material.subimgs
 
     isindir = false
-    for _ in 1:2 # Looping over isindir this way is needed to get Metal to work
+    for _ = 1:2 # Looping over isindir this way is needed to get Metal to work
         isindir ⊻= true
         for n in subimgs
             #νθ = cos(θs) < abs(cos(θo)) ? (θo > θs) ⊻ (n % 2 == 1) : !isindir
@@ -59,7 +65,8 @@ function _raytrace(observation, pix::AbstractPixel, mesh::Mesh{<:ConeGeometry{T,
                     rs, νr, νθ, _, issuccess = emission_radius(pix, θs, isindir, n)
                     Intersection(zero(rs), rs, θs, zero(rs), νr, νθ), issuccess
                 else
-                    rs, ϕs, νr, νθ, issuccess = emission_coordinates_fast_light(pix, θs, isindir, n)
+                    rs, ϕs, νr, νθ, issuccess =
+                        emission_coordinates_fast_light(pix, θs, isindir, n)
                     Intersection(zero(rs), rs, θs, ϕs, νr, νθ), issuccess
                 end
             else
@@ -80,5 +87,4 @@ function _raytrace(observation, pix::AbstractPixel, mesh::Mesh{<:ConeGeometry{T,
 end
 
 
-Disk(;attributes=nothing) = ConeGeometry(π/2;attributes=attributes)
-
+Disk(; attributes = nothing) = ConeGeometry(π / 2; attributes = attributes)
diff --git a/src/geometries/level_set_geometry.jl b/src/geometries/level_set_geometry.jl
index 401ed0b..da384af 100644
--- a/src/geometries/level_set_geometry.jl
+++ b/src/geometries/level_set_geometry.jl
@@ -3,26 +3,37 @@
 """
 abstract type AbstractLevelSetGeometry{T} <: AbstractGeometry end
 
-function _raytrace(observation::A, pixel::AbstractPixel{T}, mesh::Mesh{<:AbstractLevelSetGeometry, <:AbstractMaterial}; res=100) where {A,T}
+function _raytrace(
+    observation::A,
+    pixel::AbstractPixel{T},
+    mesh::Mesh{<:AbstractLevelSetGeometry,<:AbstractMaterial};
+    res = 100,
+) where {A,T}
     geometry = mesh.geometry
     material = mesh.material
     ray = Vector{Intersection{T}}(undef, res)
     generate_ray!(ray, pixel, res)
-    (;rs, θs, ϕs, νr, νθ) =  ray[1]
-    origin = boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(pixel.metric, rs, θs, ϕs) 
+    (; rs, θs, ϕs, νr, νθ) = ray[1]
+    origin =
+        boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(pixel.metric, rs, θs, ϕs)
     z = zero(A)
-    for i in 2:res
-        (;ts, rs, θs, ϕs, νr, νθ) =  ray[i]
+    for i = 2:res
+        (; ts, rs, θs, ϕs, νr, νθ) = ray[i]
         if rs <= Krang.horizon(pixel.metric)
             continue
         end
-        line_point_2 = boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(pixel.metric, rs, θs, ϕs)
+        line_point_2 = boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(
+            pixel.metric,
+            rs,
+            θs,
+            ϕs,
+        )
         if isinf(line_point_2[1]) || isinf(line_point_2[2]) || isinf(line_point_2[3])
             continue
         end
         didintersect, point = line_intersection(origin, line_point_2, geometry)
         rs = sqrt(sum(point .^ 2))
-        θs = acos(point[3]/rs)
+        θs = acos(point[3] / rs)
         ϕs = atan(point[2], point[1])
 
         if rs < Krang.horizon(pixel.metric)
@@ -33,18 +44,26 @@ function _raytrace(observation::A, pixel::AbstractPixel{T}, mesh::Mesh{<:Abstrac
         observation += !didintersect ? z : material(pixel, intersection)#didintersect ? 1 : 0
         origin = line_point_2
     end
-    return observation 
+    return observation
 end
 
-function line_intersection(origin::NTuple{3,T}, line_point_2, geometry::AbstractLevelSetGeometry{T}) where T
-    didintersect = geometry(origin...)*geometry(line_point_2...) < zero(T)
+function line_intersection(
+    origin::NTuple{3,T},
+    line_point_2,
+    geometry::AbstractLevelSetGeometry{T},
+) where {T}
+    didintersect = geometry(origin...) * geometry(line_point_2...) < zero(T)
     if didintersect
-        direction = (line_point_2[1] - origin[1], line_point_2[2] - origin[2], line_point_2[3] - origin[3])
+        direction = (
+            line_point_2[1] - origin[1],
+            line_point_2[2] - origin[2],
+            line_point_2[3] - origin[3],
+        )
 
-        t = find_zero((x)->geometry((origin .+ (direction .* x))...),  (zero(T), one(T)))
+        t = find_zero((x) -> geometry((origin .+ (direction .* x))...), (zero(T), one(T)))
 
         return true, origin .+ (direction .* t)
     end
     return false, (zero(T), zero(T), zero(T))
 
-end
\ No newline at end of file
+end
diff --git a/src/geometries/mesh_geometry.jl b/src/geometries/mesh_geometry.jl
index 4e4635c..819dd02 100644
--- a/src/geometries/mesh_geometry.jl
+++ b/src/geometries/mesh_geometry.jl
@@ -30,47 +30,53 @@ function rotate(mesh::MeshGeometry, angle, x, y, z)
     return GeometryBasics.Mesh([GeometryBasics.Point(x) for x in points], faces)
 end
 
-function raytrace(camera::AbstractCamera, mesh_geometry::MeshGeometry; res=100) 
+function raytrace(camera::AbstractCamera, mesh_geometry::MeshGeometry; res = 100)
     faces = begin
         temp = getfield(mesh_geometry, :faces)
         len = length(temp)
-        reshape([i[j] for i in temp for j in 1:3], 3, len)
+        reshape([i[j] for i in temp for j = 1:3], 3, len)
     end
     vertices = begin
         temp = mesh_geometry.position
         len = length(temp)
-       reshape([(i[j]) for i in temp for j in 1:3], 3, len)
+        reshape([(i[j]) for i in temp for j = 1:3], 3, len)
     end
     intersections = zeros(Int, size(camera.screen.pixels))
     Threads.@threads for I in CartesianIndices(camera.screen.pixels)
-        intersections[I] = raytrace(camera.screen.pixels[I], faces, vertices;res)
+        intersections[I] = raytrace(camera.screen.pixels[I], faces, vertices; res)
     end
     return intersections
 end
 
-function raytrace(pixel::AbstractPixel{T}, faces::Matrix{GeometryBasics.OffsetInteger{-1, UInt32}}, vertices::Matrix{T}; res=100) where T
+function raytrace(
+    pixel::AbstractPixel{T},
+    faces::Matrix{GeometryBasics.OffsetInteger{-1,UInt32}},
+    vertices::Matrix{T};
+    res = 100,
+) where {T}
     metric = pixel.metric
     intersections = 0
-    ray = Vector{Intersection{T}}(undef,res)
+    ray = Vector{Intersection{T}}(undef, res)
     generate_ray!(ray, pixel, res)
-    (;rs, θs, ϕs) = ray[1]
+    (; rs, θs, ϕs) = ray[1]
     origin = boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(metric, rs, θs, ϕs)#(rs * sin(θs)*cos(ϕs), rs * sin(θs)*sin(ϕs), rs * cos(θs))
 
-    for i in 2:res
-        (;rs, θs, ϕs) = ray[i]
+    for i = 2:res
+        (; rs, θs, ϕs) = ray[i]
         if rs <= Krang.horizon(pixel.metric)
             continue
         end
-        line_point_2 = boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(metric, rs, θs, ϕs)
-        for j in 1:(size(faces)[2])
-            f1, f2, f3 = @view faces[:,j]
-            v1 = @view vertices[:,f1]; 
-            v2 = @view vertices[:,f2]; 
-            v3 = @view vertices[:,f3];
+        line_point_2 =
+            boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(metric, rs, θs, ϕs)
+        for j = 1:(size(faces)[2])
+            f1, f2, f3 = @view faces[:, j]
+            v1 = @view vertices[:, f1]
+            v2 = @view vertices[:, f2]
+            v3 = @view vertices[:, f3]
             didintersect, point = line_intersection(origin, line_point_2, v1, v2, v3)
             intersections += didintersect ? 1 : 0
         end
         origin = line_point_2
     end
     return intersections
-end
\ No newline at end of file
+end
diff --git a/src/materials/AbstractMaterialTypes.jl b/src/materials/AbstractMaterialTypes.jl
index ccb0c01..9e0ea80 100644
--- a/src/materials/AbstractMaterialTypes.jl
+++ b/src/materials/AbstractMaterialTypes.jl
@@ -17,8 +17,8 @@ yield(material::AbstractMaterial) = 0
 abstract type AbstractReturnTrait end
 abstract type AbstractPolarizationTrait <: AbstractReturnTrait end
 struct SimplePolarizationTrait <: AbstractPolarizationTrait end
-struct SimpleIntensityTrait <: AbstractReturnTrait end 
+struct SimpleIntensityTrait <: AbstractReturnTrait end
 
 function returnTrait(mat::AbstractMaterial)
-    return SimpleIntensityTrait() 
+    return SimpleIntensityTrait()
 end
diff --git a/src/materials/ElectronSynchrotronPowerLawIntensity.jl b/src/materials/ElectronSynchrotronPowerLawIntensity.jl
index 027784f..1a24754 100644
--- a/src/materials/ElectronSynchrotronPowerLawIntensity.jl
+++ b/src/materials/ElectronSynchrotronPowerLawIntensity.jl
@@ -1,7 +1,18 @@
 """
 Calculates the intensity of a photon emitted from a fluid particle with momentum f_u and observed by an asymptotic observer.
 """
-function synchrotronIntensity(metric::Kerr{T}, α, β, ri, θs, θo, magnetic_field::SVector{3,T}, βfluid::SVector{3,T}, νr::Bool, θsign::Bool) where {T}
+function synchrotronIntensity(
+    metric::Kerr{T},
+    α,
+    β,
+    ri,
+    θs,
+    θo,
+    magnetic_field::SVector{3,T},
+    βfluid::SVector{3,T},
+    νr::Bool,
+    θsign::Bool,
+) where {T}
 
     βv = βfluid[1]
     θz = βfluid[2]
@@ -17,7 +28,11 @@ function synchrotronIntensity(metric::Kerr{T}, α, β, ri, θs, θo, magnetic_fi
     p_fluid_u = jac_fluid_u_zamo_d(metric, βv, θz, ϕz) * p_zamo_u
     magnetic_fieldx, magnetic_fieldy, magnetic_fieldz = magnetic_field
     _, p_fluid_ux, p_fluid_uy, p_fluid_uz = p_fluid_u ./ p_fluid_u[1]
-    vec = @SVector[p_fluid_uy * magnetic_fieldz - p_fluid_uz * magnetic_fieldy, p_fluid_uz * magnetic_fieldx - p_fluid_ux * magnetic_fieldz, p_fluid_ux * magnetic_fieldy - p_fluid_uy * magnetic_fieldx]
+    vec = @SVector[
+        p_fluid_uy * magnetic_fieldz - p_fluid_uz * magnetic_fieldy,
+        p_fluid_uz * magnetic_fieldx - p_fluid_ux * magnetic_fieldz,
+        p_fluid_ux * magnetic_fieldy - p_fluid_uy * magnetic_fieldx,
+    ]
     norm = √sum(vec .* vec) + eps(T)
 
     return norm, inv(p_fluid_u[1]), abs(p_fluid_u[1] / p_fluid_u[4])
@@ -38,14 +53,14 @@ end
     - `p2::T`: The second power-law index.
     - `subimgs::NTuple{N, Int}`: The sub-images to ray trace.
 """
-struct ElectronSynchrotronPowerLawIntensity{N, T} <: AbstractMaterial 
-    magnetic_field::SVector{3, T}
-    fluid_velocity::SVector{3, T}
+struct ElectronSynchrotronPowerLawIntensity{N,T} <: AbstractMaterial
+    magnetic_field::SVector{3,T}
+    fluid_velocity::SVector{3,T}
     spectral_index::T
     R::T
     p1::T
     p2::T
-    subimgs::NTuple{N, Int}
+    subimgs::NTuple{N,Int}
 
     @doc """
         ElectronSynchrotronPowerLawIntensity(
@@ -79,45 +94,57 @@ struct ElectronSynchrotronPowerLawIntensity{N, T} <: AbstractMaterial
 
     # Returns
     - `ElectronSynchrotronPowerLawIntensity{N, T}`: A new instance of `ElectronSynchrotronPowerLawPolarization`.
-"""   
+"""
     function ElectronSynchrotronPowerLawIntensity(
-        magnetic_fieldx::T, 
-        magnetic_fieldy::T, 
-        magnetic_fieldz::T, 
-        fluid_speed::T, 
-        fluid_inclination_angle::T, 
-        fluid_azimuthal_angle::T, 
+        magnetic_fieldx::T,
+        magnetic_fieldy::T,
+        magnetic_fieldz::T,
+        fluid_speed::T,
+        fluid_inclination_angle::T,
+        fluid_azimuthal_angle::T,
         spectral_index::T,
         R::T,
         p1::T,
         p2::T,
-        subimgs::NTuple{N, Int},
+        subimgs::NTuple{N,Int},
     ) where {N,T}
-        new{N, T}(SVector(magnetic_fieldx, magnetic_fieldy, magnetic_fieldz), SVector(fluid_speed, fluid_inclination_angle, fluid_azimuthal_angle), spectral_index, R, p1, p2, subimgs)
+        new{N,T}(
+            SVector(magnetic_fieldx, magnetic_fieldy, magnetic_fieldz),
+            SVector(fluid_speed, fluid_inclination_angle, fluid_azimuthal_angle),
+            spectral_index,
+            R,
+            p1,
+            p2,
+            subimgs,
+        )
     end
 end
 
 """
     Functor for the the ElectronSynchrotronPowerLawIntensity material
 """
-function (linpol::ElectronSynchrotronPowerLawIntensity{N,T})(pix::AbstractPixel, intersection) where {N,T}
-    (;magnetic_field, fluid_velocity, R, p1, p2, spectral_index) = linpol
-    (;rs, θs, νr, νθ) = intersection
+function (linpol::ElectronSynchrotronPowerLawIntensity{N,T})(
+    pix::AbstractPixel,
+    intersection,
+) where {N,T}
+    (; magnetic_field, fluid_velocity, R, p1, p2, spectral_index) = linpol
+    (; rs, θs, νr, νθ) = intersection
 
     θo = inclination(pix)
     met = metric(pix)
     α, β = screen_coordinate(pix)
 
 
-    norm, redshift, lp = synchrotronIntensity(met, α, β, rs, θs, θo, magnetic_field, fluid_velocity, νr, νθ)
+    norm, redshift, lp =
+        synchrotronIntensity(met, α, β, rs, θs, θo, magnetic_field, fluid_velocity, νr, νθ)
 
 
-    rat = (rs/R)
-    prof = rat^p1/(one(T)+rat^(p1+p2)) * redshift^(T(3) + spectral_index)
+    rat = (rs / R)
+    prof = rat^p1 / (one(T) + rat^(p1 + p2)) * redshift^(T(3) + spectral_index)
 
     # Add a clamp to lp to help remove hot pixels
-    return norm^(one(T) + spectral_index) * min(lp,T(1e2)) * prof
+    return norm^(one(T) + spectral_index) * min(lp, T(1e2)) * prof
 end
 
 isFastLight(material::ElectronSynchrotronPowerLawIntensity) = true
-isAxisymmetric(material::ElectronSynchrotronPowerLawIntensity) = true
\ No newline at end of file
+isAxisymmetric(material::ElectronSynchrotronPowerLawIntensity) = true
diff --git a/src/materials/ElectronSynchrotronPowerLawPolarization.jl b/src/materials/ElectronSynchrotronPowerLawPolarization.jl
index 39c28ac..f904857 100644
--- a/src/materials/ElectronSynchrotronPowerLawPolarization.jl
+++ b/src/materials/ElectronSynchrotronPowerLawPolarization.jl
@@ -19,7 +19,18 @@ evpa(fα, fβ) = atan(-fα, fβ)
 """
 Calculates the polarization of a photon emitted from a fluid particle with momentum f_u and observed by an asymptotic observer.
 """
-function synchrotronPolarization(metric::Kerr{T}, α, β, ri, θs, θo, magnetic_field::SVector{3,T}, βfluid::SVector{3,T}, νr::Bool, θsign::Bool) where {T}
+function synchrotronPolarization(
+    metric::Kerr{T},
+    α,
+    β,
+    ri,
+    θs,
+    θo,
+    magnetic_field::SVector{3,T},
+    βfluid::SVector{3,T},
+    νr::Bool,
+    θsign::Bool,
+) where {T}
 
     a = metric.spin
     βv = βfluid[1]
@@ -36,7 +47,11 @@ function synchrotronPolarization(metric::Kerr{T}, α, β, ri, θs, θo, magnetic
     p_fluid_u = jac_fluid_u_zamo_d(metric, βv, θz, ϕz) * p_zamo_u
     magnetic_fieldx, magnetic_fieldy, magnetic_fieldz = magnetic_field
     _, p_fluid_ux, p_fluid_uy, p_fluid_uz = p_fluid_u ./ p_fluid_u[1]
-    vec = @SVector[p_fluid_uy * magnetic_fieldz - p_fluid_uz * magnetic_fieldy, p_fluid_uz * magnetic_fieldx - p_fluid_ux * magnetic_fieldz, p_fluid_ux * magnetic_fieldy - p_fluid_uy * magnetic_fieldx]
+    vec = @SVector[
+        p_fluid_uy * magnetic_fieldz - p_fluid_uz * magnetic_fieldy,
+        p_fluid_uz * magnetic_fieldx - p_fluid_ux * magnetic_fieldz,
+        p_fluid_ux * magnetic_fieldy - p_fluid_uy * magnetic_fieldx,
+    ]
     norm = √sum(vec .* vec) + eps(T)
     f_fluid_u = SVector(zero(eltype(vec)), vec[1], vec[2], vec[3])
     f_zamo_u = jac_fluid_u_zamo_d(metric, -βv, θz, ϕz) * f_fluid_u
@@ -79,15 +94,15 @@ end
     - `p2::T`: The second power-law index.
     - `subimgs::NTuple{N, Int}`: The sub-images to ray trace.
 """
-struct ElectronSynchrotronPowerLawPolarization{N, T} <: AbstractMaterial 
-    magnetic_field::SVector{3, T}
-    fluid_velocity::SVector{3, T}
+struct ElectronSynchrotronPowerLawPolarization{N,T} <: AbstractMaterial
+    magnetic_field::SVector{3,T}
+    fluid_velocity::SVector{3,T}
     spectral_index::T
     R::T
     p1::T
     p2::T
-    subimgs::NTuple{N, Int}
-    
+    subimgs::NTuple{N,Int}
+
     @doc """
         ElectronSynchrotronPowerLawPolarization(
         magnetic_fieldx::T, 
@@ -120,21 +135,29 @@ struct ElectronSynchrotronPowerLawPolarization{N, T} <: AbstractMaterial
 
     # Returns
     - `ElectronSynchrotronPowerLawPolarization{N, T}`: A new instance of `ElectronSynchrotronPowerLawPolarization`.
-"""   
+"""
     function ElectronSynchrotronPowerLawPolarization(
-        magnetic_fieldx::T, 
-        magnetic_fieldy::T, 
-        magnetic_fieldz::T, 
-        fluid_speed::T, 
-        fluid_inclination_angle::T, 
-        fluid_azimuthal_angle::T, 
+        magnetic_fieldx::T,
+        magnetic_fieldy::T,
+        magnetic_fieldz::T,
+        fluid_speed::T,
+        fluid_inclination_angle::T,
+        fluid_azimuthal_angle::T,
         spectral_index::T,
         R::T,
         p1::T,
         p2::T,
-        subimgs::NTuple{N, Int},
-    ) where {N, T}
-    new{N, T}(SVector(magnetic_fieldx, magnetic_fieldy, magnetic_fieldz), SVector(fluid_speed, fluid_inclination_angle, fluid_azimuthal_angle), spectral_index, R, p1, p2, subimgs)
+        subimgs::NTuple{N,Int},
+    ) where {N,T}
+        new{N,T}(
+            SVector(magnetic_fieldx, magnetic_fieldy, magnetic_fieldz),
+            SVector(fluid_speed, fluid_inclination_angle, fluid_azimuthal_angle),
+            spectral_index,
+            R,
+            p1,
+            p2,
+            subimgs,
+        )
     end
 end
 
@@ -149,28 +172,43 @@ end
 """
     Functor for the ElectronSynchrotronPowerLawPolarization material
 """
-function (linpol::ElectronSynchrotronPowerLawPolarization{N,T})(pix::AbstractPixel, intersection) where {N,T}
-    (;magnetic_field, fluid_velocity, spectral_index, R, p1, p2) = linpol
-    (;rs, θs, νr, νθ) = intersection
+function (linpol::ElectronSynchrotronPowerLawPolarization{N,T})(
+    pix::AbstractPixel,
+    intersection,
+) where {N,T}
+    (; magnetic_field, fluid_velocity, spectral_index, R, p1, p2) = linpol
+    (; rs, θs, νr, νθ) = intersection
+
 
-    
     θo = inclination(pix)
     met = metric(pix)
     α, β = screen_coordinate(pix)
 
-    eα, eβ, redshift, lp = synchrotronPolarization(met, α, β, rs, θs, θo, magnetic_field, fluid_velocity, νr, νθ)
-
-    rat = (rs/R)
-    prof = rat^p1/(1+rat^(p1+p2)) * redshift^(T(3) + spectral_index)
+    eα, eβ, redshift, lp = synchrotronPolarization(
+        met,
+        α,
+        β,
+        rs,
+        θs,
+        θo,
+        magnetic_field,
+        fluid_velocity,
+        νr,
+        νθ,
+    )
+
+    rat = (rs / R)
+    prof = rat^p1 / (1 + rat^(p1 + p2)) * redshift^(T(3) + spectral_index)
     q = T(-(eα^2 - eβ^2) + eps(T))
     u = T(-2 * eα * eβ + eps(T))
 
     # Add a clamp to lp to help remove hot pixels
-    i = hypot(q, u)^(one(T) + spectral_index) * min(lp,T(1e2)) * prof
+    i = hypot(q, u)^(one(T) + spectral_index) * min(lp, T(1e2)) * prof
     return StokesParams(i, q, u, zero(T))
 end
 
 isOcclusive(::ElectronSynchrotronPowerLawPolarization) = false
 isFastLight(::ElectronSynchrotronPowerLawPolarization) = true
 isAxisymmetric(::ElectronSynchrotronPowerLawPolarization) = true
-yield(::ElectronSynchrotronPowerLawPolarization{N,T}) where {N,T} = StokesParams(zero(T), zero(T), zero(T), zero(T))
+yield(::ElectronSynchrotronPowerLawPolarization{N,T}) where {N,T} =
+    StokesParams(zero(T), zero(T), zero(T), zero(T))
diff --git a/src/materials/physicsUtils.jl b/src/materials/physicsUtils.jl
index adbd5ae..59fac68 100644
--- a/src/materials/physicsUtils.jl
+++ b/src/materials/physicsUtils.jl
@@ -1,4 +1,13 @@
-export p_bl_d, penrose_walker, screen_polarization, PowerLaw, evpa, jac_zamo_d_bl_u, jac_bl_d_zamo_u, jac_zamo_u_bl_d, jac_bl_u_zamo_d, jac_fluid_u_zamo_d
+export p_bl_d,
+    penrose_walker,
+    screen_polarization,
+    PowerLaw,
+    evpa,
+    jac_zamo_d_bl_u,
+    jac_bl_d_zamo_u,
+    jac_zamo_u_bl_d,
+    jac_bl_u_zamo_d,
+    jac_fluid_u_zamo_d
 ##----------------------------------------------------------------------------------------------------------------------
 #Polarization stuff
 ##----------------------------------------------------------------------------------------------------------------------
@@ -8,9 +17,10 @@ export p_bl_d, penrose_walker, screen_polarization, PowerLaw, evpa, jac_zamo_d_b
 function p_bl_d(metric::Kerr{T}, r, θ, η, λ, νr::Bool, νθ::Bool) where {T}
     return @SVector[
         -one(T),
-        (νr ? one(T) : -one(T))* √max(zero(T), r_potential(metric, η, λ, r)) / Δ(metric, r),
-        (νθ ? one(T) : -one(T))* √max(zero(T), θ_potential(metric, η, λ, θ)),
-        λ
+        (νr ? one(T) : -one(T)) * √max(zero(T), r_potential(metric, η, λ, r)) /
+        Δ(metric, r),
+        (νθ ? one(T) : -one(T)) * √max(zero(T), θ_potential(metric, η, λ, θ)),
+        λ,
     ]
 end
 
@@ -107,7 +117,13 @@ end
 """
 Returns the Penrose walker constant for a photon with momentum p_u emitted from a fluid particle with momentum f_u.
 """
-function penrose_walker(metric::Kerr{T}, r, θ, p_u::AbstractVector, f_u::AbstractVector) where {T}# Eq 6 arXiv:2001.08750v1
+function penrose_walker(
+    metric::Kerr{T},
+    r,
+    θ,
+    p_u::AbstractVector,
+    f_u::AbstractVector,
+) where {T}# Eq 6 arXiv:2001.08750v1
     a = metric.spin
     pt, pr, pϕ, pθ = p_u
     ft, fr, fϕ, fθ = f_u
@@ -118,4 +134,3 @@ function penrose_walker(metric::Kerr{T}, r, θ, p_u::AbstractVector, f_u::Abstra
     B = ((r * r + a * a) * (pϕ * fθ - pθ * fϕ) - a * (pt * fθ - pθ * ft)) * sinθ
     return A * r - B * a * cosθ, -(A * a * cosθ + B * r)
 end
-
diff --git a/src/metrics/AbstractMetric.jl b/src/metrics/AbstractMetric.jl
index db3cd13..4e17211 100644
--- a/src/metrics/AbstractMetric.jl
+++ b/src/metrics/AbstractMetric.jl
@@ -16,6 +16,9 @@ end
 Returns the inverse metric in some representation (usually as an nxn matrix).
 """
 function metric_uu(metric::AbstractMetric, args...)
-    try return inv(metric_dd(metric, args...)) catch err throw(err) end 
+    try
+        return inv(metric_dd(metric, args...))
+    catch err
+        throw(err)
+    end
 end
-
diff --git a/src/metrics/Kerr/Kerr.jl b/src/metrics/Kerr/Kerr.jl
index 92f0327..9b45cfd 100644
--- a/src/metrics/Kerr/Kerr.jl
+++ b/src/metrics/Kerr/Kerr.jl
@@ -9,18 +9,18 @@ $FIELDS
 """
 struct Kerr{T} <: AbstractMetric
     "M = mass"
-    mass::T  
+    mass::T
     "a = J/M, where J is the angular momentum and M is the mass of the black hole."
     spin::T
 
     @doc """
         Kerr(spin::T) where {T}
-    
+
     Constructs a `Kerr` object representing the Kerr metric.
-    
+
     # Arguments
     - `spin::T`: The spin parameter `a = J/M`, where `J` is the angular momentum and `M` is the mass of the black hole. spin ∈ (-1, 0) ∪ (0, 1).
-    
+
     # Returns
     - A `Kerr` object with the given spin and a default mass of 1.
     """
@@ -43,15 +43,15 @@ end
 function Σ(metric::Kerr{T}, r, θ) where {T}
     return r^2 + metric.spin^2 * cos(θ)^2
 end
-function A(metric::Kerr{T}, r, θ; Δ=Δ(metric, r)) where {T}
+function A(metric::Kerr{T}, r, θ; Δ = Δ(metric, r)) where {T}
     a = metric.spin
     return (r^2 + a^2)^2 - a^2 * Δ * sin(θ)^2
 end
-function Ξ(metric::Kerr{T}, r, θ; Δ=Δ(metric, r)) where {T}
+function Ξ(metric::Kerr{T}, r, θ; Δ = Δ(metric, r)) where {T}
     a = metric.spin
     return (r^2 + a^2)^2 - Δ * a^2 * sin(θ)^2
 end
-function ω(metric::Kerr{T}, r, θ; Ξ=Ξ(metric, r, θ)) where {T}
+function ω(metric::Kerr{T}, r, θ; Ξ = Ξ(metric, r, θ)) where {T}
     return T(2) * metric.spin * r / Ξ
 end
 
@@ -61,9 +61,9 @@ Inverse Kerr metric in Boyer Lindquist (BL) coordinates.
 """
 function metric_uu(metric::Kerr{T}, r, θ) where {T}
     Δt = Δ(metric, r)
-    Ξt = Ξ(metric, r, θ; Δ=Δt)
+    Ξt = Ξ(metric, r, θ; Δ = Δt)
     Σt = Σ(metric, r, θ)
-    ωt = ω(metric, r, θ; Ξ=Ξt)
+    ωt = ω(metric, r, θ; Ξ = Ξt)
     z = zero(T)
 
     return @SMatrix [ #Eq 1 2105.09440
@@ -92,9 +92,9 @@ Kerr metric in Boyer Lindquist (BL) coordinates.
 """
 function metric_dd(metric::Kerr{T}, r, θ) where {T}
     Δt = Δ(metric, r)
-    Ξt = Ξ(metric, r, θ; Δ=Δt)
+    Ξt = Ξ(metric, r, θ; Δ = Δt)
     Σt = Σ(metric, r, θ)
-    ωt = ω(metric, r, θ; Ξ=Ξt)
+    ωt = ω(metric, r, θ; Ξ = Ξt)
     sin2 = sin(θ)^2
     z = zero(T)
 
@@ -119,4 +119,3 @@ function metric_dd(metric::Kerr{T}, coordinates) where {T}
     _, r, θ, _ = coordinates
     return metric_dd(metric, r, θ)
 end
-
diff --git a/src/metrics/Kerr/emission_coordinates.jl b/src/metrics/Kerr/emission_coordinates.jl
index 31a6992..ad00ded 100644
--- a/src/metrics/Kerr/emission_coordinates.jl
+++ b/src/metrics/Kerr/emission_coordinates.jl
@@ -2,11 +2,19 @@
 # Functions generically return 0 when the emission coordinates do not exist for a given screen coordinate to play nice 
 # with Enzyme's AD.
 
-export emission_radius, emission_inclination, emission_coordinates_fast_light, emission_coordinates, emission_coordinates!
-
-function emission_coordinates!(coordinates::Array{T,3}, camera::AbstractCamera, τ::T) where {T} 
+export emission_radius,
+    emission_inclination,
+    emission_coordinates_fast_light,
+    emission_coordinates,
+    emission_coordinates!
+
+function emission_coordinates!(
+    coordinates::Array{T,3},
+    camera::AbstractCamera,
+    τ::T,
+) where {T}
     for I in CartesianIndices(camera.screen.pixels)
-        coordinates[:,I] .= emission_coordinates(camera.screen.pixels[I], τ)[1:4]
+        coordinates[:, I] .= emission_coordinates(camera.screen.pixels[I], τ)[1:4]
     end
 end
 
@@ -21,11 +29,16 @@ Returns 0 if the emission coordinates do not exist for that screen coordinate.
 - `isindir` : Is emission to observer direct or indirect
 - `n` : Image index
 """
-function emission_radius(pix::Krang.AbstractPixel, θs::T, isindir, n)::Tuple{T, Bool, Bool, Int, Bool} where {T}
+function emission_radius(
+    pix::Krang.AbstractPixel,
+    θs::T,
+    isindir,
+    n,
+)::Tuple{T,Bool,Bool,Int,Bool} where {T}
     α, β = screen_coordinate(pix)
     θo = inclination(pix)
     met = metric(pix)
-    isincone = θo ≤ θs ≤ (π-θo) || (π-θo) ≤ θs ≤ θo
+    isincone = θo ≤ θs ≤ (π - θo) || (π - θo) ≤ θs ≤ θo
     if !isincone
         αmin = αboundary(met, θs)
         βbound = (abs(α) >= (αmin + eps(T)) ? βboundary(met, α, θo, θs) : zero(T))
@@ -37,8 +50,8 @@ function emission_radius(pix::Krang.AbstractPixel, θs::T, isindir, n)::Tuple{T,
 
     # νθ is θ̇s increasing or decreasing?
     νθ = isincone ? (θo > θs) ⊻ (n % 2 == 1) : !isindir
-    if isincone 
-        νθ = (θo > θs) ⊻ (n % 2 == 1) 
+    if isincone
+        νθ = (θo > θs) ⊻ (n % 2 == 1)
     end
 
     ans, νr, numreals, issuccess = emission_radius(pix, τ)
@@ -54,7 +67,7 @@ Returns 0 if the emission coordinates do not exist for that screen coordinate.
 -`pix` : Pixel information
 - `τ` : Mino time
 """
-function emission_radius(pix::AbstractPixel, τ::T)::Tuple{T, Bool, Int, Bool} where {T}
+function emission_radius(pix::AbstractPixel, τ::T)::Tuple{T,Bool,Int,Bool} where {T}
     met = metric(pix)
     a = met.spin
     ans = zero(T)
@@ -126,7 +139,7 @@ function emission_azimuth(pix::AbstractPixel, θs, rs, τ::T, νr, isindir, n) w
     Gϕtemp, _, _, _ = Gϕ(pix, θs, isindir, n)
     (isnan(Gϕtemp) || !isfinite(Gϕtemp)) && return T(NaN)
 
-    return -(Iϕ + λtemp * Gϕtemp + T(π/2))
+    return -(Iϕ + λtemp * Gϕtemp + T(π / 2))
 end
 
 """
@@ -140,12 +153,12 @@ coordinate (`α`, `β`) for an observer located at inclination θo.
 - `isindir` : Whether emission to observer is direct or indirect
 - `n` : Image index
 """
-function emission_coordinates_fast_light(pix::AbstractPixel, θs::T, isindir, n) where T
+function emission_coordinates_fast_light(pix::AbstractPixel, θs::T, isindir, n) where {T}
     # NB: I do not return θs since it is already known, and returning it might encourage bad coding errors
     α, β = screen_coordinate(pix)
     θo = inclination(pix)
     met = metric(pix)
-    isincone = θo ≤ θs ≤ (π-θo) || (π-θo) ≤ θs ≤ θo
+    isincone = θo ≤ θs ≤ (π - θo) || (π - θo) ≤ θs ≤ θo
     if !isincone
         αmin = αboundary(met, θs)
         βbound = (abs(α) >= (αmin + eps(T)) ? βboundary(met, α, θo, θs) : zero(T))
@@ -157,8 +170,8 @@ function emission_coordinates_fast_light(pix::AbstractPixel, θs::T, isindir, n)
 
     # is θ̇s increasing or decreasing?
     νθ = isincone ? (θo > θs) ⊻ (n % 2 == 1) : !isindir
-    if isincone 
-        νθ = (θo > θs) ⊻ (n % 2 == 1) 
+    if isincone
+        νθ = (θo > θs) ⊻ (n % 2 == 1)
     end
 
     rs, νr, _, issuccess = emission_radius(pix, τ)
@@ -208,7 +221,7 @@ function emission_coordinates_fast_light(pix::AbstractPixel, τ::T) where {T}
 
     Gϕtemp, _, _, _, _ = Gϕ(pix, θs, isindir, n)
 
-    emission_azimuth = -(Iϕ + λtemp * Gϕtemp + T(π/2))
+    emission_azimuth = -(Iϕ + λtemp * Gϕtemp + T(π / 2))
 
     νθ = abs(cos(θs)) < abs(cos(θo)) ? (n % 2 == 1) ⊻ (θo > θs) : !isindir ⊻ (θs > T(π / 2))
     return rs, θs, emission_azimuth, νr, νθ, issuccess
@@ -265,13 +278,17 @@ function emission_coordinates(pix::AbstractPixel, θs::T, isindir, n) where {T}
     rp = one(T) + √(one(T) - a^2)
     rm = one(T) - √(one(T) - a^2)
 
-    It = 4 / (rp - rm) * (rp * (rp - a * λtemp / 2) * Ip - rm * (rm - a * λtemp / 2) * Im) + 4 * I0 + 2 * I1 + I2
+    It =
+        4 / (rp - rm) * (rp * (rp - a * λtemp / 2) * Ip - rm * (rm - a * λtemp / 2) * Im) +
+        4 * I0 +
+        2 * I1 +
+        I2
     Iϕ = 2a / (rp - rm) * ((rp - a * λtemp / 2) * Ip - (rm - a * λtemp / 2) * Im)
 
     Gϕtemp, _, _, _, _ = Gϕ(pix, θs, isindir, n)
     Gttemp, _, _, _, _ = Gt(pix, θs, isindir, n)
 
-    emission_azimuth = -(Iϕ + λtemp * Gϕtemp + T(π/2))
+    emission_azimuth = -(Iϕ + λtemp * Gϕtemp + T(π / 2))
     emission_time_regularized = (zero(T) + It + a^2 * Gttemp)
 
     # is θ̇s increasing or decreasing?
@@ -315,13 +332,17 @@ function emission_coordinates(pix::AbstractPixel, τ::T) where {T}
     rp = one(T) + √(one(T) - a^2)
     rm = one(T) - √(one(T) - a^2)
 
-    It = 4 / (rp - rm) * (rp * (rp - a * λtemp / 2) * Ip - rm * (rm - a * λtemp / 2) * Im) + 4 * I0 + 2 * I1 + I2
+    It =
+        4 / (rp - rm) * (rp * (rp - a * λtemp / 2) * Ip - rm * (rm - a * λtemp / 2) * Im) +
+        4 * I0 +
+        2 * I1 +
+        I2
     Iϕ = 2a / (rp - rm) * ((rp - a * λtemp / 2) * Ip - (rm - a * λtemp / 2) * Im)
 
     Gϕtemp, _, _, _, _ = Gϕ(pix, θs, isindir, n)
     Gttemp, _, _, _, _ = Gt(pix, θs, isindir, n)
 
-    emission_azimuth = -(Iϕ + λtemp * Gϕtemp + T(π/2))
+    emission_azimuth = -(Iϕ + λtemp * Gϕtemp + T(π / 2))
     emission_time_regularized = (It + a^2 * Gttemp)
 
     νθ = abs(cos(θs)) < abs(cos(θo)) ? (n % 2 == 1) ⊻ (θo > θs) : !isindir ⊻ (θs > T(π / 2))
@@ -361,7 +382,7 @@ function _θs(metric::Kerr{T}, signβ, θo, η, λ, τ) where {T}
         argr = (JacobiElliptic.sn(absτs / tempfac, k))^2
         ans = acos((θo > T(π / 2) ? -one(T) : one(T)) * √((up - um) * argr + um))
         sign = ans > T(π / 2) ? -one(T) : one(T)
-        τo = sign  * τo
+        τo = sign * τo
     else
         argo = clamp(cos(θo) / √(up), -one(T), one(T))
         k = m
@@ -371,7 +392,7 @@ function _θs(metric::Kerr{T}, signβ, θo, η, λ, τ) where {T}
         τhat = Ghat_2 + Ghat_2
         Δτtemp = (τ % τhat + signβ * τo)
         n = floor(τ / τhat)
-        sign = (n%2 == 1) ? -one(T) : one(T)
+        sign = (n % 2 == 1) ? -one(T) : one(T)
         τs = argmin(abs, (sign * signβ * (τhat - Δτtemp), sign * signβ * Δτtemp))
         newargs = JacobiElliptic.sn(τs / tempfac, k)
         ans = acos(√up * newargs)
@@ -382,17 +403,14 @@ function _θs(metric::Kerr{T}, signβ, θo, η, λ, τ) where {T}
     if isincone
         νθ = (n % 2 == 1) ⊻ (θo > ans)
         τ1 = (n + 1) * τhat - signβ * τo + (νθ ? 1 : -1) * τs
-        isindir = isapprox(τ1, τ, rtol=sqrt(eps(T)))
+        isindir = isapprox(τ1, τ, rtol = sqrt(eps(T)))
     elseif ans < T(π / 2)
         τ1 = (n + 1) * τhat - signβ * τo - τs
-        isindir = isapprox(τ1, τ, rtol=sqrt(eps(T)))
+        isindir = isapprox(τ1, τ, rtol = sqrt(eps(T)))
     else
         τ1 = (n + 1) * τhat - signβ * τo + τs
-        isindir = isapprox(τ1, τ, rtol=sqrt(eps(T)))
+        isindir = isapprox(τ1, τ, rtol = sqrt(eps(T)))
     end
 
     return ans, τs, τo, τhat, n, isindir
 end
-
-
-
diff --git a/src/metrics/Kerr/misc.jl b/src/metrics/Kerr/misc.jl
index d64969d..3d8e35e 100644
--- a/src/metrics/Kerr/misc.jl
+++ b/src/metrics/Kerr/misc.jl
@@ -1,16 +1,26 @@
 # Miscellaneous functions
 ##----------------------------------------------------------------------------------------------------------------------
-export λ, η, α, β, αboundary, βboundary, r_potential, θ_potential, get_radial_roots, 
-Ir,Gθ, total_mino_time
+export λ,
+    η,
+    α,
+    β,
+    αboundary,
+    βboundary,
+    r_potential,
+    θ_potential,
+    get_radial_roots,
+    Ir,
+    Gθ,
+    total_mino_time
 
 """
 Checks if a complex number is real to some tolerance
 """
-function _isreal2(num::Complex{T}) where T
+function _isreal2(num::Complex{T}) where {T}
     ren, imn = reim(num)
     ren2 = ren^2
     imn2 = imn^2
-    return imn2 / (imn2+ren2) < eps(T)
+    return imn2 / (imn2 + ren2) < eps(T)
 end
 
 """
@@ -35,10 +45,7 @@ end
 function regularized_R2(α::T, φ::T, j) where {T}
 
     return one(T) / (one(T) - α^2) * (
-        JacobiElliptic.F(φ, j)
-        -
-        α^2 / (j + (one(T) - j) * α^2) * (
-            JacobiElliptic.E(φ, j)
+        JacobiElliptic.F(φ, j) - α^2 / (j + (one(T) - j) * α^2) * (JacobiElliptic.E(φ, j)
         #-α*sin(φ)*sqrt(one(T)-j*sin(φ)^2)/(one(T)+α*cos(φ)) #Linear Divergent term
         )
     ) #+ 
@@ -52,8 +59,9 @@ end
 
 function regularizedS2(α, φ, j)
 
-    return -inv((1 + α^2) * (1 - j + α^2)) * ((1 - j) * JacobiElliptic.F(φ, j) + α^2 * JacobiElliptic.E(φ, j))+# + α^2*√(1-j*sin(φ)^2)-α^3) +
-    (inv(1 + α^2) + (1 - j) / (1 - j + α^2)) * regularizedS1(α, φ, j)
+    return -inv((1 + α^2) * (1 - j + α^2)) *
+           ((1 - j) * JacobiElliptic.F(φ, j) + α^2 * JacobiElliptic.E(φ, j)) +# + α^2*√(1-j*sin(φ)^2)-α^3) +
+           (inv(1 + α^2) + (1 - j) / (1 - j + α^2)) * regularizedS1(α, φ, j)
 end
 
 function p1(α::T, j::T) where {T}
@@ -68,32 +76,43 @@ end
 
 function f2(α, sinφ, j)
     p2 = √((1 + α^2) / (1 - j + α^2))
-    return p2 / 2 * log(abs((1 - p2) / (1 + p2) * (1 + p2 * √(1 - j * sinφ^2)) / (1 - p2 * √(1 - j * sinφ^2))))
+    return p2 / 2 * log(
+        abs(
+            (1 - p2) / (1 + p2) * (1 + p2 * √(1 - j * sinφ^2)) /
+            (1 - p2 * √(1 - j * sinφ^2)),
+        ),
+    )
 end
 
 function R1(α::T, φ::T, j) where {T}
-    return one(T) / (one(T) - α^2) * (JacobiElliptic.Pi(α^2 / (α^2 - one(T)), φ, j) - α * f1(α, sin(φ), j))
+    return one(T) / (one(T) - α^2) *
+           (JacobiElliptic.Pi(α^2 / (α^2 - one(T)), φ, j) - α * f1(α, sin(φ), j))
 end
 
 function R2(α::T, φ::T, j) where {T}
     return one(T) / (one(T) - α^2) * (
-        JacobiElliptic.F(φ, j)
-        -
+        JacobiElliptic.F(φ, j) -
         α^2 / (j + (one(T) - j) * α^2) * (
-            JacobiElliptic.E(φ, j) - α * sin(φ) * √(one(T) - j * sin(φ)^2) / (one(T) + α * cos(φ))
+            JacobiElliptic.E(φ, j) -
+            α * sin(φ) * √(one(T) - j * sin(φ)^2) / (one(T) + α * cos(φ))
         )
     ) + inv(j + (one(T) - j) * α^2) * (2 * j - α^2 / (α^2 - one(T))) * R1(α, φ, j)
 end
 
 function S1(α, φ, j)
     α2 = α * α
-    return inv(1 + α2) * (JacobiElliptic.F(φ, j) + α2 * JacobiElliptic.Pi(1 + α2, φ, j) - α * f2(α, sin(φ), j))
+    return inv(1 + α2) * (
+        JacobiElliptic.F(φ, j) + α2 * JacobiElliptic.Pi(1 + α2, φ, j) -
+        α * f2(α, sin(φ), j)
+    )
 end
 
 function S2(α, φ, j)
-    return -inv((1 + α^2) * (1 - j + α^2)) * ((1 - j) * JacobiElliptic.F(φ, j) + 
-        α^2 * JacobiElliptic.E(φ, j) + α^2 * √(1 - j * sin(φ)^2) * (α - tan(φ)) / (1 + α * tan(φ)) - α^3) +
-        (inv(1 + α^2) + (1 - j) / (1 - j + α^2)) * S1(α, φ, j)
+    return -inv((1 + α^2) * (1 - j + α^2)) * (
+        (1 - j) * JacobiElliptic.F(φ, j) +
+        α^2 * JacobiElliptic.E(φ, j) +
+        α^2 * √(1 - j * sin(φ)^2) * (α - tan(φ)) / (1 + α * tan(φ)) - α^3
+    ) + (inv(1 + α^2) + (1 - j) / (1 - j + α^2)) * S1(α, φ, j)
 end
 
 """
@@ -175,7 +194,10 @@ Defines a vertical boundary on the assymptotic observer's screen that emission t
 function βboundary(metric::Kerr{T}, α, θo, θs) where {T}
     a = metric.spin
     cosθs2 = cos(θs)^2
-    return √max((cos(θo)^2 - cosθs2) * (α^2 - a^2 * (one(T) - cosθs2)) / (cosθs2 - one(T)), zero(T)) #eq 15 DOI 10.3847/1538-4357/acafe3 
+    return √max(
+        (cos(θo)^2 - cosθs2) * (α^2 - a^2 * (one(T) - cosθs2)) / (cosθs2 - one(T)),
+        zero(T),
+    ) #eq 15 DOI 10.3847/1538-4357/acafe3 
 end
 
 """
@@ -191,7 +213,8 @@ Radial potential of spacetime
 function r_potential(metric::Kerr{T}, η, λ, r) where {T}
     a = metric.spin
     λ2 = λ^2
-    return a * (a * (r * (r + 2) - η) - 4 * λ * r) + r * ((η + η) + (λ2 + λ2) + r * (-η - λ2 + r^2)) # Eq 7 PhysRevD.101.044032
+    return a * (a * (r * (r + 2) - η) - 4 * λ * r) +
+           r * ((η + η) + (λ2 + λ2) + r * (-η - λ2 + r^2)) # Eq 7 PhysRevD.101.044032
 end
 
 """
@@ -237,7 +260,7 @@ function get_radial_roots(metric::Kerr{T}, η, λ) where {T}
     ωp = ^(-Q / T(2) + sqrt(-Δ3 / T(108)) + zero(T)im, T(1 / 3))
 
     #C = ((-1+0im)^(2/3), (-1+0im)^(4/3), 1) .* ωp
-    C = (-T(1/2) + T(√3/2)im, -T(1/2) - T(√3/2)im, one(T) + zero(T)im) .* ωp
+    C = (-T(1 / 2) + T(√3 / 2)im, -T(1 / 2) - T(√3 / 2)im, one(T) + zero(T)im) .* ωp
 
     v = -P .* inv.(T(3) .* C)
 
@@ -260,10 +283,10 @@ function get_radial_roots(metric::Kerr{T}, η, λ) where {T}
     if (sum(_isreal2.(roots)) == 2) && (abs(imag(roots[4])) < sqrt(eps(T)))
         roots = (roots[1], roots[4], roots[2], roots[3])
     end
-    return roots 
+    return roots
 end
 
-function _get_root_diffs(r1::T, r2::T, r3::T, r4::T) where T 
+function _get_root_diffs(r1::T, r2::T, r3::T, r4::T) where {T}
     return r2 - r1, r3 - r1, r3 - r2, r4 - r1, r4 - r2, r4 - r3
 end
 
@@ -347,7 +370,7 @@ function Ir_inf_case4(::Kerr{T}, roots::NTuple{4}) where {T}
     go = √max((T(4)a2^2 - (C - D)^2) / ((C + D)^2 - T(4)a2^2), zero(T))
     coef = 2 / (C + D)
 
-    Ir_inf = coef*JacobiElliptic.F(T(π / 2) + atan(go), k4)
+    Ir_inf = coef * JacobiElliptic.F(T(π / 2) + atan(go), k4)
     #return (C+D), coef
     return Ir_inf
 end
@@ -385,7 +408,7 @@ function Ir_s_case1_and_2(::Kerr{T}, rs, roots::NTuple{4}, νr) where {T}
     coef = 2 / √real(r31 * r42)
     Ir_s = (x2_s > one(T)) ? T(Inf) : coef * JacobiElliptic.F(asin(x2_s), k)
 
-    return -(-1)^νr*Ir_s
+    return -(-1)^νr * Ir_s
 end
 
 function Ir_s_case3(::Kerr{T}, rs, roots::NTuple{4}) where {T}
@@ -425,7 +448,7 @@ function Ir_s_case4(::Kerr{T}, rs, roots::NTuple{4}) where {T}
     go = √max((T(4)a2^2 - (C - D)^2) / ((C + D)^2 - T(4)a2^2), zero(T))
     x4_s = (rs + b1) / a2
     coef = 2 / (C + D)
-    Ir_s = coef*JacobiElliptic.F(atan(x4_s) + atan(go), k4)
+    Ir_s = coef * JacobiElliptic.F(atan(x4_s) + atan(go), k4)
     return Ir_s
 end
 
@@ -441,7 +464,7 @@ See [`r_potential(x)`](@ref) for an implementation of \$\\mathcal{R}(r)\$.
 - `roots` : Radial roots
 - `λ`  : Reduced azimuthal angular momentum
 """
-function Iϕ_inf(metric::Kerr{T}, roots, λ) where T
+function Iϕ_inf(metric::Kerr{T}, roots, λ) where {T}
     numreals = sum(_isreal2.(roots))
 
     if numreals == 4 #case 2
@@ -449,7 +472,7 @@ function Iϕ_inf(metric::Kerr{T}, roots, λ) where T
     elseif numreals == 2 #case3
         return Iϕ_inf_case3(metric, roots, λ)
     end #case4
-        return Iϕ_inf_case4(metric, roots, λ)
+    return Iϕ_inf_case4(metric, roots, λ)
 end
 
 function Iϕ_inf_case2(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
@@ -473,9 +496,10 @@ function Iϕ_inf_case2(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
     Πm_o = coef_m * JacobiElliptic.Pi(rm3 * r41 / (rm4 * r31), asin(x2_o), k)
 
     Ipo_inf_m_I0_terms = -Πp_o
-    Imo_inf_m_I0_terms = -Πm_o 
+    Imo_inf_m_I0_terms = -Πm_o
 
-    return 2a / (rp - rm) * ((rp - a * λ / 2) * Ipo_inf_m_I0_terms - (rm - a * λ / 2) * Imo_inf_m_I0_terms)
+    return 2a / (rp - rm) *
+           ((rp - a * λ / 2) * Ipo_inf_m_I0_terms - (rm - a * λ / 2) * Imo_inf_m_I0_terms)
 end
 
 function Iϕ_inf_case3(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
@@ -512,8 +536,8 @@ function Iϕ_inf_case3(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
     R1p_o = R1(αp, φ_o, k3)
     R1m_o = R1(αm, φ_o, k3)
 
-    Ip = -inv(B * rp2 + A * rp1) * ( 2 * r21 * √(A * B) / (B * rp2 - A * rp1) * R1p_o)
-    Im = -inv(B * rm2 + A * rm1) * ( 2 * r21 * √(A * B) / (B * rm2 - A * rm1) * R1m_o)
+    Ip = -inv(B * rp2 + A * rp1) * (2 * r21 * √(A * B) / (B * rp2 - A * rp1) * R1p_o)
+    Im = -inv(B * rm2 + A * rm1) * (2 * r21 * √(A * B) / (B * rm2 - A * rm1) * R1m_o)
 
     return (2a / (rp - rm) * ((rp - a * λ / 2) * Ip - (rm - a * λ / 2) * Im))
 end
@@ -549,8 +573,12 @@ function Iϕ_inf_case4(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
     S1p_o = S1(gp, T(π / 2) + atan(go), k4)
     S1m_o = S1(gm, T(π / 2) + atan(go), k4)
 
-    Ip = go / (a2 * (1 - go * x4_p)) * ( - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (S1p_o))
-    Im = go / (a2 * (1 - go * x4_m)) * ( - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (S1m_o))
+    Ip =
+        go / (a2 * (1 - go * x4_p)) *
+        (-2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (S1p_o))
+    Im =
+        go / (a2 * (1 - go * x4_m)) *
+        (-2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (S1m_o))
 
     return 2a / (rp - rm) * ((rp - a * λ / 2) * Ip - (rm - a * λ / 2) * Im)
 end
@@ -568,7 +596,7 @@ See [`r_potential(x)`](@ref) for an implementation of \$\\mathcal{R}(r)\$.
 - `νr` : Radial emission direction (Only necessary for case 1&2 geodesics)
 - `λ`  : Reduced azimuthal angular momentum
 """
-function Iϕ_w_I0_terms(metric::Kerr{T}, rs, τ, roots, νr, λ) where T
+function Iϕ_w_I0_terms(metric::Kerr{T}, rs, τ, roots, νr, λ) where {T}
     numreals = sum(_isreal2.(roots))
 
     if numreals == 4 #case 2
@@ -576,7 +604,7 @@ function Iϕ_w_I0_terms(metric::Kerr{T}, rs, τ, roots, νr, λ) where T
     elseif numreals == 2 #case3
         return Iϕ_w_I0_terms_case3(metric, rs, τ, roots, λ)
     end #case4
-        return Iϕ_w_I0_terms_case4(metric, rs, τ, roots, λ)
+    return Iϕ_w_I0_terms_case4(metric, rs, τ, roots, λ)
 end
 
 function Iϕ_w_I0_terms_case2(metric::Kerr{T}, rs, τ, roots::NTuple{4}, νr, λ) where {T}
@@ -601,8 +629,8 @@ function Iϕ_w_I0_terms_case2(metric::Kerr{T}, rs, τ, roots::NTuple{4}, νr, λ
     Πp_s = coef_p * JacobiElliptic.Pi(rp3 * r41 / (rp4 * r31), asin(x2_s), k)
     Πm_s = coef_m * JacobiElliptic.Pi(rm3 * r41 / (rm4 * r31), asin(x2_s), k)
 
-    Ip = - τ / rp3 - (-1)^νr*Πp_s
-    Im = - τ / rm3 - (-1)^νr*Πm_s
+    Ip = -τ / rp3 - (-1)^νr * Πp_s
+    Im = -τ / rm3 - (-1)^νr * Πm_s
 
     return -2a / (rp - rm) * ((rp - a * λ / 2) * Ip - (rm - a * λ / 2) * Im)
 end
@@ -642,8 +670,12 @@ function Iϕ_w_I0_terms_case3(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) whe
     R1p_s = R1(αp, φ_s, k3)
     R1m_s = R1(αm, φ_s, k3)
 
-    Ip = -inv(B * rp2 + A * rp1) * ((B + A) * τ + 2 * r21 * √(A * B) / (B * rp2 - A * rp1) * (- R1p_s))
-    Im = -inv(B * rm2 + A * rm1) * ((B + A) * τ + 2 * r21 * √(A * B) / (B * rm2 - A * rm1) * (- R1m_s))
+    Ip =
+        -inv(B * rp2 + A * rp1) *
+        ((B + A) * τ + 2 * r21 * √(A * B) / (B * rp2 - A * rp1) * (-R1p_s))
+    Im =
+        -inv(B * rm2 + A * rm1) *
+        ((B + A) * τ + 2 * r21 * √(A * B) / (B * rm2 - A * rm1) * (-R1m_s))
 
     return -(2a / (rp - rm) * ((rp - a * λ / 2) * Ip - (rm - a * λ / 2) * Im))
 end
@@ -680,8 +712,12 @@ function Iϕ_w_I0_terms_case4(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) whe
     S1p_s = S1(gp, atan(x4_s) + atan(go), k4)
     S1m_s = S1(gm, atan(x4_s) + atan(go), k4)
 
-    Ip = go / (a2 * (1 - go * x4_p)) * (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * ( - S1p_s))
-    Im = go / (a2 * (1 - go * x4_m)) * (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * ( - S1m_s))
+    Ip =
+        go / (a2 * (1 - go * x4_p)) *
+        (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (-S1p_s))
+    Im =
+        go / (a2 * (1 - go * x4_m)) *
+        (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (-S1m_s))
 
     return -2a / (rp - rm) * ((rp - a * λ / 2) * Ip - (rm - a * λ / 2) * Im)
 end
@@ -730,7 +766,8 @@ function It_inf_case2(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
     E_o = √(r31 * r42) * JacobiElliptic.E(asin(x2_o), k)
 
     #equation B37
-    I1_total = log(16 / (r31 + r42)^2) / 2 + r43 * (coef * regularized_Pi(n, asin(inv(√n)), k))# Removed the logarithmic divergence
+    I1_total =
+        log(16 / (r31 + r42)^2) / 2 + r43 * (coef * regularized_Pi(n, asin(inv(√n)), k))# Removed the logarithmic divergence
     #equation B38
     I2_total = r3 - E_o # asymptotic Divergent piece is not included
 
@@ -742,7 +779,12 @@ function It_inf_case2(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
     Ip_total = -Πp_o
     Im_total = -Πm_o
 
-    return -(4 / (rp - rm) * (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) + 2 * I1_total + I2_total)
+    return -(
+        4 / (rp - rm) *
+        (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) +
+        2 * I1_total +
+        I2_total
+    )
 end
 
 function It_inf_case3(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
@@ -781,13 +823,22 @@ function It_inf_case3(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
     Π2_o = ((2 * r21 * √(A * B) / (B2 - A2))^2) * regularized_R2(αo, φ_o, k3)# Divergence is removed, will be added back in the end
 
     # Removed logarithmic divergence
-    I1_total =  Π1_o + log(16 * r21^2 / ((A2 - B2)^2 + 4 * A * B * r21^2)) / 2 
+    I1_total = Π1_o + log(16 * r21^2 / ((A2 - B2)^2 + 4 * A * B * r21^2)) / 2
     # Removed linear divergence
-    I2_total = ( - √(A * B) * (Π2_o)) + (B * r2 + A * r1) / (A + B)
-    Ip_total = -inv(B * rp2 + A * rp1) * ( 2 * r21 * √(A * B) / (B * rp2 - A * rp1) * (R1(αp, φ_o, k3) ))
-    Im_total = -inv(B * rm2 + A * rm1) * ( 2 * r21 * √(A * B) / (B * rm2 - A * rm1) * (R1(αm, φ_o, k3) ))
-
-    return -(4 / (rp - rm) * (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) + 2 * I1_total + I2_total)
+    I2_total = (-√(A * B) * (Π2_o)) + (B * r2 + A * r1) / (A + B)
+    Ip_total =
+        -inv(B * rp2 + A * rp1) *
+        (2 * r21 * √(A * B) / (B * rp2 - A * rp1) * (R1(αp, φ_o, k3)))
+    Im_total =
+        -inv(B * rm2 + A * rm1) *
+        (2 * r21 * √(A * B) / (B * rm2 - A * rm1) * (R1(αm, φ_o, k3)))
+
+    return -(
+        4 / (rp - rm) *
+        (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) +
+        2 * I1_total +
+        I2_total
+    )
 end
 
 function It_inf_case4(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
@@ -821,23 +872,38 @@ function It_inf_case4(metric::Kerr{T}, roots::NTuple{4}, λ) where {T}
     S1m_o = S1(gm, T(π / 2) + atan(go), k4)
 
     Π1_o = 2 / (C + D) * (a2 / go * (1 + go^2)) * regularizedS1(go, T(π / 2) + atan(go), k4) # Divergence is removed, will be added back in the end
-    Π2_o = 2 / (C + D) * (a2 / go * (1 + go^2))^2 * regularizedS2(go, T(π / 2) + atan(go), k4) # Divergence is removed, will be added back in the end
-
-    I1_total = - (Π1_o) + 1/2*log(
-        (16*(1 + go^2 - sqrt((1 + go^2)*(1 + go^2 - k4)))*(1 + go^2 - 
-   k4))/((C + D)^2*((1 + go^2)^2 - k4)*k4*(1 + sqrt(
-   1 - k4/(1 + go^2)))))
+    Π2_o =
+        2 / (C + D) * (a2 / go * (1 + go^2))^2 * regularizedS2(go, T(π / 2) + atan(go), k4) # Divergence is removed, will be added back in the end
+
+    I1_total =
+        -(Π1_o) +
+        1 / 2 * log(
+            (16 * (1 + go^2 - sqrt((1 + go^2) * (1 + go^2 - k4))) * (1 + go^2 - k4)) /
+            ((C + D)^2 * ((1 + go^2)^2 - k4) * k4 * (1 + sqrt(1 - k4 / (1 + go^2)))),
+        )
 
     # Removed linear divergence
-    I2_total = - 2(a2 / go - b1) * (Π1_o) + (Π2_o) -
-        (
-            (16 * a2^4 + (C^2 - D^2)^2 - 8 * (a2^2) * (C^2 + D^2) + 8 * (a2^3) * (C + D - 2 * b1 * go) + 
-            2 * a2 * (C + D) * (-(C - D)^2 + 2 * b1 * (C + D) * go)) / (4 * a2 * (4 * a2^2 - (C + D)^2) * go)
+    I2_total =
+        -2(a2 / go - b1) * (Π1_o) + (Π2_o) - (
+            (
+                16 * a2^4 + (C^2 - D^2)^2 - 8 * (a2^2) * (C^2 + D^2) +
+                8 * (a2^3) * (C + D - 2 * b1 * go) +
+                2 * a2 * (C + D) * (-(C - D)^2 + 2 * b1 * (C + D) * go)
+            ) / (4 * a2 * (4 * a2^2 - (C + D)^2) * go)
         )
-    Ip_total = go / (a2 * (1 - go * x4_p)) * ( - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (S1p_o))
-    Im_total = go / (a2 * (1 - go * x4_m)) * ( - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (S1m_o))
+    Ip_total =
+        go / (a2 * (1 - go * x4_p)) *
+        (-2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (S1p_o))
+    Im_total =
+        go / (a2 * (1 - go * x4_m)) *
+        (-2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (S1m_o))
 
-    return -(4 / (rp - rm) * (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total)  + 2 * I1_total + I2_total)# + (logdiv + lineardiv)
+    return -(
+        4 / (rp - rm) *
+        (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) +
+        2 * I1_total +
+        I2_total
+    )# + (logdiv + lineardiv)
 end
 
 """
@@ -901,17 +967,23 @@ function It_w_I0_terms_case2(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ, νr)
     I1_total = r3 * I0_total + r43 * (-1)^νr * Π1_s# Removed the logarithmic divergence
     #equation B38
     I2_s = √(evalpoly(rs, poly_coefs)) / (rs - r3) - E_s
-    I2_total = - (r1 * r4 + r2 * r3) / 2 * τ + (-1)^νr * I2_s# asymptotic Divergent piece is not included
+    I2_total = -(r1 * r4 + r2 * r3) / 2 * τ + (-1)^νr * I2_s# asymptotic Divergent piece is not included
 
     coef_p = 2 / √(r31 * r42) * r43 / (rp3 * rp4)
     coef_m = 2 / √(r31 * r42) * r43 / (rm3 * rm4)
     Πp_s = coef_p * JacobiElliptic.Pi(rp3 * r41 / (rp4 * r31), asin(x2_s), k)
     Πm_s = coef_m * JacobiElliptic.Pi(rm3 * r41 / (rm4 * r31), asin(x2_s), k)
 
-    Ip_total = - τ / rp3 - (-1)^νr*Πp_s
-    Im_total = - τ / rm3 - (-1)^νr*Πm_s 
+    Ip_total = -τ / rp3 - (-1)^νr * Πp_s
+    Im_total = -τ / rm3 - (-1)^νr * Πm_s
 
-    return (4 / (rp - rm) * (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) + 4 * I0_total + 2 * I1_total + I2_total)
+    return (
+        4 / (rp - rm) *
+        (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) +
+        4 * I0_total +
+        2 * I1_total +
+        I2_total
+    )
 end
 
 function It_w_I0_terms_case3(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) where {T}
@@ -952,13 +1024,26 @@ function It_w_I0_terms_case3(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) wher
 
     I0_total = τ
     # Removed logarithmic divergence
-    I1_total = (B * r2 + A * r1) / (B + A) * I0_total - Π1_s 
+    I1_total = (B * r2 + A * r1) / (B + A) * I0_total - Π1_s
     # Removed linear divergence
-    I2_total = ((((B * r2 + A * r1) / (B + A))^2) * I0_total - 2 * (B * r2 + A * r1) / (B + A) * (-Π1_s) - √(A * B) * ( - Π2_s)) #+ (B * r2 + A * r1) / (A + B)
-    Ip_total = -inv(B * rp2 + A * rp1) * ((B + A) * I0_total + 2 * r21 * √(A * B) / (B * rp2 - A * rp1) * (- R1(αp, φ_s, k3)))
-    Im_total = -inv(B * rm2 + A * rm1) * ((B + A) * I0_total + 2 * r21 * √(A * B) / (B * rm2 - A * rm1) * (- R1(αm, φ_s, k3)))
-
-    return (4 / (rp - rm) * (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) + 4 * I0_total + 2 * I1_total + I2_total)
+    I2_total = (
+        (((B * r2 + A * r1) / (B + A))^2) * I0_total -
+        2 * (B * r2 + A * r1) / (B + A) * (-Π1_s) - √(A * B) * (-Π2_s)
+    ) #+ (B * r2 + A * r1) / (A + B)
+    Ip_total =
+        -inv(B * rp2 + A * rp1) *
+        ((B + A) * I0_total + 2 * r21 * √(A * B) / (B * rp2 - A * rp1) * (-R1(αp, φ_s, k3)))
+    Im_total =
+        -inv(B * rm2 + A * rm1) *
+        ((B + A) * I0_total + 2 * r21 * √(A * B) / (B * rm2 - A * rm1) * (-R1(αm, φ_s, k3)))
+
+    return (
+        4 / (rp - rm) *
+        (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) +
+        4 * I0_total +
+        2 * I1_total +
+        I2_total
+    )
 end
 
 function It_w_I0_terms_case4(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) where {T}
@@ -997,18 +1082,27 @@ function It_w_I0_terms_case4(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) wher
 
     I0_total = τ
     # Removed logarithmic divergence
-    I1_total = (a2 / go - b1) * τ - ( - Π1_s) 
+    I1_total = (a2 / go - b1) * τ - (-Π1_s)
 
     # Removed linear divergence
-    I2_total =
-        ((a2 / go - b1)^2) * τ - 2(a2 / go - b1) * (- Π1_s) + (- Π2_s) 
-    Ip_total = go / (a2 * (1 - go * x4_p)) * (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * ( - S1p_s))
-    Im_total = go / (a2 * (1 - go * x4_m)) * (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * ( - S1m_s))
-
-    return (4 / (rp - rm) * (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) + 4 * I0_total + 2 * I1_total + I2_total)# + (logdiv + lineardiv)
-end
-
-function radial_inf_integrals(met::Kerr{T}, roots::NTuple{4}) where T
+    I2_total = ((a2 / go - b1)^2) * τ - 2(a2 / go - b1) * (-Π1_s) + (-Π2_s)
+    Ip_total =
+        go / (a2 * (1 - go * x4_p)) *
+        (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (-S1p_s))
+    Im_total =
+        go / (a2 * (1 - go * x4_m)) *
+        (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (-S1m_s))
+
+    return (
+        4 / (rp - rm) *
+        (rp * (rp - a * λ / 2) * Ip_total - rm * (rm - a * λ / 2) * Im_total) +
+        4 * I0_total +
+        2 * I1_total +
+        I2_total
+    )# + (logdiv + lineardiv)
+end
+
+function radial_inf_integrals(met::Kerr{T}, roots::NTuple{4}) where {T}
     numreals = sum(_isreal2.(roots))
     if numreals == 4
         I1, I2, Ip, Im = Krang.radial_inf_integrals_case2(met, roots)
@@ -1044,7 +1138,8 @@ function radial_inf_integrals_case2(metric::Kerr{T}, roots::NTuple{4}) where {T}
     E_o = √(r31 * r42) * JacobiElliptic.E(asin(x2_o), k)
 
     #equation B37
-    I1o_m_I0_terms =  log(16 / (r31 + r42)^2) / 2 + r43 * (coef * regularized_Pi(n, asin(inv(√n)), k) )
+    I1o_m_I0_terms =
+        log(16 / (r31 + r42)^2) / 2 + r43 * (coef * regularized_Pi(n, asin(inv(√n)), k))
     #equation B38
     I2o_m_I0_terms = r3 - E_o# asymptotic Divergent piece is not included
 
@@ -1088,11 +1183,15 @@ function radial_inf_integrals_case3(metric::Kerr{T}, roots::NTuple{4}) where {T}
     Π1_o = 2 * r21 * √real(A * B) / (B2 - A2) * regularized_R1(αo, φ_o, k3)
     Π2_o = ((2 * r21 * √(A * B) / (B2 - A2))^2) * regularized_R2(αo, φ_o, k3)
 
-    I1o_m_I0_terms =  Π1_o + log(16 * r21^2 / ((A2 - B2)^2 + 4 * A * B * r21^2)) / 2
+    I1o_m_I0_terms = Π1_o + log(16 * r21^2 / ((A2 - B2)^2 + 4 * A * B * r21^2)) / 2
     # Removed linear divergence
-    I2o_m_I0_terms = - √(A * B) * Π2_o + (B * r2 + A * r1) / (A + B)
-    Ipo_m_I0_terms = -inv(B * rp2 + A * rp1) * (2 * r21 * √(A * B) / (B * rp2 - A * rp1) * (R1(αp, φ_o, k3)))
-    Imo_m_I0_terms = -inv(B * rm2 + A * rm1) * (2 * r21 * √(A * B) / (B * rm2 - A * rm1) * (R1(αm, φ_o, k3)))
+    I2o_m_I0_terms = -√(A * B) * Π2_o + (B * r2 + A * r1) / (A + B)
+    Ipo_m_I0_terms =
+        -inv(B * rp2 + A * rp1) *
+        (2 * r21 * √(A * B) / (B * rp2 - A * rp1) * (R1(αp, φ_o, k3)))
+    Imo_m_I0_terms =
+        -inv(B * rm2 + A * rm1) *
+        (2 * r21 * √(A * B) / (B * rm2 - A * rm1) * (R1(αm, φ_o, k3)))
 
     return I1o_m_I0_terms, I2o_m_I0_terms, Ipo_m_I0_terms, Imo_m_I0_terms
 end
@@ -1123,32 +1222,41 @@ function radial_inf_integrals_case4(metric::Kerr{T}, roots::NTuple{4}) where {T}
     S1p_o = S1(gp, T(π / 2) + atan(go), k4)
     S1m_o = S1(gm, T(π / 2) + atan(go), k4)
     Π1_o = 2 / (C + D) * (a2 / go * (1 + go^2)) * regularizedS1(go, T(π / 2) + atan(go), k4) # Divergence is removed, will be added back in the end
-    Π2_o = 2 / (C + D) * (a2 / go * (1 + go^2))^2 * regularizedS2(go, T(π / 2) + atan(go), k4) # Divergence is removed, will be added back in the end
+    Π2_o =
+        2 / (C + D) * (a2 / go * (1 + go^2))^2 * regularizedS2(go, T(π / 2) + atan(go), k4) # Divergence is removed, will be added back in the end
 
     # Removed logarithmic divergence
-    I1o_m_I0_terms =  -Π1_o + 1/2*log(
-        (16*(1 + go^2 - sqrt((1 + go^2)*(1 + go^2 - k4)))*(1 + go^2 - 
-        k4))/((C + D)^2*((1 + go^2)^2 - k4)*k4*(1 + sqrt(
-        1 - k4/(1 + go^2))))
-    )
+    I1o_m_I0_terms =
+        -Π1_o +
+        1 / 2 * log(
+            (16 * (1 + go^2 - sqrt((1 + go^2) * (1 + go^2 - k4))) * (1 + go^2 - k4)) /
+            ((C + D)^2 * ((1 + go^2)^2 - k4) * k4 * (1 + sqrt(1 - k4 / (1 + go^2)))),
+        )
 
     # Removed linear divergence
-    I2o_m_I0_terms = - 2(a2 / go - b1) * Π1_o + Π2_o -
-        ((16 * a2^4 + 
-            (C^2 - D^2)^2 - 8 * (a2^2) * (C^2 + D^2) + 
-                8 * (a2^3) * (C + D - 2 * b1 * go) + 
-                    2 * a2 * (C + D) * (-(C - D)^2 + 2 * b1 * (C + D) * go)
-        ) / (4 * a2 * (4 * a2^2 - (C + D)^2) * go))
-    Ipo_m_I0_terms = go / (a2 * (1 - go * x4_p)) * ( - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * S1p_o)
-    Imo_m_I0_terms = go / (a2 * (1 - go * x4_m)) * ( - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * S1m_o)
+    I2o_m_I0_terms =
+        -2(a2 / go - b1) * Π1_o + Π2_o - (
+            (
+                16 * a2^4 + (C^2 - D^2)^2 - 8 * (a2^2) * (C^2 + D^2) +
+                8 * (a2^3) * (C + D - 2 * b1 * go) +
+                2 * a2 * (C + D) * (-(C - D)^2 + 2 * b1 * (C + D) * go)
+            ) / (4 * a2 * (4 * a2^2 - (C + D)^2) * go)
+        )
+    Ipo_m_I0_terms =
+        go / (a2 * (1 - go * x4_p)) *
+        (-2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * S1p_o)
+    Imo_m_I0_terms =
+        go / (a2 * (1 - go * x4_m)) *
+        (-2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * S1m_o)
 
     return I1o_m_I0_terms, I2o_m_I0_terms, Ipo_m_I0_terms, Imo_m_I0_terms
 end
 
-function radial_w_I0_terms_integrals(met::Kerr{T}, rs, roots::NTuple{4}, τ, νr) where T
+function radial_w_I0_terms_integrals(met::Kerr{T}, rs, roots::NTuple{4}, τ, νr) where {T}
     numreals = sum(_isreal2.(roots))
     if numreals == 4
-        I1, I2, Ip, Im = Krang.radial_w_I0_terms_integrals_case2(met, rs, real.(roots), τ, νr)
+        I1, I2, Ip, Im =
+            Krang.radial_w_I0_terms_integrals_case2(met, rs, real.(roots), τ, νr)
     elseif numreals == 2
         I1, I2, Ip, Im = Krang.radial_w_I0_terms_integrals_case3(met, rs, roots, τ)
     else
@@ -1160,7 +1268,13 @@ end
 """
 Returns the radial integrals for the case where there are four real roots in the radial potential, with roots outside the horizon.
 """
-function radial_w_I0_terms_integrals_case2(metric::Kerr{T}, rs, roots::NTuple{4}, τ, νr) where {T}
+function radial_w_I0_terms_integrals_case2(
+    metric::Kerr{T},
+    rs,
+    roots::NTuple{4},
+    τ,
+    νr,
+) where {T}
     r1, r2, r3, r4 = roots
     _, r31, r32, r41, r42, _ = _get_root_diffs(roots...)
     r43 = r4 - r3
@@ -1192,7 +1306,7 @@ function radial_w_I0_terms_integrals_case2(metric::Kerr{T}, rs, roots::NTuple{4}
     )
 
     #equation B37
-    I1_total = - r3 * I0_total - r43 * (-1)^νr * Π1_s# Removed the logarithmic divergence
+    I1_total = -r3 * I0_total - r43 * (-1)^νr * Π1_s# Removed the logarithmic divergence
     #equation B38
     I2_s = √abs(evalpoly(rs, poly_coefs)) / (rs - r3) - E_s
     I2_total = (r1 * r4 + r2 * r3) / 2 * τ - (-1)^νr * I2_s# asymptotic Divergent piece is not included
@@ -1202,8 +1316,8 @@ function radial_w_I0_terms_integrals_case2(metric::Kerr{T}, rs, roots::NTuple{4}
     Πp_s = coef_p * JacobiElliptic.Pi(rp3 * r41 / (rp4 * r31), asin(x2_s), k)
     Πm_s = coef_m * JacobiElliptic.Pi(rm3 * r41 / (rm4 * r31), asin(x2_s), k)
 
-    Ip_total = τ / rp3 + (-1)^νr*Πp_s
-    Im_total = τ / rm3 + (-1)^νr*Πm_s 
+    Ip_total = τ / rp3 + (-1)^νr * Πp_s
+    Im_total = τ / rm3 + (-1)^νr * Πm_s
 
     return I1_total, I2_total, Ip_total, Im_total
 
@@ -1212,7 +1326,12 @@ end
 """
 Returns the radial integrals for the case where there are two real roots in the radial potential
 """
-function radial_w_I0_terms_integrals_case3(metric::Kerr{T}, rs, roots::NTuple{4}, τ) where {T}
+function radial_w_I0_terms_integrals_case3(
+    metric::Kerr{T},
+    rs,
+    roots::NTuple{4},
+    τ,
+) where {T}
     r1, r2, _, _ = roots
     r21, r31, r32, r41, r42, _ = _get_root_diffs(roots...)
 
@@ -1247,11 +1366,18 @@ function radial_w_I0_terms_integrals_case3(metric::Kerr{T}, rs, roots::NTuple{4}
 
     I0_total = τ
     # Removed logarithmic divergence
-    I1_total = - (B * r2 + A * r1) / (B + A) * I0_total + Π1_s 
+    I1_total = -(B * r2 + A * r1) / (B + A) * I0_total + Π1_s
     # Removed linear divergence
-    I2_total =  -((((B * r2 + A * r1) / (B + A))^2) * I0_total - 2 * (B * r2 + A * r1) / (B + A) * (-Π1_s) - √(A * B) * (- Π2_s)) 
-    Ip_total =  inv(B * rp2 + A * rp1) * ((B + A) * I0_total + 2 * r21 * √(A * B) / (B * rp2 - A * rp1) * (- R1(αp, φ_s, k3)))
-    Im_total =  inv(B * rm2 + A * rm1) * ((B + A) * I0_total + 2 * r21 * √(A * B) / (B * rm2 - A * rm1) * (- R1(αm, φ_s, k3)))
+    I2_total = -(
+        (((B * r2 + A * r1) / (B + A))^2) * I0_total -
+        2 * (B * r2 + A * r1) / (B + A) * (-Π1_s) - √(A * B) * (-Π2_s)
+    )
+    Ip_total =
+        inv(B * rp2 + A * rp1) *
+        ((B + A) * I0_total + 2 * r21 * √(A * B) / (B * rp2 - A * rp1) * (-R1(αp, φ_s, k3)))
+    Im_total =
+        inv(B * rm2 + A * rm1) *
+        ((B + A) * I0_total + 2 * r21 * √(A * B) / (B * rm2 - A * rm1) * (-R1(αm, φ_s, k3)))
 
     return I1_total, I2_total, Ip_total, Im_total
 end
@@ -1259,7 +1385,12 @@ end
 """
 Returns the radial integrals for the case where there are no real roots in the radial potential
 """
-function radial_w_I0_terms_integrals_case4(metric::Kerr{T}, rs, roots::NTuple{4}, τ) where {T}
+function radial_w_I0_terms_integrals_case4(
+    metric::Kerr{T},
+    rs,
+    roots::NTuple{4},
+    τ,
+) where {T}
     r1, _, _, r4 = roots
     _, r31, r32, r41, r42, _ = _get_root_diffs(roots...)
     a = metric.spin
@@ -1289,12 +1420,16 @@ function radial_w_I0_terms_integrals_case4(metric::Kerr{T}, rs, roots::NTuple{4}
     Π2_s = 2 / (C + D) * (a2 / go * (1 + go^2))^2 * S2(go, atan(x4_s) + atan(go), k4)
 
     # Removed logarithmic divergence
-    I1_total = - (a2 / go - b1) * τ + (- Π1_s) 
+    I1_total = -(a2 / go - b1) * τ + (-Π1_s)
 
     # Removed linear divergence
-    I2_total = -((a2 / go - b1)^2) * τ + 2(a2 / go - b1) * ( - Π1_s) - ( - Π2_s)
-    Ip_total =  -go / (a2 * (1 - go * x4_p)) * (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (- S1p_s))
-    Im_total =  -go / (a2 * (1 - go * x4_m)) * (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (- S1m_s))
+    I2_total = -((a2 / go - b1)^2) * τ + 2(a2 / go - b1) * (-Π1_s) - (-Π2_s)
+    Ip_total =
+        -go / (a2 * (1 - go * x4_p)) *
+        (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (-S1p_s))
+    Im_total =
+        -go / (a2 * (1 - go * x4_m)) *
+        (τ - 2 / (C + D) * ((1 + go^2) / (go * (go + x4_p))) * (-S1m_s))
 
     return I1_total, I2_total, Ip_total, Im_total
 end
@@ -1308,11 +1443,11 @@ Returns the maximum mino time that can be accrued along a ray.
 - `roots` : Roots of the radial potential
 - `I0_inf` : Mino time at infinity
 """
-function total_mino_time(metric::Kerr{T}, roots::NTuple{4}) where T
+function total_mino_time(metric::Kerr{T}, roots::NTuple{4}) where {T}
     numreals = unsafe_trunc(Int, sum((Krang._isreal2.(roots))))
     I0_inf = Ir_inf(metric, roots)
 
-    if numreals == 4 
+    if numreals == 4
         τf = 2I0_inf
     else
         rh = Krang.horizon(metric)
@@ -1371,7 +1506,7 @@ See [`r_potential(x)`](@ref) for an implementation of \$\\mathcal{R}(r)\$.
 - `νr` : Sign of radial velocity direction at emission. This is always positive for case 3 and case 4 geodesics.
 """
 function It(pix::AbstractPixel, rs, τ, νr)
-    metric=pix.metric
+    metric = pix.metric
     λtemp = λ(metric, pix.screen_coordinate[1], pix.θo)
     tempIt_inf = It_inf(pix)
 
@@ -1390,10 +1525,10 @@ function radial_integrals(pix::AbstractPixel, rs, τ, νr)
     met = metric(pix)
     I1_o, I2_o, Ip_o, Im_o = radial_inf_integrals_m_I0_terms(pix)
     I1_s, I2_s, Ip_s, Im_s = radial_w_I0_terms_integrals(met, rs, pix.roots, τ, νr)
-    return τ, I1_o-I1_s, I2_o-I2_s, Ip_o-Ip_s, Im_o-Im_s
+    return τ, I1_o - I1_s, I2_o - I2_s, Ip_o - Ip_s, Im_o - Im_s
 end
 
-function _rs_case1_and_2(pix::AbstractPixel, rh, τ::T)::Tuple{T, Bool, Bool} where {T}
+function _rs_case1_and_2(pix::AbstractPixel, rh, τ::T)::Tuple{T,Bool,Bool} where {T}
     radial_roots = real.(roots(pix))
     _, _, r3, r4 = radial_roots
     _, r31, r32, r41, r42, _ = _get_root_diffs(radial_roots...)
@@ -1415,7 +1550,7 @@ function _rs_case1_and_2(pix::AbstractPixel, rh, τ::T)::Tuple{T, Bool, Bool} wh
     return (r31 * r4 - r3 * sn) / (r31 - sn), X2 > zero(T), true
 end
 
-function _rs_case3(pix::AbstractPixel, rh, τ::T)::Tuple{T, Bool, Bool} where {T}
+function _rs_case3(pix::AbstractPixel, rh, τ::T)::Tuple{T,Bool,Bool} where {T}
     radial_roots = roots(pix)
     r1, r2, _, _ = radial_roots
     r21, r31, r32, r41, r42, _ = _get_root_diffs(radial_roots...)
@@ -1443,7 +1578,7 @@ function _rs_case3(pix::AbstractPixel, rh, τ::T)::Tuple{T, Bool, Bool} where {T
     return real(num / den), X3 > zero(T), true
 end
 
-function _rs_case4(pix::AbstractPixel, rh, τ::T)::Tuple{T, Bool, Bool} where {T}
+function _rs_case4(pix::AbstractPixel, rh, τ::T)::Tuple{T,Bool,Bool} where {T}
     radial_roots = roots(pix)
     r1, _, _, r4 = radial_roots
     root_diffs = _get_root_diffs(radial_roots...)
@@ -1460,10 +1595,10 @@ function _rs_case4(pix::AbstractPixel, rh, τ::T)::Tuple{T, Bool, Bool} where {T
     go = √max((T(4)a2^2 - (C - D)^2) / ((C + D)^2 - T(4)a2^2), zero(T))
     x4_s = (rh + b1) / a2
     coef = 2 / (C + D)
-    Ir_s = coef*JacobiElliptic.F(atan(x4_s) + atan(go), k4)
+    Ir_s = coef * JacobiElliptic.F(atan(x4_s) + atan(go), k4)
 
     fo = I0_inf(pix)
-    τ > (fo-Ir_s) && return zero(T), false, false
+    τ > (fo - Ir_s) && return zero(T), false, false
 
     X4 = (C + D) / T(2) * (fo - τ)
     num = go - JacobiElliptic.sc(X4, k4)
@@ -1520,7 +1655,8 @@ function Gθ(pix::AbstractPixel, θs::T, isindir, n)::Tuple{T,T,T,T,Bool,Bool} w
     if isincone && (isindir != ((signβ > 0) ⊻ (θo > T(π / 2))))
         return minotime, Gs, Go, Ghat, isvortical, false
     end
-    if ((((signβ < 0) ⊻ (θs > T(π / 2))) ⊻ (n % 2 == 1)) && !isincone && !isvortical) || (isvortical && ((θo >= T(π / 2)) ⊻ (θs > T(π / 2))))
+    if ((((signβ < 0) ⊻ (θs > T(π / 2))) ⊻ (n % 2 == 1)) && !isincone && !isvortical) ||
+       (isvortical && ((θo >= T(π / 2)) ⊻ (θs > T(π / 2))))
         return minotime, Gs, Go, Ghat, isvortical, false
     end
 
@@ -1558,7 +1694,10 @@ function Gθ(pix::AbstractPixel, θs::T, isindir, n)::Tuple{T,T,T,T,Bool,Bool} w
     end
 
     νθ = isincone ? (n % 2 == 1) ⊻ (θo > θs) : !isindir ⊻ (θs > T(π / 2))
-    minotime = (isindir ? (n + 1) * Ghat - signβ * Go - (-1)^νθ * Gs : n * Ghat - signβ * Go - (-1)^νθ * Gs) #Sign of Go indicates whether the ray is from the forward cone or the rear cone
+    minotime = (
+        isindir ? (n + 1) * Ghat - signβ * Go - (-1)^νθ * Gs :
+        n * Ghat - signβ * Go - (-1)^νθ * Gs
+    ) #Sign of Go indicates whether the ray is from the forward cone or the rear cone
     return minotime, Gs, Go, Ghat, isvortical, true
 end
 
@@ -1592,14 +1731,23 @@ function Gs(pix::AbstractPixel, τ::T) where {T}
         tempfac = one(T) / √abs(um * a^2)
         Δτtemp = (τ % Ghat + (θo > T(π / 2) ? -one(T) : one(T)) * signβ * Go)
         n = floor(τ / Ghat)
-        Δτ = (θo > T(π / 2) ? -one(T) : one(T)) * abs(argmin(abs, [(-one(T))^n * signβ * (Ghat - Δτtemp), (-one(T))^n * signβ * Δτtemp]))
+        Δτ =
+            (θo > T(π / 2) ? -one(T) : one(T)) * abs(
+                argmin(
+                    abs,
+                    [(-one(T))^n * signβ * (Ghat - Δτtemp), (-one(T))^n * signβ * Δτtemp],
+                ),
+            )
     else
         argo = cos(θo) / √(up)
         k = m
         tempfac = inv(√abs(um * a^2))
         Δτtemp = (τ % Ghat + signβ * Go)
         n = floor(τ / Ghat)
-        Δτ = argmin(abs, [(-one(T))^n * signβ * (Ghat - Δτtemp), (-one(T))^n * signβ * Δτtemp])
+        Δτ = argmin(
+            abs,
+            [(-one(T))^n * signβ * (Ghat - Δτtemp), (-one(T))^n * signβ * Δτtemp],
+        )
     end
 
     return Δτ
@@ -1622,7 +1770,9 @@ function Gϕ(pix::AbstractPixel, θs::T, isindir, n) where {T}
     if isincone && (isindir != ((signβ > zero(T)) ⊻ (θo > T(π / 2))))
         return ans, Gs, Go, Ghat, isvortical, false
     end
-    if ((((signβ < zero(T)) ⊻ (θs > T(π / 2))) ⊻ (n % 2 == 1)) && !isincone && !isvortical) || (isvortical && ((θo >= T(π / 2)) ⊻ (θs > T(π / 2))))
+    if (
+        (((signβ < zero(T)) ⊻ (θs > T(π / 2))) ⊻ (n % 2 == 1)) && !isincone && !isvortical
+    ) || (isvortical && ((θo >= T(π / 2)) ⊻ (θs > T(π / 2))))
         return ans, Gs, Go, Ghat, isvortical, false
     end
 
@@ -1659,7 +1809,10 @@ function Gϕ(pix::AbstractPixel, θs::T, isindir, n) where {T}
     end
 
     νθ = isincone ? (n % 2 == 1) ⊻ (θo > θs) : !isindir ⊻ (θs > T(π / 2))
-    ans = (isindir ? (n + 1) * Ghat - signβ * Go + (νθ ? 1 : -1) * Gs : n * Ghat - signβ * Go + (νθ ? 1 : -1) * Gs) #Sign of Go indicates whether the ray is from the forward cone or the rear cone
+    ans = (
+        isindir ? (n + 1) * Ghat - signβ * Go + (νθ ? 1 : -1) * Gs :
+        n * Ghat - signβ * Go + (νθ ? 1 : -1) * Gs
+    ) #Sign of Go indicates whether the ray is from the forward cone or the rear cone
     return ans, Gs, Go, Ghat, isvortical, true
 end
 
@@ -1679,7 +1832,9 @@ function Gt(pix::AbstractPixel, θs::T, isindir, n) where {T}
     if isincone && (isindir != ((signβ > zero(T)) ⊻ (θo > T(π / 2))))
         return ans, Gs, Go, Ghat, isvortical, false
     end
-    if ((((signβ < zero(T)) ⊻ (θs > T(π / 2))) ⊻ (n % 2 == 1)) && !isincone && !isvortical) || (isvortical && ((θo >= T(π / 2)) ⊻ (θs > T(π / 2))))
+    if (
+        (((signβ < zero(T)) ⊻ (θs > T(π / 2))) ⊻ (n % 2 == 1)) && !isincone && !isvortical
+    ) || (isvortical && ((θo >= T(π / 2)) ⊻ (θs > T(π / 2))))
         return ans, Gs, Go, Ghat, isvortical, false
     end
 
@@ -1702,7 +1857,7 @@ function Gt(pix::AbstractPixel, θs::T, isindir, n) where {T}
             return ans, Gs, Go, Ghat, isvortical, false
         end
         tempfac = √abs(um / a^2)
-        Go *= ((θs > T(π / 2)) ? -1 : 1) 
+        Go *= ((θs > T(π / 2)) ? -1 : 1)
         Gs = ((θs > T(π / 2)) ? -1 : 1) * tempfac * JacobiElliptic.E(asin(√args), k)
 
     else
@@ -1713,11 +1868,16 @@ function Gt(pix::AbstractPixel, θs::T, isindir, n) where {T}
             return ans, Gs, Go, Ghat, isvortical, false
         end
         tempfac = -2 * up * inv(√abs(um * a^2))
-        Gs = tempfac * (JacobiElliptic.E(asin(args), k) - JacobiElliptic.F(asin(args), k)) / (2k)
+        Gs =
+            tempfac * (JacobiElliptic.E(asin(args), k) - JacobiElliptic.F(asin(args), k)) /
+            (2k)
     end
 
     νθ = isincone ? (n % 2 == 1) ⊻ (θo > θs) : !isindir ⊻ (θs > T(π / 2))
-    ans = (isindir ? (n + 1) * Ghat - signβ * Go + (νθ ? 1 : -1) * Gs : n * Ghat - signβ * Go + (νθ ? 1 : -1) * Gs) #Sign of Go indicates whether the ray is from the forward cone or the rear cone
+    ans = (
+        isindir ? (n + 1) * Ghat - signβ * Go + (νθ ? 1 : -1) * Gs :
+        n * Ghat - signβ * Go + (νθ ? 1 : -1) * Gs
+    ) #Sign of Go indicates whether the ray is from the forward cone or the rear cone
     return ans, Gs, Go, Ghat, isvortical, true
 end
 
@@ -1808,9 +1968,9 @@ function _absGto_Gthat(metric::Kerr{T}, θo, η, λ)::NTuple{2,T} where {T}
     a = metric.spin
     Go, Ghat, isvortical = zero(T), zero(T), η < zero(T)
 
-    Δθ = (1 - (η+ λ^2) / a^2) / T(2)
-    up = Δθ + √(Δθ^2 + η/ a^2)
-    um = Δθ - √(Δθ^2 + η/ a^2)
+    Δθ = (1 - (η + λ^2) / a^2) / T(2)
+    up = Δθ + √(Δθ^2 + η / a^2)
+    um = Δθ - √(Δθ^2 + η / a^2)
     m = up / um
     k = m
 
@@ -1825,7 +1985,7 @@ function _absGto_Gthat(metric::Kerr{T}, θo, η, λ)::NTuple{2,T} where {T}
             return Go, Ghat
         end
         tempfac = √abs(um / a^2)
-        Go =  tempfac * JacobiElliptic.E(asin(√argo), k)
+        Go = tempfac * JacobiElliptic.E(asin(√argo), k)
         Ghat = 2tempfac * JacobiElliptic.E(k)
 
     else
@@ -1835,10 +1995,11 @@ function _absGto_Gthat(metric::Kerr{T}, θo, η, λ)::NTuple{2,T} where {T}
             return Go, Ghat
         end
         tempfac = -2 * up * inv(√abs(um * a^2))
-        Go = tempfac * (JacobiElliptic.E(asin(argo), k) - JacobiElliptic.F(asin(argo), k)) / (2k)
+        Go =
+            tempfac * (JacobiElliptic.E(asin(argo), k) - JacobiElliptic.F(asin(argo), k)) /
+            (2k)
         Ghat = tempfac * (JacobiElliptic.E(k) - JacobiElliptic.K(k)) / k
     end
 
     return Go, Ghat
 end
-
diff --git a/src/schemes/RayTrace.jl b/src/schemes/RayTrace.jl
index d31195a..0b16fb7 100644
--- a/src/schemes/RayTrace.jl
+++ b/src/schemes/RayTrace.jl
@@ -2,41 +2,55 @@ export generate_ray, raytrace
 
 struct RayTrace <: AbstractScheme end
 
-function generate_ray!(ray::Vector{Intersection{T}}, pixel::Krang.AbstractPixel, res::Int) where T
-    actual_res = unsafe_trunc(Int, sum((Krang._isreal2.(Krang.roots(pixel))))) == 4 ? res+1 : res 
+function generate_ray!(
+    ray::Vector{Intersection{T}},
+    pixel::Krang.AbstractPixel,
+    res::Int,
+) where {T}
+    actual_res =
+        unsafe_trunc(Int, sum((Krang._isreal2.(Krang.roots(pixel))))) == 4 ? res + 1 : res
     τf = total_mino_time(pixel)
 
-    Δτ = τf/actual_res
+    Δτ = τf / actual_res
     for I in range(1, res)
-        ts, rs, θs, ϕs, νr, νθ, _ = emission_coordinates(pixel, Δτ*I)
+        ts, rs, θs, ϕs, νr, νθ, _ = emission_coordinates(pixel, Δτ * I)
         ϕs = ϕs % T(2π)
         ray[I] = Intersection(ts, rs, θs, ϕs, νr, νθ)
     end
 end
 
-function generate_ray(pixel::AbstractPixel{T}, res::Int) where T
-    ray = Vector{Intersection{T}}(undef,res)#zeros(T, 3, res)
+function generate_ray(pixel::AbstractPixel{T}, res::Int) where {T}
+    ray = Vector{Intersection{T}}(undef, res)#zeros(T, 3, res)
     generate_ray!(ray, pixel, res)
     return ray
 end
 
 #TODO: debug this
-@kernel function _generate_rays!(rays::AbstractArray{Intersection{T}}, pixels::AbstractMatrix{P}, res::Int) where {T, P<:AbstractPixel}
-    (I,J,K) = @index(Global, NTuple)
+@kernel function _generate_rays!(
+    rays::AbstractArray{Intersection{T}},
+    pixels::AbstractMatrix{P},
+    res::Int,
+) where {T,P<:AbstractPixel}
+    (I, J, K) = @index(Global, NTuple)
     k = Int(K)
-    pixel = pixels[I,J]
-    actual_res = unsafe_trunc(Int, sum((Krang._isreal2.(Krang.roots(pixel))))) == 4 ? res+1 : res 
-    τf =  total_mino_time(pixel)
+    pixel = pixels[I, J]
+    actual_res =
+        unsafe_trunc(Int, sum((Krang._isreal2.(Krang.roots(pixel))))) == 4 ? res + 1 : res
+    τf = total_mino_time(pixel)
 
-    Δτ = τf/actual_res
+    Δτ = τf / actual_res
 
-    ts, rs, θs, ϕs, νr, νθ = (2f0, 0f0, 0f0, 0f0, true, true)
-    rays[I,J,K] = Intersection(ts, rs, θs, ϕs, νr, νθ)
-    emission_coordinates(pixel, Δτ*k)
+    ts, rs, θs, ϕs, νr, νθ = (2.0f0, 0.0f0, 0.0f0, 0.0f0, true, true)
+    rays[I, J, K] = Intersection(ts, rs, θs, ϕs, νr, νθ)
+    emission_coordinates(pixel, Δτ * k)
 
 
 end
-function generate_rays(pixels::AbstractMatrix{P}, res::Int; A=Array)  where {P <: AbstractPixel{T}} where T
+function generate_rays(
+    pixels::AbstractMatrix{P},
+    res::Int;
+    A = Array,
+) where {P<:AbstractPixel{T}} where {T}
     dims = (size(pixels)..., res)
     rays = A{Intersection{T}}(undef, dims...)
     backend = get_backend(rays)
@@ -44,79 +58,112 @@ function generate_rays(pixels::AbstractMatrix{P}, res::Int; A=Array)  where {P <
     return rays
 end
 
-function line_intersection(origin::NTuple{3,T}, line_point_2, t_a, t_b, t_c) where T
-    e1 = (t_b[1]-t_a[1], t_b[2]-t_a[2], t_b[3]-t_a[3])
-    e2 = (t_c[1]-t_a[1], t_c[2]-t_a[2], t_c[3]-t_a[3])
-	s = (origin[1] - t_a[1], origin[2] - t_a[2], origin[3] - t_a[3])
-    direction = (line_point_2[1] - origin[1], line_point_2[2] - origin[2], line_point_2[3] - origin[3])
+function line_intersection(origin::NTuple{3,T}, line_point_2, t_a, t_b, t_c) where {T}
+    e1 = (t_b[1] - t_a[1], t_b[2] - t_a[2], t_b[3] - t_a[3])
+    e2 = (t_c[1] - t_a[1], t_c[2] - t_a[2], t_c[3] - t_a[3])
+    s = (origin[1] - t_a[1], origin[2] - t_a[2], origin[3] - t_a[3])
+    direction = (
+        line_point_2[1] - origin[1],
+        line_point_2[2] - origin[2],
+        line_point_2[3] - origin[3],
+    )
     dir_len = sqrt(dot(direction, direction))
     direction = direction ./ dir_len
 
-	ray_cross_e2 = cross(direction, e2)
-	det = dot(e1, ray_cross_e2)
+    ray_cross_e2 = cross(direction, e2)
+    det = dot(e1, ray_cross_e2)
 
     outarr = (zero(T), zero(T), zero(T))
 
-	if ( -eps(T) < det < eps(T)) return false, outarr end#This ray is parallel to this triangle.
+    if (-eps(T) < det < eps(T))
+        return false, outarr
+    end#This ray is parallel to this triangle.
+
+    inv_det = inv(det)
+    u = inv_det * dot(s, ray_cross_e2)
+    if !(one(T) > u > zero(T))
+        return false, outarr
+    end
 
-	inv_det = inv(det)
-	u = inv_det * dot(s, ray_cross_e2)
-	if !(one(T) > u > zero(T)) return false, outarr end
+    s_cross_e1 = cross(s, e1)
+    v = inv_det * dot(direction, s_cross_e1)
+    if ((v < zero(T)) || (u + v > one(T)))
+        return false, outarr
+    end
+    # At this stage we can compute t to find out where the intersection point is on the line.
+    t = inv_det * dot(e2, s_cross_e1)
 
-	s_cross_e1 = cross(s, e1)
-	v = inv_det * dot(direction, s_cross_e1)
-	if ((v < zero(T)) || (u + v > one(T) )) return false, outarr end
-	# At this stage we can compute t to find out where the intersection point is on the line.
-	t = inv_det * dot(e2, s_cross_e1)
+    intersection_point = Tuple(origin .+ (direction .* t))
+    #return Some(intersection_point);
 
-	intersection_point = Tuple(origin .+ (direction .* t))
-	#return Some(intersection_point);
+    if t > zero(T)
+        return t < dir_len, intersection_point
+    end# ray intersection
 
-	if t > zero(T) return t < dir_len, intersection_point end# ray intersection
-        
-	# This means that there is a line intersection but not a ray intersection.
-	return false, outarr
+    # This means that there is a line intersection but not a ray intersection.
+    return false, outarr
 end
 
-function cross(a, b) 
-    return a[2]*b[3] - a[3]*b[2], a[3]*b[1] - a[1]*b[3], a[1]*b[2] - a[2]*b[1]
+function cross(a, b)
+    return a[2] * b[3] - a[3] * b[2], a[3] * b[1] - a[1] * b[3], a[1] * b[2] - a[2] * b[1]
 end
-function dot(a,b)
-    return a[1]*b[1] + a[2]*b[2] + a[3]*b[3]
+function dot(a, b)
+    return a[1] * b[1] + a[2] * b[2] + a[3] * b[3]
 end
 
-function raytrace(pixel::AbstractPixel, mesh::Mesh; res=100)
+function raytrace(pixel::AbstractPixel, mesh::Mesh; res = 100)
     return_trait = returnTrait(mesh.material)
-    return raytrace(return_trait, pixel, mesh; res=res)
+    return raytrace(return_trait, pixel, mesh; res = res)
 end
 
-function raytrace(::AbstractReturnTrait, pix::AbstractPixel{T}, mesh::Mesh; res=100) where {T}
+function raytrace(
+    ::AbstractReturnTrait,
+    pix::AbstractPixel{T},
+    mesh::Mesh;
+    res = 100,
+) where {T}
     observation = zero(T)
-    return _raytrace(observation, pix, mesh; res=res)
+    return _raytrace(observation, pix, mesh; res = res)
 end
 
-function raytrace(::SimplePolarizationTrait, pix::AbstractPixel{T}, mesh::Mesh; res=100) where {T}
+function raytrace(
+    ::SimplePolarizationTrait,
+    pix::AbstractPixel{T},
+    mesh::Mesh;
+    res = 100,
+) where {T}
     observation = StokesParams(zero(T), zero(T), zero(T), zero(T))
-    return _raytrace(observation, pix, mesh; res=res)
+    return _raytrace(observation, pix, mesh; res = res)
 end
 
-function raytrace(camera::AbstractCamera, mesh::Mesh{A}; res=100) where A
+function raytrace(camera::AbstractCamera, mesh::Mesh{A}; res = 100) where {A}
     return_trait = returnTrait(mesh.material)
-    return raytrace(return_trait, camera, mesh; res=res)
+    return raytrace(return_trait, camera, mesh; res = res)
 end
 
-function raytrace(::AbstractReturnTrait, camera::AbstractCamera, mesh::Mesh{A}; res=100) where {A}
-    intersections = Array{typeof(camera.metric).parameters[1]}(undef, size(camera.screen.pixels)...)#zeros(yield(mesh.material), size(camera.screen.pixels))
+function raytrace(
+    ::AbstractReturnTrait,
+    camera::AbstractCamera,
+    mesh::Mesh{A};
+    res = 100,
+) where {A}
+    intersections =
+        Array{typeof(camera.metric).parameters[1]}(undef, size(camera.screen.pixels)...)#zeros(yield(mesh.material), size(camera.screen.pixels))
     for I in CartesianIndices(camera.screen.pixels)
-        intersections[I] = raytrace(camera.screen.pixels[I], mesh; res=res)
+        intersections[I] = raytrace(camera.screen.pixels[I], mesh; res = res)
     end
     return intersections
 end
 
-function raytrace(::SimplePolarizationTrait, camera::AbstractCamera, mesh::Mesh{A}; res=100) where {A}
+function raytrace(
+    ::SimplePolarizationTrait,
+    camera::AbstractCamera,
+    mesh::Mesh{A};
+    res = 100,
+) where {A}
     intersections = Array{StokesParams}(undef, size(camera.screen.pixels)...)#zeros(yield(mesh.material), size(camera.screen.pixels))
     for I in CartesianIndices(camera.screen.pixels)
-        intersections[I] = raytrace(camera.screen.pixels[I], mesh; res=res)
+        intersections[I] = raytrace(camera.screen.pixels[I], mesh; res = res)
     end
     return intersections
 end
@@ -124,30 +171,30 @@ end
 
 
 #https://journals.aps.org/prd/abstract/10.1103/PhysRevD.86.084049
-function t_kerr_schild(metric::Kerr{T}, tBL, rBL) where T
+function t_kerr_schild(metric::Kerr{T}, tBL, rBL) where {T}
     a = metric.spin
-    temp = sqrt(1-a^2)
+    temp = sqrt(1 - a^2)
     rp = 1 + temp
     rm = 1 - temp
     num = rBL - rp
     den = rBL - rm
 
-    term1 = (rp^2+a^2)*log(abs(num/den))/(2temp)
-    term2 = -2*log(2/(rBL-rm))
+    term1 = (rp^2 + a^2) * log(abs(num / den)) / (2temp)
+    term2 = -2 * log(2 / (rBL - rm))
     return tBL + term1 + term2
 end
 
-function ϕ_kerr_schild(metric::Kerr{T}, rBL, ϕBL) where T
+function ϕ_kerr_schild(metric::Kerr{T}, rBL, ϕBL) where {T}
     a = metric.spin
-    temp = sqrt(1-a^2)
+    temp = sqrt(1 - a^2)
     rp = 1 + temp
     rm = 1 - temp
     num = rBL - rp
     den = rBL - rm
 
-    term1 = a/(2temp)*log(abs(num/den))
-    term2 = -atan(a/rBL)
-    return ϕBL -  term1 - term2
+    term1 = a / (2temp) * log(abs(num / den))
+    term2 = -atan(a / rBL)
+    return ϕBL - term1 - term2
 end
 
 
@@ -155,29 +202,32 @@ end
     Transforms from Boyer-Lindquist to Kerr-Schild coordinates
 """
 function boyer_lindquist_to_quasi_cartesian_kerr_schild(metric::Kerr, tBL, rBL, θBL, ϕBL)
-    a = metric.spin
     t = t_kerr_schild(metric, tBL, rBL)
     r = rBL
     θ = θBL
-    ϕ = ϕ_kerr_schild(metric, rBL, ϕBL) 
+    ϕ = ϕ_kerr_schild(metric, rBL, ϕBL)
     sθ = sin(θ)
     sϕ = sin(ϕ)
     cϕ = cos(ϕ)
-    return (t, sθ*(r*cϕ), sθ*(r*sϕ), r*cos(θ))
+    return (t, sθ * (r * cϕ), sθ * (r * sϕ), r * cos(θ))
 end
 
 """    
     Transforms from Boyer-Lindquist to Kerr-Schild coordinates ignoring time
 """
-function boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(metric::Kerr{T}, rBL, θBL, ϕBL) where {T}
+function boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(
+    metric::Kerr{T},
+    rBL,
+    θBL,
+    ϕBL,
+) where {T}
     a = metric.spin
     r = rBL
     θ = θBL
-    ϕ = ϕ_kerr_schild(metric,rBL, ϕBL)
+    ϕ = ϕ_kerr_schild(metric, rBL, ϕBL)
     sθ = sin(θ)
     sϕ = sin(ϕ)
     cϕ = cos(ϕ)
 
-    return r .* (sθ*cϕ, sθ*sϕ, cos(θ))
+    return r .* (sθ * cϕ, sθ * sϕ, cos(θ))
 end
-
diff --git a/src/schemes/schemes.jl b/src/schemes/schemes.jl
index b1e25d9..2df85c1 100644
--- a/src/schemes/schemes.jl
+++ b/src/schemes/schemes.jl
@@ -1,26 +1,30 @@
 export render, render!, render_cpu_threaded!
 abstract type AbstractScheme end
 
-function render(camera::AbstractCamera, scene::Scene)  
+function render(camera::AbstractCamera, scene::Scene)
     return render.(camera.screen.pixels, Ref(scene))
 end
 
-function render(pixel::AbstractPixel, scene::Scene) 
-     render(returnTrait(scene[1].material), pixel, scene)
+function render(pixel::AbstractPixel, scene::Scene)
+    render(returnTrait(scene[1].material), pixel, scene)
 end
 
 function render(::AbstractReturnTrait, pixel::AbstractPixel{T}, scene::Scene) where {T}
     return _render(zero(T), pixel, scene)
 end
 
-function render(::AbstractPolarizationTrait, pixel::AbstractPixel{T}, scene::Scene) where {T}
+function render(
+    ::AbstractPolarizationTrait,
+    pixel::AbstractPixel{T},
+    scene::Scene,
+) where {T}
     return _render(StokesParams(zero(T), zero(T), zero(T), zero(T)), pixel, scene)
 end
 
 function _render(observation, pixel::AbstractPixel{T}, scene::Scene) where {T}
     mesh = scene[1]
 
-    for itr in 1:length(scene)
+    for itr = 1:length(scene)
         mesh = scene[itr]
         observation += raytrace(pixel, mesh)
     end
@@ -32,18 +36,21 @@ end
     @inbounds store[I] = mesh.material(pixels[I], mesh.geometry)
 end
 
-function render!(store, camera::AbstractCamera, scene::Scene)  
+function render!(store, camera::AbstractCamera, scene::Scene)
     backend = get_backend(store)
     @assert backend == get_backend(camera.screen.pixels)
     @assert size(store) == size(camera.screen.pixels)
     mapreduce(
-        mesh -> begin _render!(backend)(store, camera.screen.pixels, mesh, ndrange=size(store)); store end,
+        mesh -> begin
+            _render!(backend)(store, camera.screen.pixels, mesh, ndrange = size(store))
+            store
+        end,
         +,
-        scene
+        scene,
     )
 end
 
-function render_cpu_threaded!(store, camera::AbstractCamera, scene::Scene)  
+function render_cpu_threaded!(store, camera::AbstractCamera, scene::Scene)
     @assert size(store) == size(camera.screen.pixels)
     mapreduce(
         mesh -> begin
@@ -53,6 +60,6 @@ function render_cpu_threaded!(store, camera::AbstractCamera, scene::Scene)
             store
         end,
         +,
-        scene
+        scene,
     )
 end
diff --git a/test/camera_tests.jl b/test/camera_tests.jl
index 49e7124..fbe24ac 100644
--- a/test/camera_tests.jl
+++ b/test/camera_tests.jl
@@ -1,22 +1,18 @@
 @testset "Camera" begin
     met = Krang.Kerr(0.0)
-    α = √27*cos(π/4)
-    β = √27*sin(π/4)
-    θo = π/4
+    α = √27 * cos(π / 4)
+    β = √27 * sin(π / 4)
+    θo = π / 4
     crit_roots = (-6.0 + 0im, 0.0 + 0im, 3.0 + 0im, 3.0 + 0im)
     @testset "Abstract Camera" begin
 
         @testset "Abstract Pixel" begin
             struct TestPixel{T} <: Krang.AbstractPixel{T}
                 metric::Kerr{T}
-                screen_coordinate::NTuple{2, T}
+                screen_coordinate::NTuple{2,T}
                 θo::T
                 function TestPixel(met::Kerr{T}, α, β, θo) where {T}
-                    new{T}(
-                        met,
-                        (α, β), 
-                        θo
-                    )
+                    new{T}(met, (α, β), θo)
                 end
             end
             abspix = TestPixel(met, α, β, θo)
@@ -52,4 +48,4 @@
         end
 
     end
-end
\ No newline at end of file
+end
diff --git a/test/emission_coordinates_tests.jl b/test/emission_coordinates_tests.jl
index cb1981a..9196de7 100644
--- a/test/emission_coordinates_tests.jl
+++ b/test/emission_coordinates_tests.jl
@@ -5,16 +5,28 @@
         met = Krang.Kerr(a)
         α = 10.0
         β = 1.0
-        θo = π/2
-        @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+        θo = π / 2
+        @testset "$pixtype" for (pixtype, pix) in [
+            ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+            (
+                "Cached Slow Light Intensity Pixel",
+                Krang.SlowLightIntensityPixel(met, α, β, θo),
+            ),
+        ]
 
-            @test !emission_radius(pix, π/2, true, 2)[5]
+            @test !emission_radius(pix, π / 2, true, 2)[5]
 
             @testset "Case 2" begin
                 α = 10.0
                 β = 10.0
                 θo = π / 4
-                @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+                @testset "$pixtype" for (pixtype, pix) in [
+                    ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                    (
+                        "Cached Slow Light Intensity Pixel",
+                        Krang.SlowLightIntensityPixel(met, α, β, θo),
+                    ),
+                ]
 
                     ηcase1 = η(met, α, β, θo)
                     λcase1 = λ(met, α, θo)
@@ -32,7 +44,13 @@
                 β = 1.0
                 θo = π / 4
 
-                @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+                @testset "$pixtype" for (pixtype, pix) in [
+                    ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                    (
+                        "Cached Slow Light Intensity Pixel",
+                        Krang.SlowLightIntensityPixel(met, α, β, θo),
+                    ),
+                ]
 
                     ηcase3 = η(met, α, β, θo)
                     λcase3 = λ(met, α, θo)
@@ -47,7 +65,13 @@
                 α = 0.1
                 β = 0.1
                 θo = π / 4
-                @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+                @testset "$pixtype" for (pixtype, pix) in [
+                    ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                    (
+                        "Cached Slow Light Intensity Pixel",
+                        Krang.SlowLightIntensityPixel(met, α, β, θo),
+                    ),
+                ]
 
                     ηcase4 = η(met, α, β, θo)
                     λcase4 = λ(met, α, θo)
@@ -64,15 +88,34 @@
         a = 0.99
         met = Krang.Kerr(a)
 
-        @test !emission_radius(Krang.SlowLightIntensityPixel(met, 10.0, 1.0, π/2), π/2, true, 2)[5]
+        @test !emission_radius(
+            Krang.SlowLightIntensityPixel(met, 10.0, 1.0, π / 2),
+            π / 2,
+            true,
+            2,
+        )[5]
         @testset "Ordinary Geodesics" begin
-            α = √27*cos(π/4)
-            β = √27*sin(π/4)
+            α = √27 * cos(π / 4)
+            β = √27 * sin(π / 4)
             θo = π / 4
-            @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+            @testset "$pixtype" for (pixtype, pix) in [
+                ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                (
+                    "Cached Slow Light Intensity Pixel",
+                    Krang.SlowLightIntensityPixel(met, α, β, θo),
+                ),
+            ]
 
-                @testset "n:$n" for n in 0:2
-                    @testset "θs:$θs" for θs in [π/5, π/4, π/3, π/2, 2π/3, 3π/4, 4π/5]
+                @testset "n:$n" for n = 0:2
+                    @testset "θs:$θs" for θs in [
+                        π / 5,
+                        π / 4,
+                        π / 3,
+                        π / 2,
+                        2π / 3,
+                        3π / 4,
+                        4π / 5,
+                    ]
                         @testset "isindir:$isindir" for isindir in [true, false]
                             ηcase1 = η(met, α, β, θo)
                             λcase1 = λ(met, α, θo)
@@ -82,7 +125,7 @@
                             if issuccess
                                 testθs = Krang.emission_inclination(pix, τ1)[1]
                                 if !isnan(testθs)
-                                    @test  testθs/ θs ≈ 1 atol = 1e-5
+                                    @test testθs / θs ≈ 1 atol = 1e-5
                                 end
                             end
                         end
@@ -95,10 +138,24 @@
             α = 0.1
             β = 0.1
             θo = π / 4
-            @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+            @testset "$pixtype" for (pixtype, pix) in [
+                ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                (
+                    "Cached Slow Light Intensity Pixel",
+                    Krang.SlowLightIntensityPixel(met, α, β, θo),
+                ),
+            ]
 
-                @testset "n:$n" for n in 0:2
-                    @testset "θs:$θs" for θs in [π/5, π/4, π/3, π/2, 2π/3, 3π/4, 4π/5]
+                @testset "n:$n" for n = 0:2
+                    @testset "θs:$θs" for θs in [
+                        π / 5,
+                        π / 4,
+                        π / 3,
+                        π / 2,
+                        2π / 3,
+                        3π / 4,
+                        4π / 5,
+                    ]
                         @testset "isindir:$isindir" for isindir in [true, false]
                             ηcase1 = η(met, α, β, θo)
                             λcase1 = λ(met, α, θo)
@@ -108,7 +165,7 @@
                             if issuccess
                                 testθs = Krang.emission_inclination(pix, τ1)[1]
                                 if !isnan(testθs)
-                                    @test  testθs/ θs ≈ 1 atol = 1e-5
+                                    @test testθs / θs ≈ 1 atol = 1e-5
                                 end
                             end
                         end
@@ -118,9 +175,10 @@
         end
     end
     @testset "Emission Coordinates" begin
-        @testset "α:$α, β:$β, isindir:$isindir" for (α, β, isindir) in [(10.0, 10.0, true), (1.0, -1.0, false)]
-            θo = π/4
-            θs = π/3
+        @testset "α:$α, β:$β, isindir:$isindir" for (α, β, isindir) in
+                                                    [(10.0, 10.0, true), (1.0, -1.0, false)]
+            θo = π / 4
+            θs = π / 3
             a = 0.99
             met = Krang.Kerr(a)
             ηtemp = η(met, α, β, θo)
@@ -131,47 +189,57 @@
             desc = √(Δθ2 + ηtemp / a2)
             up = Δθ + desc
             θturning = acos(√up) * (1 + 1e-10)
-            @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+            @testset "$pixtype" for (pixtype, pix) in [
+                ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                (
+                    "Cached Slow Light Intensity Pixel",
+                    Krang.SlowLightIntensityPixel(met, α, β, θo),
+                ),
+            ]
 
                 τ = Krang.Gθ(pix, θs, isindir, 0)[1]
                 ts, testrs, ϕs, νr, νθ = Krang.emission_coordinates(pix, θs, isindir, 0)
-                testrs2, ϕs2, νr2, νθ2 = Krang.emission_coordinates_fast_light(pix, θs, isindir, 0)
+                testrs2, ϕs2, νr2, νθ2 =
+                    Krang.emission_coordinates_fast_light(pix, θs, isindir, 0)
                 ts3, testrs3, testθs, ϕs3, νr3, νθ3 = Krang.emission_coordinates(pix, τ)
                 testrs4, θs4, ϕs4, νr4, νθ4 = Krang.emission_coordinates_fast_light(pix, τ)
 
                 testτ = Krang.Ir(pix, isindir, testrs)[1]
                 @testset "Consistency between raytracing methods" begin
-                    @test testrs/testrs2 ≈ 1.0 atol = 1e-5
-                    @test ϕs2/ϕs ≈ 1.0 atol = 1e-5
+                    @test testrs / testrs2 ≈ 1.0 atol = 1e-5
+                    @test ϕs2 / ϕs ≈ 1.0 atol = 1e-5
                     @test νr == νr2
                     @test νθ == νθ2
-                    @test testrs/testrs3 ≈ 1.0 atol = 1e-5
-                    @test ϕs/ϕs3 ≈ 1.0 atol = 1e-5
+                    @test testrs / testrs3 ≈ 1.0 atol = 1e-5
+                    @test ϕs / ϕs3 ≈ 1.0 atol = 1e-5
                     @test νr == νr3
                     @test νθ == νθ3
-                    @test testrs/testrs4 ≈ 1.0 atol = 1e-5
-                    @test θs4/θs ≈ 1.0 atol = 1e-5
-                    @test ϕs4/ϕs ≈ 1.0 atol = 1e-5
+                    @test testrs / testrs4 ≈ 1.0 atol = 1e-5
+                    @test θs4 / θs ≈ 1.0 atol = 1e-5
+                    @test ϕs4 / ϕs ≈ 1.0 atol = 1e-5
                     @test νr == νr4
                     @test νθ == νθ4
 
 
-                    @test testθs/θs ≈ 1.0 atol = 1e-5
-                    @test testτ/τ ≈ 1.0 atol = 1e-5
+                    @test testθs / θs ≈ 1.0 atol = 1e-5
+                    @test testτ / τ ≈ 1.0 atol = 1e-5
                 end
 
-                fϕ(r, p) = a * (2r - a * λtemp) * inv((r^2 - 2r + a^2) * √(r_potential(met, ηtemp, λtemp, r)))
-                probϕ = IntegralProblem(fϕ, testrs, Inf; nout=1)
-                solϕ = solve(probϕ, HCubatureJL(); reltol=1e-8, abstol=1e-8)
+                fϕ(r, p) =
+                    a *
+                    (2r - a * λtemp) *
+                    inv((r^2 - 2r + a^2) * √(r_potential(met, ηtemp, λtemp, r)))
+                probϕ = IntegralProblem(fϕ, testrs, Inf; nout = 1)
+                solϕ = solve(probϕ, HCubatureJL(); reltol = 1e-8, abstol = 1e-8)
                 Iϕ = solϕ.u
 
                 gϕ(θ, p) = csc(θ)^2 * inv(√(Krang.θ_potential(met, ηtemp, λtemp, θ)))
-                probϕs = IntegralProblem(gϕ, θs, π / 2; nout=1)
-                solθs = solve(probϕs, HCubatureJL(); reltol=1e-12, abstol=1e-12)
-                probϕo = IntegralProblem(gϕ, θo, π / 2; nout=1)
-                solθo = solve(probϕo, HCubatureJL(); reltol=1e-12, abstol=1e-12)
-                probϕ = IntegralProblem(gϕ, θturning, π / 2; nout=1)
-                solθt = solve(probϕ, HCubatureJL(); reltol=1e-12, abstol=1e-12)
+                probϕs = IntegralProblem(gϕ, θs, π / 2; nout = 1)
+                solθs = solve(probϕs, HCubatureJL(); reltol = 1e-12, abstol = 1e-12)
+                probϕo = IntegralProblem(gϕ, θo, π / 2; nout = 1)
+                solθo = solve(probϕo, HCubatureJL(); reltol = 1e-12, abstol = 1e-12)
+                probϕ = IntegralProblem(gϕ, θturning, π / 2; nout = 1)
+                solθt = solve(probϕ, HCubatureJL(); reltol = 1e-12, abstol = 1e-12)
 
                 Gϕ = 0
                 if isindir
@@ -188,20 +256,23 @@
                     end
                 end
 
-                @test -(Iϕ + λtemp*Gϕ+ π/2)/ϕs ≈ 1.0 atol = 1e-3
+                @test -(Iϕ + λtemp * Gϕ + π / 2) / ϕs ≈ 1.0 atol = 1e-3
 
-                ft(r, p) = ((r^2 * (r^2 - 2r + a^2) + 2r * (r^2 + a^2 - a * λtemp)) * inv((r^2 - 2r + a^2) * √(r_potential(met, ηtemp, λtemp, r))))
-                probt = IntegralProblem(ft, testrs, 1e6; nout=1)
-                solt = solve(probt, HCubatureJL(); reltol=1e-8, abstol=1e-8)
+                ft(r, p) = (
+                    (r^2 * (r^2 - 2r + a^2) + 2r * (r^2 + a^2 - a * λtemp)) *
+                    inv((r^2 - 2r + a^2) * √(r_potential(met, ηtemp, λtemp, r)))
+                )
+                probt = IntegralProblem(ft, testrs, 1e6; nout = 1)
+                solt = solve(probt, HCubatureJL(); reltol = 1e-8, abstol = 1e-8)
                 It = solt.u
 
                 gt(θ, p) = cos(θ)^2 * inv(√(Krang.θ_potential(met, ηtemp, λtemp, θ)))
-                probts = IntegralProblem(gt, θs, π / 2; nout=1)
-                solθs = solve(probts, HCubatureJL(); reltol=1e-12, abstol=1e-12)
-                probto = IntegralProblem(gt, θo, π / 2; nout=1)
-                solθo = solve(probto, HCubatureJL(); reltol=1e-12, abstol=1e-12)
-                probt = IntegralProblem(gt, θturning, π / 2; nout=1)
-                solθt = solve(probt, HCubatureJL(); reltol=1e-12, abstol=1e-12)
+                probts = IntegralProblem(gt, θs, π / 2; nout = 1)
+                solθs = solve(probts, HCubatureJL(); reltol = 1e-12, abstol = 1e-12)
+                probto = IntegralProblem(gt, θo, π / 2; nout = 1)
+                solθo = solve(probto, HCubatureJL(); reltol = 1e-12, abstol = 1e-12)
+                probt = IntegralProblem(gt, θturning, π / 2; nout = 1)
+                solθt = solve(probt, HCubatureJL(); reltol = 1e-12, abstol = 1e-12)
 
                 Gt = 0
                 if isindir
@@ -218,8 +289,8 @@
                     end
                 end
 
-                @test ((It - 1e6 - 2log(1e6)) + a^2*Gt)/ts ≈ 1.0 atol = 1e-3
+                @test ((It - 1e6 - 2log(1e6)) + a^2 * Gt) / ts ≈ 1.0 atol = 1e-3
             end
         end
     end
-end
\ No newline at end of file
+end
diff --git a/test/kerr_misc_tests.jl b/test/kerr_misc_tests.jl
index ac17e7b..684230e 100644
--- a/test/kerr_misc_tests.jl
+++ b/test/kerr_misc_tests.jl
@@ -3,7 +3,7 @@
     ηtemp = η(met, 5.0, 5.0, π / 4)
     λtemp = λ(met, 5.0, π / 4)
     @test ^(-4.0 + 5.0im, 1 / 5) ≈ (-4.0 + 5.0im)^(1 / 5)
-    @test abs(^(-4.0+0im, 1 / 5)) ≈ abs((-4.0 + 0.0im)^(1 / 5))
+    @test abs(^(-4.0 + 0im, 1 / 5)) ≈ abs((-4.0 + 0.0im)^(1 / 5))
     @test Krang._isreal2(1.0 + 0.0im) == true
     @test λtemp == -5.0 * sin(π / 4)
     @test ηtemp == (25.0 - 0.9^2) * cos(π / 4)^2 + 25.0
@@ -12,7 +12,11 @@
 
     ssmet = Kerr(0.0)
     @test r_potential(ssmet, 27.0, 0.0, 5.0) ≈ 5.0 * (5.0 - 3.0) * (5.0 - 3.0) * (5.0 + 6.0)
-    @test maximum(abs, Krang.get_radial_roots(ssmet, 27.0, 0.0) .- (-6.0 + 0im, 0.0 + 0im, 3.0 + 0im, 3.0 + 0im)) ≈ 0.0 atol = 1e-5
+    @test maximum(
+        abs,
+        Krang.get_radial_roots(ssmet, 27.0, 0.0) .-
+        (-6.0 + 0im, 0.0 + 0im, 3.0 + 0im, 3.0 + 0im),
+    ) ≈ 0.0 atol = 1e-5
 
 
     @testset "Radial Integrals 1" begin
@@ -22,7 +26,13 @@
             α = 10.0
             β = 10.0
             θo = π / 4
-            @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+            @testset "$pixtype" for (pixtype, pix) in [
+                ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                (
+                    "Cached Slow Light Intensity Pixel",
+                    Krang.SlowLightIntensityPixel(met, α, β, θo),
+                ),
+            ]
                 ηcase1 = η(met, α, β, θo)
                 λcase1 = λ(met, α, θo)
                 roots = get_radial_roots(met, ηcase1, λcase1)
@@ -35,26 +45,26 @@
                 I0, I1, I2, Ip, Im = Krang.radial_integrals(pix, rs, τ1, true)
                 @testset "I0/Ir" begin
                     f(r, p) = inv(√(r_potential(met, ηcase1, λcase1, r)))
-                    prob = IntegralProblem(f, rs, Inf; nout=1)
-                    sol = solve(prob, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob = IntegralProblem(f, rs, Inf; nout = 1)
+                    sol = solve(prob, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test τ1 / sol.u ≈ 1.0 atol = 1e-5
 
-                    prob = IntegralProblem(f, real(root), Inf; nout=1)
-                    sol = solve(prob, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob = IntegralProblem(f, real(root), Inf; nout = 1)
+                    sol = solve(prob, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test τf / (2sol.u) ≈ 1.0 atol = 1e-5
                 end
                 @testset "I1" begin
                     f1(r, p) = r * inv(√(r_potential(met, ηcase1, λcase1, r)))
-                    prob1 = IntegralProblem(f1, rs, 1e6; nout=1)
-                    sol1 = solve(prob1, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob1 = IntegralProblem(f1, rs, 1e6; nout = 1)
+                    sol1 = solve(prob1, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test I1 / (sol1.u - log(1e6)) ≈ 1.0 atol = 1e-5
                 end
 
                 @testset "I2" begin
                     #Regularized Time            
                     f2(r, p) = r^2 * inv(√(r_potential(met, ηcase1, λcase1, r)))
-                    prob2 = IntegralProblem(f2, rs, 1e6; nout=1)
-                    sol2 = solve(prob2, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob2 = IntegralProblem(f2, rs, 1e6; nout = 1)
+                    sol2 = solve(prob2, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test I2 / (sol2.u - 1e6) ≈ 1.0 atol = 1e-3
                 end
             end
@@ -64,7 +74,13 @@
             α = 1.0
             β = 3.0
             θo = π / 4
-            @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+            @testset "$pixtype" for (pixtype, pix) in [
+                ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                (
+                    "Cached Slow Light Intensity Pixel",
+                    Krang.SlowLightIntensityPixel(met, α, β, θo),
+                ),
+            ]
 
                 ηcase3 = η(met, α, β, θo)
                 λcase3 = λ(met, α, θo)
@@ -77,25 +93,25 @@
                 I0, I1, I2, Ip, Im = Krang.radial_integrals(pix, rs, τ3, true)
                 @testset "I0/Ir" begin
                     f(r, p) = inv(√(r_potential(met, ηcase3, λcase3, r)))
-                    prob = IntegralProblem(f, rs, Inf; nout=1)
-                    sol = solve(prob, HCubatureJL(); reltol=1e-5, abstol=1e-5)
+                    prob = IntegralProblem(f, rs, Inf; nout = 1)
+                    sol = solve(prob, HCubatureJL(); reltol = 1e-5, abstol = 1e-5)
                     @test τ3 / sol.u ≈ 1.0 atol = 1e-5
 
-                    prob = IntegralProblem(f, rh, Inf; nout=1)
-                    sol = solve(prob, HCubatureJL(); reltol=1e-5, abstol=1e-5)
+                    prob = IntegralProblem(f, rh, Inf; nout = 1)
+                    sol = solve(prob, HCubatureJL(); reltol = 1e-5, abstol = 1e-5)
                     @test τf / sol.u ≈ 1.0 atol = 1e-5
                 end
                 @testset "I1" begin
                     f1(r, p) = r * inv(√(r_potential(met, ηcase3, λcase3, r)))
-                    prob1 = IntegralProblem(f1, rs, 1e6; nout=1)
-                    sol1 = solve(prob1, HCubatureJL(); reltol=1e-8, abstol=1e-8)
+                    prob1 = IntegralProblem(f1, rs, 1e6; nout = 1)
+                    sol1 = solve(prob1, HCubatureJL(); reltol = 1e-8, abstol = 1e-8)
                     @test I1 / (sol1.u - log(1e6)) ≈ 1.0 atol = 1e-5
                 end
                 @testset "I2" begin
                     #Regularized Time            
                     f2(r, p) = r^2 * inv(√(r_potential(met, ηcase3, λcase3, r)))
-                    prob2 = IntegralProblem(f2, rs, 1e6; nout=1)
-                    sol2 = solve(prob2, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob2 = IntegralProblem(f2, rs, 1e6; nout = 1)
+                    sol2 = solve(prob2, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test I2 / (sol2.u - 1e6) ≈ 1.0 atol = 1e-5
                 end
             end
@@ -104,7 +120,13 @@
             α = 0.1
             β = 0.1
             θo = π / 4
-            @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+            @testset "$pixtype" for (pixtype, pix) in [
+                ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                (
+                    "Cached Slow Light Intensity Pixel",
+                    Krang.SlowLightIntensityPixel(met, α, β, θo),
+                ),
+            ]
 
                 ηcase4 = η(met, α, β, θo)
                 λcase4 = λ(met, α, θo)
@@ -112,30 +134,30 @@
                 @test sum(Krang._isreal2.(roots)) == 0
                 rs = horizon(met)
                 rh = horizon(met)
-                τ4 = Ir( pix, true, rs)[1]
+                τ4 = Ir(pix, true, rs)[1]
                 τf = total_mino_time(met, roots)
                 I0, I1, I2, Ip, Im = Krang.radial_integrals(pix, rs, τ4, true)
                 @testset "I0/Ir" begin
                     f(r, p) = inv(√(r_potential(met, ηcase4, λcase4, r)))
-                    prob = IntegralProblem(f, rs, Inf; nout=1)
-                    sol = solve(prob, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob = IntegralProblem(f, rs, Inf; nout = 1)
+                    sol = solve(prob, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test τ4 / sol.u ≈ 1.0 atol = 1e-5
 
-                    prob = IntegralProblem(f, rh, Inf; nout=1)
-                    sol = solve(prob, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob = IntegralProblem(f, rh, Inf; nout = 1)
+                    sol = solve(prob, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test τf / sol.u ≈ 1.0 atol = 1e-5
                 end
                 @testset "I1" begin
                     f1(r, p) = r * inv(√(r_potential(met, ηcase4, λcase4, r)))
-                    prob1 = IntegralProblem(f1, rs, 1e10; nout=1)
-                    sol1 = solve(prob1, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob1 = IntegralProblem(f1, rs, 1e10; nout = 1)
+                    sol1 = solve(prob1, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test I1 / (sol1.u - log(1e10)) ≈ 1.0 atol = 1e-5
                 end
                 @testset "I2" begin
                     #Regularized Time            
                     f2(r, p) = r^2 * inv(√(r_potential(met, ηcase4, λcase4, r)))
-                    prob2 = IntegralProblem(f2, rs, 1e6; nout=1)
-                    sol2 = solve(prob2, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob2 = IntegralProblem(f2, rs, 1e6; nout = 1)
+                    sol2 = solve(prob2, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test I2 / (sol2.u - 1e6) ≈ 1.0 atol = 1e-5
                 end
             end
@@ -148,7 +170,13 @@
             α = 10.0
             β = 10.0
             θo = π / 4
-            @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+            @testset "$pixtype" for (pixtype, pix) in [
+                ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                (
+                    "Cached Slow Light Intensity Pixel",
+                    Krang.SlowLightIntensityPixel(met, α, β, θo),
+                ),
+            ]
 
                 ηcase1 = η(met, α, β, θo)
                 λcase1 = λ(met, α, θo)
@@ -157,28 +185,34 @@
                 @test sum(Krang._isreal2.(roots)) == 4
 
                 rs = 1.1 * real(root)
-                τ1 = Ir( pix, true, rs)[1]
+                τ1 = Ir(pix, true, rs)[1]
                 @testset "Ir" begin
                     f(r, p) = inv(√(r_potential(met, ηcase1, λcase1, r)))
-                    prob = IntegralProblem(f, rs, Inf; nout=1)
-                    sol = solve(prob, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob = IntegralProblem(f, rs, Inf; nout = 1)
+                    sol = solve(prob, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test τ1 / sol.u ≈ 1.0 atol = 1e-5
                 end
 
                 @testset "Iϕ" begin
-                    fϕ(r, p) = a * (2r - a * λcase1) * inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase1, λcase1, r)))
-                    probϕ = IntegralProblem(fϕ, rs, Inf; nout=1)
-                    solϕ = solve(probϕ, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    fϕ(r, p) =
+                        a *
+                        (2r - a * λcase1) *
+                        inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase1, λcase1, r)))
+                    probϕ = IntegralProblem(fϕ, rs, Inf; nout = 1)
+                    solϕ = solve(probϕ, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     Iϕ = Krang.Iϕ(pix, rs, τ1, true)
                     @test Iϕ / solϕ.u ≈ 1.0 atol = 1e-5
                 end
 
                 @testset "It" begin
                     #Regularized Time            
-                    ft(r, p) = -((r^2 * (r^2 - 2r + a^2) + 2r * (r^2 + a^2 - a * λcase1)) * inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase1, λcase1, r))))
+                    ft(r, p) = -(
+                        (r^2 * (r^2 - 2r + a^2) + 2r * (r^2 + a^2 - a * λcase1)) *
+                        inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase1, λcase1, r)))
+                    )
 
-                    probt = IntegralProblem(ft, rs, 1e6; nout=1)
-                    solt = solve(probt, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    probt = IntegralProblem(ft, rs, 1e6; nout = 1)
+                    solt = solve(probt, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     It = Krang.It(pix, rs, τ1, true)
                     @test It / (solt.u + 1e6 + 2log(1e6)) ≈ 1.0 atol = 1e-3
                 end
@@ -189,32 +223,44 @@
             α = 1.0
             β = 3.0
             θo = π / 4
-            @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+            @testset "$pixtype" for (pixtype, pix) in [
+                ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                (
+                    "Cached Slow Light Intensity Pixel",
+                    Krang.SlowLightIntensityPixel(met, α, β, θo),
+                ),
+            ]
 
                 ηcase3 = η(met, α, β, θo)
                 λcase3 = λ(met, α, θo)
                 roots = get_radial_roots(met, ηcase3, λcase3)
                 @test sum(Krang._isreal2.(roots)) == 2
                 rs = 1.1horizon(met)
-                τ3 = Ir( pix, true, rs)[1]
+                τ3 = Ir(pix, true, rs)[1]
                 @testset "Ir" begin
                     f(r, p) = inv(√(r_potential(met, ηcase3, λcase3, r)))
-                    prob = IntegralProblem(f, rs, Inf; nout=1)
-                    sol = solve(prob, HCubatureJL(); reltol=1e-5, abstol=1e-5)
+                    prob = IntegralProblem(f, rs, Inf; nout = 1)
+                    sol = solve(prob, HCubatureJL(); reltol = 1e-5, abstol = 1e-5)
                     @test τ3 / sol.u ≈ 1.0 atol = 1e-5
                 end
                 @testset "Iϕ" begin
-                    fϕ(r, p) = a * (2r - a * λcase3) * inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase3, λcase3, r)))
-                    probϕ = IntegralProblem(fϕ, rs, Inf; nout=1)
-                    solϕ = solve(probϕ, HCubatureJL(); reltol=1e-8, abstol=1e-8)
+                    fϕ(r, p) =
+                        a *
+                        (2r - a * λcase3) *
+                        inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase3, λcase3, r)))
+                    probϕ = IntegralProblem(fϕ, rs, Inf; nout = 1)
+                    solϕ = solve(probϕ, HCubatureJL(); reltol = 1e-8, abstol = 1e-8)
                     Iϕ = Krang.Iϕ(pix, rs, τ3, true)
                     @test Iϕ / (solϕ.u) ≈ 1.0 atol = 1e-5
                 end
                 @testset "It" begin
                     #Regularized Time            
-                    ft(r, p) = -((r^2 * (r^2 - 2r + a^2) + 2r * (r^2 + a^2 - a * λcase3)) * inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase3, λcase3, r))))
-                    probt = IntegralProblem(ft, rs, 1e6; nout=1)
-                    solt = solve(probt, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    ft(r, p) = -(
+                        (r^2 * (r^2 - 2r + a^2) + 2r * (r^2 + a^2 - a * λcase3)) *
+                        inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase3, λcase3, r)))
+                    )
+                    probt = IntegralProblem(ft, rs, 1e6; nout = 1)
+                    solt = solve(probt, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     It = Krang.It(pix, rs, τ3, true)
                     #@test It ≈ solt.u + 1e6 + 2log(1e6) atol = 1e-3
                     @test It / (solt.u + 1e6 + 2log(1e6)) ≈ 1.0 atol = 1e-3
@@ -226,7 +272,13 @@
             α = 0.05
             β = 0.1
             θo = π / 4
-            @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+            @testset "$pixtype" for (pixtype, pix) in [
+                ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                (
+                    "Cached Slow Light Intensity Pixel",
+                    Krang.SlowLightIntensityPixel(met, α, β, θo),
+                ),
+            ]
 
                 ηcase4 = η(met, α, β, θo)
                 λcase4 = λ(met, α, θo)
@@ -236,22 +288,28 @@
                 τ4 = Ir(pix, true, rs)[1]
                 @testset "Ir" begin
                     f(r, p) = inv(√(r_potential(met, ηcase4, λcase4, r)))
-                    prob = IntegralProblem(f, rs, Inf; nout=1)
-                    sol = solve(prob, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    prob = IntegralProblem(f, rs, Inf; nout = 1)
+                    sol = solve(prob, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     @test τ4 / sol.u ≈ 1.0 atol = 1e-5
                 end
                 @testset "Iϕ" begin
-                    fϕ(r, p) = a * (2r - a * λcase4) * inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase4, λcase4, r)))
-                    probϕ = IntegralProblem(fϕ, rs, Inf; nout=1)
-                    solϕ = solve(probϕ, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    fϕ(r, p) =
+                        a *
+                        (2r - a * λcase4) *
+                        inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase4, λcase4, r)))
+                    probϕ = IntegralProblem(fϕ, rs, Inf; nout = 1)
+                    solϕ = solve(probϕ, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     Iϕ = Krang.Iϕ(pix, rs, τ4, true)
                     @test Iϕ / solϕ.u ≈ 1.0 atol = 2e-2
                 end
                 @testset "It" begin
                     #Regularized Time            
-                    ft(r, p) = -((r^2 * (r^2 - 2r + a^2) + 2r * (r^2 + a^2 - a * λcase4)) * inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase4, λcase4, r))))
-                    probt = IntegralProblem(ft, rs, 1e6; nout=1)
-                    solt = solve(probt, HCubatureJL(); reltol=1e-10, abstol=1e-10)
+                    ft(r, p) = -(
+                        (r^2 * (r^2 - 2r + a^2) + 2r * (r^2 + a^2 - a * λcase4)) *
+                        inv((r^2 - 2r + a^2) * √(r_potential(met, ηcase4, λcase4, r)))
+                    )
+                    probt = IntegralProblem(ft, rs, 1e6; nout = 1)
+                    solt = solve(probt, HCubatureJL(); reltol = 1e-10, abstol = 1e-10)
                     It = Krang.It(pix, rs, τ4, true)
 
                     @test It / (solt.u + 1e6 + 2log(1e6)) ≈ 1.0 atol = 1e-2
@@ -268,8 +326,14 @@
                 @testset "isindir: $isindir" for isindir in [true, false]
                     β = isindir ? 10.0 : -10.0
                     β *= sign(cos(θo))
-                   
-                    @testset "$pixtype" for (pixtype, pix) in [("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)), ("Cached Slow Light Intensity Pixel", Krang.SlowLightIntensityPixel(met, α, β, θo))] 
+
+                    @testset "$pixtype" for (pixtype, pix) in [
+                        ("Intensity Pixel", Krang.IntensityPixel(met, α, β, θo)),
+                        (
+                            "Cached Slow Light Intensity Pixel",
+                            Krang.SlowLightIntensityPixel(met, α, β, θo),
+                        ),
+                    ]
 
 
                         tempη = η(met, α, β, θo)
@@ -282,14 +346,34 @@
 
                         θturning = acos(√up) * (1 + 1e-10)
                         @testset "Gθ" begin
-                            @testset "θs:$θs" for θs in clamp.((1, 3) .* (acos(-√up) - acos(√up)) ./ 4 .+ acos(√up), 0.0, π)
+                            @testset "θs:$θs" for θs in
+                                                  clamp.(
+                                (1, 3) .* (acos(-√up) - acos(√up)) ./ 4 .+ acos(√up),
+                                0.0,
+                                π,
+                            )
                                 fθ(θ, p) = inv(√(Krang.θ_potential(met, tempη, tempλ, θ)))
-                                probts = IntegralProblem(fθ, θs, π / 2; nout=1)
-                                solθs = solve(probts, HCubatureJL(); reltol=1e-12, abstol=1e-12)
-                                probto = IntegralProblem(fθ, θo, π / 2; nout=1)
-                                solθo = solve(probto, HCubatureJL(); reltol=1e-12, abstol=1e-12)
-                                probt = IntegralProblem(fθ, θturning, π / 2; nout=1)
-                                solθt = solve(probt, HCubatureJL(); reltol=1e-12, abstol=1e-12)
+                                probts = IntegralProblem(fθ, θs, π / 2; nout = 1)
+                                solθs = solve(
+                                    probts,
+                                    HCubatureJL();
+                                    reltol = 1e-12,
+                                    abstol = 1e-12,
+                                )
+                                probto = IntegralProblem(fθ, θo, π / 2; nout = 1)
+                                solθo = solve(
+                                    probto,
+                                    HCubatureJL();
+                                    reltol = 1e-12,
+                                    abstol = 1e-12,
+                                )
+                                probt = IntegralProblem(fθ, θturning, π / 2; nout = 1)
+                                solθt = solve(
+                                    probt,
+                                    HCubatureJL();
+                                    reltol = 1e-12,
+                                    abstol = 1e-12,
+                                )
 
                                 τ = 0
                                 if isindir
@@ -297,7 +381,8 @@
                                 else
                                     τ = abs(solθo.u - solθs.u)
                                 end
-                                tempGθ, tempGs, tempGo, tempGhat, _ = Krang.Gθ(pix, θs, isindir, 0)
+                                tempGθ, tempGs, tempGo, tempGhat, _ =
+                                    Krang.Gθ(pix, θs, isindir, 0)
 
                                 @test tempGθ / τ ≈ 1.0 atol = 1e-3
                                 @test tempGs / solθs.u ≈ 1.0 atol = 1e-3
@@ -308,14 +393,36 @@
                             end
                         end
                         @testset "Gϕ" begin
-                            @testset "θs:$θs" for θs in clamp.((1, 3) .* (acos(-√up) - acos(√up)) ./ 4 .+ acos(√up), 0.0, π)
-                                fθ(θ, p) = csc(θ)^2 * inv(√(Krang.θ_potential(met, tempη, tempλ, θ)))
-                                probts = IntegralProblem(fθ, θs, π / 2; nout=1)
-                                solθs = solve(probts, HCubatureJL(); reltol=1e-12, abstol=1e-12)
-                                probto = IntegralProblem(fθ, θo, π / 2; nout=1)
-                                solθo = solve(probto, HCubatureJL(); reltol=1e-12, abstol=1e-12)
-                                probt = IntegralProblem(fθ, θturning, π / 2; nout=1)
-                                solθt = solve(probt, HCubatureJL(); reltol=1e-12, abstol=1e-12)
+                            @testset "θs:$θs" for θs in
+                                                  clamp.(
+                                (1, 3) .* (acos(-√up) - acos(√up)) ./ 4 .+ acos(√up),
+                                0.0,
+                                π,
+                            )
+                                fθ(θ, p) =
+                                    csc(θ)^2 *
+                                    inv(√(Krang.θ_potential(met, tempη, tempλ, θ)))
+                                probts = IntegralProblem(fθ, θs, π / 2; nout = 1)
+                                solθs = solve(
+                                    probts,
+                                    HCubatureJL();
+                                    reltol = 1e-12,
+                                    abstol = 1e-12,
+                                )
+                                probto = IntegralProblem(fθ, θo, π / 2; nout = 1)
+                                solθo = solve(
+                                    probto,
+                                    HCubatureJL();
+                                    reltol = 1e-12,
+                                    abstol = 1e-12,
+                                )
+                                probt = IntegralProblem(fθ, θturning, π / 2; nout = 1)
+                                solθt = solve(
+                                    probt,
+                                    HCubatureJL();
+                                    reltol = 1e-12,
+                                    abstol = 1e-12,
+                                )
 
                                 τ = 0
                                 if isindir
@@ -323,7 +430,8 @@
                                 else
                                     τ = abs(solθo.u - solθs.u)
                                 end
-                                tempGϕ, tempGs, tempGo, tempGhat,_ = Krang.Gϕ(pix, θs, isindir, 0)
+                                tempGϕ, tempGs, tempGo, tempGhat, _ =
+                                    Krang.Gϕ(pix, θs, isindir, 0)
 
                                 @test tempGϕ / τ ≈ 1.0 atol = 1e-3
                                 @test tempGs / solθs.u ≈ 1.0 atol = 1e-3
@@ -333,14 +441,36 @@
                             end
                         end
                         @testset "Gt" begin
-                            @testset "θs:$θs" for θs in clamp.((1, 3) .* (acos(-√up) - acos(√up)) ./ 4 .+ acos(√up), 0.0, π)
-                                ft(θ, p) = cos(θ)^2 * inv(√(Krang.θ_potential(met, tempη, tempλ, θ)))
-                                probts = IntegralProblem(ft, θs, π / 2; nout=1)
-                                solθs = solve(probts, HCubatureJL(); reltol=1e-12, abstol=1e-12)
-                                probto = IntegralProblem(ft, θo, π / 2; nout=1)
-                                solθo = solve(probto, HCubatureJL(); reltol=1e-12, abstol=1e-12)
-                                probt = IntegralProblem(ft, θturning, π / 2; nout=1)
-                                solθt = solve(probt, HCubatureJL(); reltol=1e-12, abstol=1e-12)
+                            @testset "θs:$θs" for θs in
+                                                  clamp.(
+                                (1, 3) .* (acos(-√up) - acos(√up)) ./ 4 .+ acos(√up),
+                                0.0,
+                                π,
+                            )
+                                ft(θ, p) =
+                                    cos(θ)^2 *
+                                    inv(√(Krang.θ_potential(met, tempη, tempλ, θ)))
+                                probts = IntegralProblem(ft, θs, π / 2; nout = 1)
+                                solθs = solve(
+                                    probts,
+                                    HCubatureJL();
+                                    reltol = 1e-12,
+                                    abstol = 1e-12,
+                                )
+                                probto = IntegralProblem(ft, θo, π / 2; nout = 1)
+                                solθo = solve(
+                                    probto,
+                                    HCubatureJL();
+                                    reltol = 1e-12,
+                                    abstol = 1e-12,
+                                )
+                                probt = IntegralProblem(ft, θturning, π / 2; nout = 1)
+                                solθt = solve(
+                                    probt,
+                                    HCubatureJL();
+                                    reltol = 1e-12,
+                                    abstol = 1e-12,
+                                )
 
                                 τ = 0
                                 if isindir
@@ -349,7 +479,8 @@
                                     τ = abs(solθo.u - solθs.u)
                                 end
 
-                                tempGt, tempGs, tempGo, tempGhat, _ = Krang.Gt(pix, θs, isindir, 0)
+                                tempGt, tempGs, tempGo, tempGhat, _ =
+                                    Krang.Gt(pix, θs, isindir, 0)
 
                                 @test tempGt / τ ≈ 1.0 atol = 1e-3
                                 @test tempGs / solθs.u ≈ 1.0 atol = 1e-3
@@ -362,4 +493,4 @@
             end
         end
     end
-end
\ No newline at end of file
+end
diff --git a/test/mesh_geometry_tests.jl b/test/mesh_geometry_tests.jl
index 0947ee6..d83bd19 100644
--- a/test/mesh_geometry_tests.jl
+++ b/test/mesh_geometry_tests.jl
@@ -1,10 +1,15 @@
 @testset "Mesh Geometry Manipulation" begin
-    bunny = load(download("https://graphics.stanford.edu/~mdfisher/Data/Meshes/bunny.obj", "bunny.obj"))
+    bunny = load(
+        download(
+            "https://graphics.stanford.edu/~mdfisher/Data/Meshes/bunny.obj",
+            "bunny.obj",
+        ),
+    )
     rot = Rotations.AngleAxis(π / 2, 0.0, 0.0, 1.0)
-    points = (Ref(rot) .* (bunny.position .* 150)) .+   Ref([15, 3, 0])
+    points = (Ref(rot) .* (bunny.position .* 150)) .+ Ref([15, 3, 0])
     faces = getfield(bunny, :faces)
     bunny_mesh1 = GeometryBasics.Mesh([Point(x) for x in points], faces)
-    
+
     bunny_mesh = translate(rotate(scale(bunny, 150), π / 2, 0.0, 0.0, 1.0), 15, 3, 0)
-    @test all(map(x-> x[1] == x[2],zip(bunny_mesh.position, bunny_mesh1.position)))
-end
\ No newline at end of file
+    @test all(map(x -> x[1] == x[2], zip(bunny_mesh.position, bunny_mesh1.position)))
+end
diff --git a/test/polarization_tests.jl b/test/polarization_tests.jl
index bb99178..aed6dee 100644
--- a/test/polarization_tests.jl
+++ b/test/polarization_tests.jl
@@ -3,125 +3,195 @@
         metric = Kerr(0.99)
         α = 5.0
         β = 5.0
-        θo = π/4
+        θo = π / 4
         pix = Krang.SlowLightIntensityPixel(metric, α, β, θo)
-        λtemp = λ(metric, 5.0, π/4)
-        ηtemp = η(metric, 5.0, 5.0, π/4)
-        emrs,_,_,issuccess = emission_radius(pix, π/5, true, 0)
-        p_d = p_bl_d(metric, emrs, π/4, ηtemp, λtemp, true, true)
-        met_uu = metric_uu(metric, emrs, π/4)
+        λtemp = λ(metric, 5.0, π / 4)
+        ηtemp = η(metric, 5.0, 5.0, π / 4)
+        emrs, _, _, issuccess = emission_radius(pix, π / 5, true, 0)
+        p_d = p_bl_d(metric, emrs, π / 4, ηtemp, λtemp, true, true)
+        met_uu = metric_uu(metric, emrs, π / 4)
 
         # Null momentum magnitude should be 0.
-        @test (p_d' * met_uu * p_d) ≈ 0.0 atol=1e-10 
+        @test (p_d' * met_uu * p_d) ≈ 0.0 atol = 1e-10
     end
 
     @testset "Tetrad transformation" begin
         metric = Kerr(0.99)
         rs = 4.0
-        θs = π/3
+        θs = π / 3
         jacbz_d_u = jac_bl_d_zamo_u(metric, rs, θs)
         jacbz_u_d = jac_bl_u_zamo_d(metric, rs, θs)
         jaczb_d_u = jac_zamo_d_bl_u(metric, rs, θs)
         jaczb_u_d = jac_zamo_u_bl_d(metric, rs, θs)
-        jacfz_u_d = jac_fluid_u_zamo_d(metric, 0.5, π/4, π/4)
-        jaczf_u_d = jac_fluid_u_zamo_d(metric, -0.5, π/4, π/4)
+        jacfz_u_d = jac_fluid_u_zamo_d(metric, 0.5, π / 4, π / 4)
+        jaczf_u_d = jac_fluid_u_zamo_d(metric, -0.5, π / 4, π / 4)
         trivjacfz_u_d = jac_fluid_u_zamo_d(metric, 0.0, 0.0, 0.0)
         met_uu = metric_uu(metric, rs, θs)
         met_dd = metric_dd(metric, rs, θs)
         minkowski = [-1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
         identity = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
 
-        @test maximum(abs, trivjacfz_u_d .- identity) ≈ 0 atol=1e-10
+        @test maximum(abs, trivjacfz_u_d .- identity) ≈ 0 atol = 1e-10
         @testset "Inversion" begin
-            @test maximum(abs, (jacbz_u_d * jaczb_u_d) .- identity) ≈ 0 atol=1e-10
-            @test maximum(abs, (jacbz_d_u * jaczb_d_u) .- identity) ≈ 0 atol=1e-10
-            @test maximum(abs, (jacfz_u_d * jaczf_u_d) .- identity) ≈ 0 atol=1e-10
+            @test maximum(abs, (jacbz_u_d * jaczb_u_d) .- identity) ≈ 0 atol = 1e-10
+            @test maximum(abs, (jacbz_d_u * jaczb_d_u) .- identity) ≈ 0 atol = 1e-10
+            @test maximum(abs, (jacfz_u_d * jaczf_u_d) .- identity) ≈ 0 atol = 1e-10
         end
 
         # Check that the local metric is Minkowski in the ZAMO basis.
         @testset "Metric transformation" begin
-            @test maximum(abs, (jaczb_u_d * met_uu * jaczb_u_d') .- minkowski) ≈ 0 atol=1e-10
-            @test maximum(abs, (jaczb_d_u * met_dd * jaczb_d_u') .- minkowski) ≈ 0 atol=1e-10
+            @test maximum(abs, (jaczb_u_d * met_uu * jaczb_u_d') .- minkowski) ≈ 0 atol =
+                1e-10
+            @test maximum(abs, (jaczb_d_u * met_dd * jaczb_d_u') .- minkowski) ≈ 0 atol =
+                1e-10
         end
 
         @testset "Polarization routine" begin
             metric = Kerr(0.01)
-            θs = π/2
-            θo = 0.01/180*π
+            θs = π / 2
+            θo = 0.01 / 180 * π
             α = 5.0
             β = 5.0
             pix = Krang.SlowLightIntensityPixel(metric, α, β, θo)
             λtemp = λ(metric, α, θo)
             ηtemp = η(metric, α, β, θo)
-            rs,_,_,issuccess = emission_radius(pix, θs, true, 0)
-
-            magfield = @SVector[0.,0.,1.0]
-            βfluid = @SVector[0.,0.,.0]
-
-            eα, eβ, redshift1, lp1 = Krang.synchrotronPolarization(metric, α, β, rs, θs, θo, magfield, βfluid, true, false)
-            i, redshift2, lp2 = Krang.synchrotronIntensity(metric, α, β, rs, θs, θo, magfield, βfluid, true, false)
-            @test evpa(eα, eβ) ≈ -3π/4 atol = 1e-3
-            @test hypot(eα, eβ)/i ≈ 1 atol = 1e-3
-            @test redshift1/redshift2 ≈ 1.0 atol = 1e-3
-            @test lp1/lp2 ≈ 1.0 atol = 1e-3
+            rs, _, _, issuccess = emission_radius(pix, θs, true, 0)
+
+            magfield = @SVector[0.0, 0.0, 1.0]
+            βfluid = @SVector[0.0, 0.0, 0.0]
+
+            eα, eβ, redshift1, lp1 = Krang.synchrotronPolarization(
+                metric,
+                α,
+                β,
+                rs,
+                θs,
+                θo,
+                magfield,
+                βfluid,
+                true,
+                false,
+            )
+            i, redshift2, lp2 = Krang.synchrotronIntensity(
+                metric,
+                α,
+                β,
+                rs,
+                θs,
+                θo,
+                magfield,
+                βfluid,
+                true,
+                false,
+            )
+            @test evpa(eα, eβ) ≈ -3π / 4 atol = 1e-3
+            @test hypot(eα, eβ) / i ≈ 1 atol = 1e-3
+            @test redshift1 / redshift2 ≈ 1.0 atol = 1e-3
+            @test lp1 / lp2 ≈ 1.0 atol = 1e-3
         end
 
-        @testset "Polarization pixel interface"  begin
+        @testset "Polarization pixel interface" begin
             metric = Kerr(0.01)
-            θs = π/2
-            θo = 0.01/180*π
+            θs = π / 2
+            θo = 0.01 / 180 * π
             α = 5.0
             β = 5.0
             pix = Krang.SlowLightIntensityPixel(metric, α, β, θo)
             λtemp = λ(metric, α, θo)
             ηtemp = η(metric, α, β, θo)
-            rs,_,_,issuccess = emission_radius(pix, θs, true, 0)
+            rs, _, _, issuccess = emission_radius(pix, θs, true, 0)
 
-            magfield = @SVector[0.,0.,1.0]
-            βfluid = @SVector[0.,0.,0.]
+            magfield = @SVector[0.0, 0.0, 1.0]
+            βfluid = @SVector[0.0, 0.0, 0.0]
 
             @testset "Intensity" begin
-                subimgs = (0, )
+                subimgs = (0,)
                 spec = 1.0
                 rpeak = 5.0
                 p1 = 1.0
                 p2 = 1.0
 
-                mat = Krang.ElectronSynchrotronPowerLawIntensity(magfield..., βfluid..., spec, rpeak, p1, p2, subimgs)
+                mat = Krang.ElectronSynchrotronPowerLawIntensity(
+                    magfield...,
+                    βfluid...,
+                    spec,
+                    rpeak,
+                    p1,
+                    p2,
+                    subimgs,
+                )
                 geometry = Krang.ConeGeometry((θs))
 
                 int = render(pix, (Krang.Mesh(geometry, mat),))
 
-                profile(r) = (r/rpeak)^p1/(1+(r/rpeak)^(p1+p2))
-                norm, redshift, lp = Krang.synchrotronIntensity(metric, α, β, rs, θs, θo, magfield, βfluid, true, false)
-                int2 = norm^(1 + spec) * lp * profile(rs)*redshift^(3 + spec)
+                profile(r) = (r / rpeak)^p1 / (1 + (r / rpeak)^(p1 + p2))
+                norm, redshift, lp = Krang.synchrotronIntensity(
+                    metric,
+                    α,
+                    β,
+                    rs,
+                    θs,
+                    θo,
+                    magfield,
+                    βfluid,
+                    true,
+                    false,
+                )
+                int2 = norm^(1 + spec) * lp * profile(rs) * redshift^(3 + spec)
                 @test int ≈ int2 atol = 1e-3
             end
 
             @testset "Polarization" begin
-                @inline profile(r) = 1               
-                subimgs = (0, )
+                @inline profile(r) = 1
+                subimgs = (0,)
                 spec = 1.0
                 rpeak = 5.0
                 p1 = 1.0
                 p2 = 1.0
 
-                mat = Krang.ElectronSynchrotronPowerLawPolarization(magfield..., βfluid..., spec, rpeak, p1, p2, subimgs)
+                mat = Krang.ElectronSynchrotronPowerLawPolarization(
+                    magfield...,
+                    βfluid...,
+                    spec,
+                    rpeak,
+                    p1,
+                    p2,
+                    subimgs,
+                )
                 geometry = Krang.ConeGeometry((θs))
 
                 stokes = render(pix, (Krang.Mesh(geometry, mat),))
 
-                eα, eβ, redshift, lp = Krang.synchrotronPolarization(metric, α, β, rs, θs, θo, magfield, βfluid, true, false)
-                q = (-(eα^2 - eβ^2) )
-                u = (-2 * eα * eβ )
-
-                profile(r) = (r/rpeak)^p1/(1+(r/rpeak)^(p1+p2))
-                int2 = hypot(q,u)^(1 + spec) * lp * profile(rs)*redshift^(3 + spec)
+                eα, eβ, redshift, lp = Krang.synchrotronPolarization(
+                    metric,
+                    α,
+                    β,
+                    rs,
+                    θs,
+                    θo,
+                    magfield,
+                    βfluid,
+                    true,
+                    false,
+                )
+                q = (-(eα^2 - eβ^2))
+                u = (-2 * eα * eβ)
+
+                profile(r) = (r / rpeak)^p1 / (1 + (r / rpeak)^(p1 + p2))
+                int2 = hypot(q, u)^(1 + spec) * lp * profile(rs) * redshift^(3 + spec)
                 @test stokes ≈ Krang.StokesParams(int2, q, u, 0.0)
             end
         end
     end
     @testset "Penrose Walker Constant" begin
-        @test all(penrose_walker(Krang.Kerr(0.0), 1.0, 0.0, [1.0,0.0,0.0,0.0], [0.0,1.0,0.0,0.0]) .== (1.0, 0.0))
+        @test all(
+            penrose_walker(
+                Krang.Kerr(0.0),
+                1.0,
+                0.0,
+                [1.0, 0.0, 0.0, 0.0],
+                [0.0, 1.0, 0.0, 0.0],
+            ) .== (1.0, 0.0),
+        )
     end
-end
\ No newline at end of file
+end
diff --git a/test/raytracer_tests.jl b/test/raytracer_tests.jl
index f769dda..79382fc 100644
--- a/test/raytracer_tests.jl
+++ b/test/raytracer_tests.jl
@@ -2,60 +2,68 @@
     @testset "Boyer-Lindquist to Kerr-Schild transformations" begin
         met = Krang.Kerr(0.99)
         rs = 1e10
-        θs = π/4.0
-        ϕs = π/4.0
+        θs = π / 4.0
+        ϕs = π / 4.0
 
         x = rs * sin(θs) * cos(ϕs)
         y = rs * sin(θs) * sin(ϕs)
         z = rs * cos(θs)
 
         # Boyer-Lindquist and Kerr-Schild coordinates are the same far from the black hole
-        trat, xrat, yrat, zrat =  Krang.boyer_lindquist_to_quasi_cartesian_kerr_schild(met, 1e10, rs, θs, ϕs) ./ (1e10, x, y, z) 
-        @test trat ≈ 1.0 atol = 1e-5 
-        @test xrat ≈ 1.0 atol = 1e-5 
-        @test yrat ≈ 1.0 atol = 1e-5 
-        @test zrat ≈ 1.0 atol = 1e-5 
+        trat, xrat, yrat, zrat =
+            Krang.boyer_lindquist_to_quasi_cartesian_kerr_schild(met, 1e10, rs, θs, ϕs) ./
+            (1e10, x, y, z)
+        @test trat ≈ 1.0 atol = 1e-5
+        @test xrat ≈ 1.0 atol = 1e-5
+        @test yrat ≈ 1.0 atol = 1e-5
+        @test zrat ≈ 1.0 atol = 1e-5
 
         # Boyer-Lindquist and Kerr-Schild coordinates are the same far from the black hole
-        xrat, yrat, zrat =  Krang.boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(met, rs, θs, ϕs) ./ (x, y, z) 
-        @test trat ≈ 1.0 atol = 1e-5 
-        @test xrat ≈ 1.0 atol = 1e-5 
-        @test yrat ≈ 1.0 atol = 1e-5 
-        @test zrat ≈ 1.0 atol = 1e-5 
+        xrat, yrat, zrat =
+            Krang.boyer_lindquist_to_quasi_cartesian_kerr_schild_fast_light(
+                met,
+                rs,
+                θs,
+                ϕs,
+            ) ./ (x, y, z)
+        @test trat ≈ 1.0 atol = 1e-5
+        @test xrat ≈ 1.0 atol = 1e-5
+        @test yrat ≈ 1.0 atol = 1e-5
+        @test zrat ≈ 1.0 atol = 1e-5
 
     end
 
     @testset "Level Set" begin
         metric = Krang.Kerr(0.999)
-        θo = 89/180*π 
+        θo = 89 / 180 * π
         ρmax = 10.0
-        sze = 200 
-            
+        sze = 200
+
         camera = Krang.SlowLightIntensityCamera(metric, θo, -ρmax, ρmax, -ρmax, ρmax, sze)
 
         struct Cone{T} <: Krang.AbstractLevelSetGeometry{T}
             θo::T
         end
-        
-        function (cone::Cone)(x,y,z)
-            r = hypot(x,y,z)
-            return z/r - cos(cone.θo)
+
+        function (cone::Cone)(x, y, z)
+            r = hypot(x, y, z)
+            return z / r - cos(cone.θo)
         end
-        struct XMaterial{N} <: Krang.AbstractMaterial 
+        struct XMaterial{N} <: Krang.AbstractMaterial
             subimgs::NTuple{N,Int}
         end
         Krang.horizon(metric)
-        function (mat::XMaterial)(pix, intersection) 
-            (;rs,) = intersection
-            return horizon(pix.metric)<rs <10.0
+        function (mat::XMaterial)(pix, intersection)
+            (; rs,) = intersection
+            return horizon(pix.metric) < rs < 10.0
         end
 
-        mesh1 = Krang.Mesh(Cone(π/4), XMaterial((0,1,2,3)))
-        intersections1 = raytrace(camera, mesh1, res=1_000)
+        mesh1 = Krang.Mesh(Cone(π / 4), XMaterial((0, 1, 2, 3)))
+        intersections1 = raytrace(camera, mesh1, res = 1_000)
 
-        mesh2 = Krang.Mesh(Krang.ConeGeometry(π/4), XMaterial((0,1,2,3)))
+        mesh2 = Krang.Mesh(Krang.ConeGeometry(π / 4), XMaterial((0, 1, 2, 3)))
         intersections2 = raytrace(camera, mesh2)
 
         sum(intersections1 .- intersections2) .≈ 0.0
     end
-end
\ No newline at end of file
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 51199ae..1b75844 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -13,4 +13,4 @@ include("polarization_tests.jl")
 include("emission_coordinates_tests.jl")
 include("camera_tests.jl")
 include("raytracer_tests.jl")
-include("mesh_geometry_tests.jl")
\ No newline at end of file
+include("mesh_geometry_tests.jl")