diff --git a/examples/Project.toml b/examples/Project.toml index 550dda08..14006a36 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -2,10 +2,8 @@ Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" -ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" -Comrade = "99d987ce-9a1e-4df8-bc0b-1ea019aa547b" DisplayAs = "0b91fe84-8a4c-11e9-3e1d-67c38462b6d6" Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" diff --git a/gpu_examples/Project.toml b/gpu_examples/Project.toml index add6012e..77a7a635 100644 --- a/gpu_examples/Project.toml +++ b/gpu_examples/Project.toml @@ -1,8 +1,6 @@ [deps] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" -FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" -GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Krang = "54806c32-d51a-438d-8447-e0041be2fbfb" Metal = "dde4c033-4e86-420c-a63e-0dd931031962" diff --git a/gpu_examples/coordinate-example-metal.jl b/gpu_examples/coordinate-example-metal.jl index 3dd513dc..0aee7725 100644 --- a/gpu_examples/coordinate-example-metal.jl +++ b/gpu_examples/coordinate-example-metal.jl @@ -4,105 +4,84 @@ using Metal 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())) metric = Krang.Kerr(0.99f0); # Kerr metric with a spin of 0.99 -θo = 45.0f0 * π / 180; # inclination angle of the observer. θo ∈ (0, π) +θo = 45f0 * π / 180; # inclination angle of the observer. θo ∈ (0, π) sze = 400; # resolution of the screen is sze x sze rmin = Krang.horizon(metric); # minimum radius to be ray traced -rmax = 10.0f0; # maximum radius to be ray traced -ρmax = 15.0f0; # horizontal and vertical limits of the screen +rmax = 10f0; # maximum radius to be ray traced +ρmax = 15f0; # horizontal and vertical limits of the screen # 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 = 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 in 1:3 ] # Initialize Camera and Pre-Allocate Memory for data to be plotted -coordinates = (zeros(sze, sze) for _ = 1:3) +coordinates = (zeros(sze, sze) for _ in 1:3) camera = Krang.SlowLightIntensityCamera(metric, θo, -ρmax, ρmax, -ρmax, ρmax, sze); colormaps = (:afmhot, :afmhot, :hsv) colorrange = ((-20, 20), (0, rmax), (0, 2π)) -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 _ = 1:2 # Looping over isindir this way is needed to get Metal to work + isindir = false + for _ in 1:2 # Looping over isindir this way is needed to get Metal to work isindir ⊻= true - ts, rs, θs, ϕs = emission_coordinates(pix, geometry.opening_angle, isindir, n) + ts, rs, θs, ϕs = emission_coordinates(pix, geometry.opening_angle, isindir, n) if rs ≤ rmin || rs ≥ rmax continue end coords = (ts, rs, θs, ϕs) end - return coords + return coords end -θs = (75 * π / 180) # opening angle of the cone -geometry = Krang.ConeGeometry(θs, (0, rmin, rmax)) -coordinate_point.(MtlArray(camera.screen.pixels), Ref(geometry)) - # Draw Function # This function draws 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 = 0:2) - + geometries = (Krang.ConeGeometry(θs, (i, rmin, rmax)) for i in 0:2) + for (i, geometry) in enumerate(geometries) - rendered_scene = - Array(coordinate_point.(MtlArray(camera.screen.pixels), Ref(geometry))) + rendered_scene = Array(coordinate_point.(MtlArray(camera.screen.pixels), Ref(geometry))) for I in CartesianIndices(rendered_scene) times[I] = rendered_scene[I][1] radii[I] = rendered_scene[I][2] azimuths[I] = rendered_scene[I][4] end - coordinates = (times, radii, mod2pi.(azimuths)) - for j = 1:3 - heatmap!( - axes_list[i][j], - coordinates[j], - colormap = colormaps[j], - colorrange = colorrange[j], - ) + coordinates = (times, radii, mod2pi.(azimuths )) + for j in 1:3 + heatmap!(axes_list[i][j], coordinates[j], colormap = colormaps[j], colorrange=colorrange[j]) end end end # Create the animation of Cone of Emission Coordinates -recording = CairoMakie.record( - fig, - "coordinate.gif", - range(0.0f0, 1.0f0π, length = 180), - framerate = 12, -) do θs +recording = CairoMakie.record(fig, "coordinate.gif", range(0f0, 1f0π, length=180), framerate=12) do θs draw!(axes_list, camera, coordinates, rmin, rmax, θs) end + diff --git a/gpu_examples/mino-time-example-metal.jl b/gpu_examples/mino-time-example-metal.jl index b72c76c3..131c2294 100644 --- a/gpu_examples/mino-time-example-metal.jl +++ b/gpu_examples/mino-time-example-metal.jl @@ -4,101 +4,83 @@ using Metal GLMk.Makie.inline!(true) metric = Krang.Kerr(0.99f0); # Kerr spacetime with 0.99 spin -θo = 85.0f0 * π / 180; # Observer inclination angle with respect to spin axis +θo = 85f0 * π / 180; # Observer inclination angle with respect to spin axis sze = 200; # Number of pixels along each axis of the screen -ρmax = 5.0f0; # Size of the screen +ρmax = 5f0; # Size of the screen camera = Krang.SlowLightIntensityCamera(metric, θo, -ρmax, ρmax, -ρmax, ρmax, sze); 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, ), ) GLMk.set_theme!(GLMk.merge(curr_theme, GLMk.theme_latexfonts())) -fig = GLMk.Figure(resolution = (500, 600)); - -recording = - GLMk.record(fig, "raytrace.gif", range(0.1f0, 3.0f0, length = 290), framerate = 15) do τ - GLMk.empty!(fig) - - coordinates = Array( - ((x, τ) -> emission_coordinates(x, τ)[1:4]).(MtlArray(camera.screen.pixels), τ), - ) - time = [i[1] for i in coordinates] - radius = [i[2] for i in coordinates] - inclination = [i[3] for i in coordinates] - azimuth = mod2pi.([i[4] for i in coordinates]) - - 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)", +fig = GLMk.Figure(resolution=(500, 600)); + +recording = GLMk.record(fig, "raytrace.gif", range(0.1f0, 3f0, length=290), framerate=15) do τ + GLMk.empty!(fig) + + coordinates = Array(((x, τ)->emission_coordinates(x, τ)[1:4]).(MtlArray(camera.screen.pixels), τ)) + time = [i[1] for i in coordinates] + radius = [i[2] for i in coordinates] + inclination = [i[3] for i in coordinates] + azimuth = mod2pi.([i[4] for i in coordinates]) + + 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, 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] ) - colormaps = (:afmhot, :afmhot, :afmhot, :hsv) - colorrange = ((-20, 20), (0, 10), (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)) + 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 + # ![image](raytrace.gif) 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)) -lines_to_plot = - Krang.generate_rays(MtlArray(camera.screen.pixels), 5_000; A = MtlArray) |> Array - -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 = 1:4 - for j = 1:4 - GLMk.lines!(ax, lines_to_plot[i, j, :, :]) +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 = Krang.generate_rays(MtlArray(camera.screen.pixels), 5_000; A=MtlArray) |> Array + +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 1:4; + for j in 1:4; + GLMk.lines!(ax, lines_to_plot[i,j,:,:]) end end fig GLMk.save("mino_time_rays.png", fig) -# ![Photons trajectories around Kerr black hole in Boyer-Lindquist Coordinates](mino_time_rays.png) +# ![Photons trajectories around Kerr black hole in Boyer-Lindquist Coordinates](mino_time_rays.png) \ No newline at end of file diff --git a/notebook_examples/mino-time-example.ipynb b/notebook_examples/mino-time-example.ipynb index 0dd396ef..32288771 100644 --- a/notebook_examples/mino-time-example.ipynb +++ b/notebook_examples/mino-time-example.ipynb @@ -181,7 +181,7 @@ "GLMk.ylims!(ax, (-3, 3))\n", "GLMk.zlims!(ax, (-3, 3))\n", "lines_to_plot = []\n", - "lines_to_plot = Krang.generate_ray_cartesian.(camera.screen.pixels, 5_000)\n", + "lines_to_plot = Krang.generate_ray.(camera.screen.pixels, 5_000)\n", "\n", "sphere = GLMk.Sphere(GLMk.Point(0.0,0.0,0.0), horizon(metric))\n", "GLMk.mesh!(ax, sphere, color=:black) # Sphere to represent black hole\n", diff --git a/notebook_examples/raytracing-mesh-example.ipynb b/notebook_examples/raytracing-mesh-example.ipynb index c9710287..96360024 100644 --- a/notebook_examples/raytracing-mesh-example.ipynb +++ b/notebook_examples/raytracing-mesh-example.ipynb @@ -145,7 +145,7 @@ "\n", "GLMk.hidedecorations!(ax)\n", "sphere = GLMk.Sphere(GLMk.Point(0.0,0.0,0.0), horizon(metric)) # Sphere to represent black hole\n", - "lines_to_plot = Krang.generate_ray_cartesian.(camera.screen.pixels, 100) # 100 is the number of steps to take along the ray\n", + "lines_to_plot = Krang.generate_ray.(camera.screen.pixels, 100) # 100 is the number of steps to take along the ray\n", "\n", "img = zeros(sze, sze)\n", "recording = GLMk.record(fig, \"mesh.mp4\", 1:sze*sze, framerate=120) do i\n", diff --git a/src/cameras/IntensityCamera.jl b/src/cameras/IntensityCamera.jl index 651b26e5..ce2c6332 100644 --- a/src/cameras/IntensityCamera.jl +++ b/src/cameras/IntensityCamera.jl @@ -3,14 +3,18 @@ export IntensityCamera $TYPEDEF Intensity Pixel Type. -Each Pixel is associated with a single ray, and caches some information about the ray. +Each Pixel is associated with a single ray. +Pixels cache some information about the ray. """ struct IntensityPixel{T} <: AbstractPixel{T} metric::Kerr{T} + "Bardeen coordiantes if ro is \`Inf`, otherwise longitude and latitude" screen_coordinate::NTuple{2,T} "Radial roots" roots::NTuple{4,Complex{T}} - "Radial antiderivative" + "Radial anti-derivative at observer" + I0_o::T + "Radial anti-derivative at infinity" I0_inf::T "Total possible Mino time" total_mino_time::T @@ -18,8 +22,35 @@ struct IntensityPixel{T} <: AbstractPixel{T} absGθo_Gθhat::NTuple{2,T} "Inclination" θo::T + "radius" + ro::T η::T λ::T + function _IntensityPixel(met::Kerr{T}, tempη, tempλ, θo::T, ro::T, x, y) where {T} + roots = Krang.get_radial_roots(met, tempη, tempλ) + numreals = sum(_isreal2.(roots)) + if (numreals == 2) && (abs(imag(roots[4]) / real(roots[4])) < eps(T)) + roots = (roots[1], roots[4], roots[2], roots[3]) + end + I0_o = Krang.Ir_s(met, ro, roots, true) + I0_inf = Krang.Ir_inf(met, roots) + total_mino_time = numreals == 4 ? I0_o + I0_inf : I0_o - Krang.Ir_s(met, Krang.horizon(met), roots, true) + + new{T}( + met, + (x, y), + roots, + I0_o, + I0_inf, + total_mino_time, + Krang._absGθo_Gθhat(met, θo, tempη, tempλ), + θo, + ro, + tempη, + tempλ, + ) + end + @doc """ IntensityPixel(met::Kerr{T}, α::T, β::T, θo::T) where {T} @@ -42,24 +73,28 @@ struct IntensityPixel{T} <: AbstractPixel{T} function IntensityPixel(met::Kerr{T}, α::T, β::T, θo::T) where {T} tempη = Krang.η(met, α, β, θo) tempλ = Krang.λ(met, α, θo) - roots = Krang.get_radial_roots(met, tempη, tempλ) - numreals = sum(_isreal2.(roots)) - if (numreals == 2) && (abs(imag(roots[4]) / real(roots[4])) < eps(T)) - roots = (roots[1], roots[4], roots[2], roots[3]) - end - I0_inf = Krang.Ir_inf(met, roots) - new{T}( - met, - (α, β), - roots, - I0_inf, - total_mino_time(met, roots), - Krang._absGθo_Gθhat(met, θo, tempη, tempλ), - θo, - tempη, - tempλ, - ) + _IntensityPixel(met, tempη, tempλ, θo, T(Inf), α, β) end + + function IntensityPixel(met::Kerr{T}, p_local_u::Vector{T}, θo::T, ro::T; x=p_local_u[3], y=p_local_u[4]) where {T} + a = met.spin + p_bl_u = jac_bl_u_zamo_d(met, ro, θo) * p_local_u + E, _, _, L = metric_dd(met, ro, θo) * p_bl_u + tempλ = L/E + tempη = (Σ(met, ro, θo)/E*p_bl_u[3])^2 - (a*cos(θo))^2 + (tempλ*cot(θo))^2 + + _IntensityPixel(met, tempη, tempλ, θo, ro, x, y) + end + + function IntensityPixel(met::Kerr{T}, longitude::T, latitude::T, θo::T, ro::T) where {T} + @assert latitude >= -π/2 && latitude <= π/2 "latitude must be in [-π/2, π/2]" + red_α = sin(longitude) + red_β = sin(latitude) + p_local_u = [1, √abs(1-(red_α^2+red_β^2)), -red_α, -red_β] + IntensityPixel(met, p_local_u, θo, ro; x=red_α, y=red_β) + end + + end """ @@ -68,11 +103,11 @@ end Screen made of `IntensityPixel`s. """ struct IntensityScreen{T,A<:AbstractMatrix} <: AbstractScreen - "Minimum and Maximum Bardeen α values" - αrange::NTuple{2,T} + "Minimum and Maximum x coordinate values" + xrange::NTuple{2,T} - "Minimum and Maximum Bardeen β values" - βrange::NTuple{2,T} + "Minimum and Maximum y coordinate values" + yrange::NTuple{2,T} "Data type that stores screen pixel information" pixels::A @@ -93,6 +128,23 @@ struct IntensityScreen{T,A<:AbstractMatrix} <: AbstractScreen screen[I, J] = IntensityPixel(met, α, β, θo) end + @kernel function _generate_screen!( + screen, + met::Kerr{T}, + longitude_min, + longitude_max, + latitude_min, + latitude_max, + θo, + ro, + res, + ) where {T} + I, J = @index(Global, NTuple) + long = latitude_min + (latitude_max - latitude_min) * (T(I) - 1) / (res - 1) + lat = longitude_min + (longitude_max - longitude_min) * (T(J) - 1) / (res - 1) + screen[I, J] = IntensityPixel(met, long, lat, θo, ro) + end + @doc """ IntensityScreen(met::Kerr{T}, αmin::T, αmax::T, βmin::T, βmax::T, θo::T, res::Int; A=Matrix) where {T} @@ -138,7 +190,37 @@ struct IntensityScreen{T,A<:AbstractMatrix} <: AbstractScreen ndrange = (res, res), ) - new{T,A}((αmin, αmax), (βmin, βmax), screen) + new{T,typeof(screen)}((αmin, αmax), (βmin, βmax), screen) + end + function IntensityScreen( + met::Kerr{T}, + longitude_min::T, + longitude_max::T, + latitude_min::T, + latitude_max::T, + θo::T, + ro::T, + res; + A = Matrix, + ) where {T} + screen = A(Matrix{IntensityPixel{T}}(undef, res, res)) + + backend = get_backend(screen) + + _generate_screen!(backend)( + screen, + met, + longitude_min, + longitude_max, + latitude_min, + latitude_max, + θo, + ro, + res, + ndrange = (res, res), + ) + + new{T,typeof(screen)}((longitude_min, longitude_max), (latitude_min, latitude_max), screen) end end @@ -152,8 +234,8 @@ struct IntensityCamera{T,A} <: AbstractCamera metric::Kerr{T} "Data type that stores screen pixel information" screen::IntensityScreen{T,A} - "Observer screen_coordinate" - screen_coordinate::NTuple{2,T} + "camera location" + camera_location::NTuple{2,T} @doc """ IntensityCamera(met::Kerr{T}, θo, αmin, αmax, βmin, βmax, res::Int; A=Matrix) where {T} @@ -183,12 +265,34 @@ struct IntensityCamera{T,A} <: AbstractCamera res::Int; A = Matrix, ) where {T} - new{T,A}( + + screen = IntensityScreen(met, αmin, αmax, βmin, βmax, θo, res; A = A) + new{T,typeof(screen.pixels)}( met, - IntensityScreen(met, αmin, αmax, βmin, βmax, θo, res; A = A), + screen, (T(Inf), θo), ) end + function IntensityCamera( + met::Kerr{T}, + θo, + ro, + longitude_min, + longitude_max, + latitude_min, + latitude_max, + res::Int; + A = Matrix, + ) where {T} + + screen = IntensityScreen(met, longitude_min, longitude_max, latitude_min, latitude_max, θo, ro, res; A = A) + new{T,typeof(screen.pixels)}( + met, + screen, + (ro, θo), + ) + end + end function η(pix::IntensityPixel) @@ -206,15 +310,13 @@ end function inclination(pix::IntensityPixel) return pix.θo end -function I0_inf(pix::IntensityPixel) - return pix.I0_inf +function I0_o(pix::IntensityPixel) + return pix.I0_o 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 5987fdca..ef8a8617 100644 --- a/src/cameras/SlowLightIntensityCamera.jl +++ b/src/cameras/SlowLightIntensityCamera.jl @@ -11,67 +11,53 @@ struct SlowLightIntensityPixel{T} <: AbstractPixel{T} screen_coordinate::NTuple{2,T} "Radial roots" roots::NTuple{4,Complex{T}} - "Radial antiderivative" + "Radial anti-derivative at observer" + I0_o::T + "Radial anti-derivative at infinity" I0_inf::T "Total possible Mino time" total_mino_time::T - "Radial phi antiderivative" + "Radial phi anti-derivative" Iϕ_inf::T - "Radial time antiderivative" + "Radial time anti-derivative" It_inf::T - I1_inf_m_I0_terms::T - I2_inf_m_I0_terms::T - Ip_inf_m_I0_terms::T - Im_inf_m_I0_terms::T - "Angular antiderivative" + I1_o_m_I0_terms::T + I2_o_m_I0_terms::T + Ip_o_m_I0_terms::T + Im_o_m_I0_terms::T + "Angular anti-derivative" absGθo_Gθhat::NTuple{2,T} - "Angular ϕ antiderivative" + "Angular ϕ anti-derivative" absGϕo_Gϕhat::NTuple{2,T} - "Angular t antiderivative" + "Angular t anti-derivative" absGto_Gthat::NTuple{2,T} - "Half orbit of angular t antiderivative" + "Half orbit of angular t anti-derivative" θo::T + ro::T η::T λ::T - @doc """ - SlowLightIntensityPixel(met::Kerr{T}, α::T, β::T, θo::T) where {T} - - Construct a `SlowLightIntensityPixel` object with the given Kerr metric, screen coordinates, and inclination. - - # Arguments - - `met::Kerr{T}`: The Kerr metric. - - `α::T`: The Bardeen α value (screen coordinate). - - `β::T`: The Bardeen β value (screen coordinate). - - `θo::T`: The inclination angle. - # Returns - - A `SlowLightIntensityPixel` object initialized with the given parameters. - - # Details - This function calculates the η and λ values using the provided Kerr metric and screen coordinates. - It then computes the radial roots and adjusts them if necessary. - It also calculates the radial and angular antiderivatives. - Finally, it initializes a `SlowLightIntensityPixel` object with the calculated values and the provided parameters. - """ - function SlowLightIntensityPixel(met::Kerr{T}, α::T, β::T, θo) where {T} - tempη = Krang.η(met, α, β, θo) - tempλ = Krang.λ(met, α, θo) + function _SlowLightIntensityPixel(met::Kerr{T}, tempη, tempλ, θo, ro, x, y) where T roots = Krang.get_radial_roots(met, tempη, tempλ) numreals = sum(_isreal2.(roots)) if (numreals == 2) && (abs(imag(roots[4]) / real(roots[4])) < eps(T)) roots = (roots[1], roots[4], roots[2], roots[3]) end - I1, I2, Ip, Im = radial_inf_integrals(met, roots) + I1, I2, Ip, Im = radial_o_m_I0_terms_integrals(met, ro, roots, true) + I0_o = Krang.Ir_s(met, ro, roots, true) I0_inf = Krang.Ir_inf(met, roots) + total_mino_time = numreals == 4 ? I0_o + I0_inf : I0_o - Krang.Ir_s(met, Krang.horizon(met), roots, true) + new{T}( met, - (α, β), + (x, y), roots, + I0_o, I0_inf, - total_mino_time(met, roots), - Krang.Iϕ_inf(met, roots, tempλ), - Krang.It_inf(met, roots, tempλ), + total_mino_time, + Krang.Iϕ_o_m_I0_terms(met, ro, roots, tempλ, true), + Krang.It_o_m_I0_terms(met, ro, roots, tempλ, true), I1, I2, Ip, @@ -80,10 +66,58 @@ struct SlowLightIntensityPixel{T} <: AbstractPixel{T} Krang._absGϕo_Gϕhat(met, θo, tempη, tempλ), Krang._absGto_Gthat(met, θo, tempη, tempλ), θo, + ro, tempη, tempλ, ) end + + @doc """ + SlowLightIntensityPixel(met::Kerr{T}, α::T, β::T, θo::T) where {T} + + Construct a `SlowLightIntensityPixel` object with the given Kerr metric, screen coordinates, and inclination. + + # Arguments + - `met::Kerr{T}`: The Kerr metric. + - `α::T`: The Bardeen α value (screen coordinate). + - `β::T`: The Bardeen β value (screen coordinate). + - `θo::T`: The inclination angle. + + # Returns + - A `SlowLightIntensityPixel` object initialized with the given parameters. + + # Details + This function calculates the η and λ values using the provided Kerr metric and screen coordinates. + It then computes the radial roots and adjusts them if necessary. + It also calculates the radial and angular anti-derivatives. + Finally, it initializes a `SlowLightIntensityPixel` object with the calculated values and the provided parameters. + """ + function SlowLightIntensityPixel(met::Kerr{T}, α::T, β::T, θo) where {T} + @assert α != 0 "Bardeen α = 0 is not implemented" + @assert β != 0 "Bardeen β = 0 is not implemented" + tempη = Krang.η(met, α, β, θo) + tempλ = Krang.λ(met, α, θo) + _SlowLightIntensityPixel(met, tempη, tempλ, θo, T(Inf), α, β) + end + + function SlowLightIntensityPixel(met::Kerr{T}, p_local_u::Vector{T}, θo::T, ro::T; x=p_local_u[3], y=p_local_u[4]) where {T} + a = met.spin + p_bl_u = jac_bl_u_zamo_d(met, ro, θo) * p_local_u + E, _, _, L = metric_dd(met, ro, θo) * p_bl_u + tempλ = -L/E + tempη = (Σ(met, ro, θo)/E*p_bl_u[3])^2 - (a*cos(θo))^2 + (tempλ*cot(θo))^2 + + _SlowLightIntensityPixel(met, tempη, tempλ, θo, ro, x, y) + end + + function SlowLightIntensityPixel(met::Kerr{T}, longitude::T, latitude::T, θo::T, ro::T) where {T} + @assert longitude >= -π/2 && longitude <= π/2 "Longitude must be in [-π/2, π/2]" + red_α = sin(longitude) + red_β = sin(latitude) + p_local_u = [1, √abs(1-(red_α^2+red_β^2)), -red_α, -red_β] + SlowLightIntensityPixel(met, p_local_u, θo, ro; x=red_α, y=red_β) + end + end """ @@ -116,6 +150,23 @@ struct SlowLightIntensityScreen{T,A<:AbstractMatrix} <: AbstractScreen β = βmin + (βmax - βmin) * (T(J) - 1) / (res - 1) screen[I, J] = SlowLightIntensityPixel(met, α, β, θo) end + + @kernel function _generate_screen!( + screen, + met::Kerr{T}, + longitude_min, + longitude_max, + latitude_min, + latitude_max, + θo, + ro, + res, + ) where {T} + I, J = @index(Global, NTuple) + long = longitude_min + (longitude_max - longitude_min) * (T(I) - 1) / (res - 1) + lat = latitude_min + (latitude_max - latitude_min) * (T(J) - 1) / (res - 1) + screen[I, J] = SlowLightIntensityPixel(met, long, lat, θo, ro) + end @doc """ SlowLightIntensityScreen(met::Kerr{T}, αmin, αmax, βmin, βmax, θo, res; A=Matrix) where {T} @@ -162,6 +213,37 @@ struct SlowLightIntensityScreen{T,A<:AbstractMatrix} <: AbstractScreen new{T,typeof(screen)}((αmin, αmax), (βmin, βmax), screen) end + function SlowLightIntensityScreen( + met::Kerr{T}, + longitude_min, + longitude_max, + latitude_min, + latitude_max, + θo, + ro, + res; + A = Matrix, + ) where {T} + screen = A(Matrix{SlowLightIntensityPixel{T}}(undef, res, res)) + + backend = get_backend(screen) + + _generate_screen!(backend)( + screen, + met, + longitude_min, + longitude_max, + latitude_min, + latitude_max, + θo, + ro, + res, + ndrange = (res, res), + ) + + new{T,typeof(screen)}((longitude_min, longitude_max), (latitude_min, latitude_max), screen) + end + end """ @@ -207,6 +289,21 @@ struct SlowLightIntensityCamera{T,A} <: AbstractCamera screen = SlowLightIntensityScreen(met, αmin, αmax, βmin, βmax, θo, res; A) new{T,typeof(screen.pixels)}(met, screen, (T(Inf), θo)) end + + function SlowLightIntensityCamera( + met::Kerr{T}, + θo, + ro, + longitude_min, + longitude_max, + latitude_min, + latitude_max, + res; + A = Matrix, + ) where {T} + screen = SlowLightIntensityScreen(met, longitude_min, longitude_max, latitude_min, latitude_max, θo, ro, res; A) + new{T,typeof(screen.pixels)}(met, screen, (ro, θo)) + end end function η(pix::SlowLightIntensityPixel) @@ -233,23 +330,23 @@ 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 +function I1_o_m_I0_terms(pix::SlowLightIntensityPixel) + return pix.I1_o_m_I0_terms end -function I2_inf_m_I0_terms(pix::SlowLightIntensityPixel) - return pix.I2_inf_m_I0_terms +function I2_o_m_I0_terms(pix::SlowLightIntensityPixel) + return pix.I2_o_m_I0_terms end -function Ip_inf_m_I0_terms(pix::SlowLightIntensityPixel) - return pix.Ip_inf_m_I0_terms +function Ip_o_m_I0_terms(pix::SlowLightIntensityPixel) + return pix.Ip_o_m_I0_terms end -function Im_inf_m_I0_terms(pix::SlowLightIntensityPixel) - return pix.Im_inf_m_I0_terms +function Im_o_m_I0_terms(pix::SlowLightIntensityPixel) + return pix.Im_o_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) +function radial_o_integrals_m_I0_terms(pix::SlowLightIntensityPixel) + return I1_o_m_I0_terms(pix), + I2_o_m_I0_terms(pix), + Ip_o_m_I0_terms(pix), + Im_o_m_I0_terms(pix) end function Iϕ_inf(pix::SlowLightIntensityPixel) return pix.Iϕ_inf diff --git a/src/cameras/camera_types.jl b/src/cameras/camera_types.jl index d5a8677c..efdaf688 100644 --- a/src/cameras/camera_types.jl +++ b/src/cameras/camera_types.jl @@ -76,6 +76,8 @@ Calculate the η value for a given pixel. - The η value of the pixel. """ function η(pix::AbstractPixel) + @warn "η is not stored in pixel, calculating η under infinite observer assumption instead" + met = metric(pix) α, β = screen_coordinate(pix) θo = inclination(pix) @@ -94,6 +96,7 @@ Calculate the λ value for a given pixel. - The λ value of the pixel. """ function λ(pix::AbstractPixel) + @warn "λ is not stored in pixel, calculating λ under infinite observer assumption instead" met = metric(pix) α, _ = screen_coordinate(pix) θo = inclination(pix) @@ -130,6 +133,25 @@ function I0_inf(pix::AbstractPixel) return Ir_inf(metric(pix), roots(pix)) end + +""" + I0_o(pix::AbstractPixel) + +Calculate the I0 anti-derivative value for a given pixel at radius ro. + +# Arguments +- `pix::AbstractPixel`: The pixel of a screen. + +# Returns +- The I0 anti-derivative of the pixel. +""" +function I0_o(pix::AbstractPixel) + if isinf(pix.ro) + return I0_inf(pix) + end + return Ir_s(metric(pix), pix.ro, roots(pix), true) +end + """ total_mino_time(pix::AbstractPixel) @@ -171,9 +193,27 @@ Calculate the I1 infinity minus I0 terms for a given pixel. - The I1 infinity minus I0 terms of the pixel. """ function I1_inf_m_I0_terms(pix::AbstractPixel) + @warn "I1_inf_m_I0_terms are not stored in pixel, calculating I1_inf_integrals under infinite observer assumption instead" return radial_inf_integrals(metric(pix), roots(pix))[1] end +""" + I1_o_m_I0_terms(pix::AbstractPixel) + +Calculate the I1 anti-derivative minus I0 terms for a given pixel at radius ro. + +# Arguments +- `pix::AbstractPixel`: The pixel of a screen. + +# Returns +- The I1 anti-derivative minus I0 terms of the pixel. +""" +function I1_o_m_I0_terms(pix::AbstractPixel) + @warn "I1_o_m_I0_terms are not stored in pixel, calculating I1_o_integrals under infinite observer assumption instead" + return radial_o_integrals(metric(pix), roots(pix))[1] +end + + """ I2_inf_m_I0_terms(pix::AbstractPixel) @@ -186,8 +226,27 @@ Calculate the I2 infinity minus I0 terms for a given pixel. - The I2 infinity minus I0 terms of the pixel. """ function I2_inf_m_I0_terms(pix::AbstractPixel) + @warn "I2_inf_m_I0_terms are not stored in pixel, calculating I2_inf_integrals under infinite observer assumption instead" return radial_inf_integrals(metric(pix), roots(pix))[2] end + +""" + I2_o_m_I0_terms(pix::AbstractPixel) + +Calculate the I2 anti-derivative minus I0 terms for a given pixel at radius ro. + +# Arguments +- `pix::AbstractPixel`: The pixel of a screen. + +# Returns +- The I2 anti-derivative minus I0 terms of the pixel. +""" +function I2_o_m_I0_terms(pix::AbstractPixel) + @warn "I2_o_m_I0_terms are not stored in pixel, calculating I2_o_integrals under infinite observer assumption instead" + return radial_o_integrals(metric(pix), roots(pix))[2] +end + + """ Ip_inf_m_I0_terms(pix::AbstractPixel) @@ -203,6 +262,22 @@ function Ip_inf_m_I0_terms(pix::AbstractPixel) return radial_inf_integrals(metric(pix), roots(pix))[3] end +""" + Ip_o_m_I0_terms(pix::AbstractPixel) + +Calculate the Ip anti-derivative minus I0 terms for a given pixel at radius ro. + +# Arguments +- `pix::AbstractPixel`: The pixel of a screen. + +# Returns +- The Ip anti-derivative minus I0 terms of the pixel. +""" +function Ip_o_m_I0_terms(pix::AbstractPixel) + return radial_o_integrals(metric(pix), roots(pix))[3] +end + + """ Im_inf_m_I0_terms(pix::AbstractPixel) @@ -218,6 +293,22 @@ function Im_inf_m_I0_terms(pix::AbstractPixel) return radial_inf_integrals(metric(pix), roots(pix))[4] end +""" + Im_o_m_I0_terms(pix::AbstractPixel) + +Calculate the Im anti-derivative minus I0 terms for a given pixel at radius ro. + +# Arguments +- `pix::AbstractPixel`: The pixel of a screen. + +# Returns +- The Im anti-derivative minus I0 terms of the pixel. +""" +function Im_o_m_I0_terms(pix::AbstractPixel) + return radial_o_integrals(metric(pix), roots(pix))[4] +end + + """ radial_inf_integrals_m_I0_terms(pix::AbstractPixel) @@ -233,6 +324,22 @@ function radial_inf_integrals_m_I0_terms(pix::AbstractPixel) return radial_inf_integrals(metric(pix), roots(pix)) end +""" + radial_o_integrals(pix::AbstractPixel) + +Calculate the radial anti-derivative minus I0 terms for a given pixel at radius ro. + +# Arguments +- `pix::AbstractPixel`: The pixel of a screen. + +# Returns +- The radial anti-derivative minus I0 terms of the pixel. +""" +function radial_o_integrals(pix::AbstractPixel) + return radial_o_integrals_m_I0_terms(metric(pix), roots(pix)) +end + + """ Iϕ_inf(pix::AbstractPixel) diff --git a/src/materials/physicsUtils.jl b/src/materials/physicsUtils.jl index 59fac68f..a436baad 100644 --- a/src/materials/physicsUtils.jl +++ b/src/materials/physicsUtils.jl @@ -2,117 +2,11 @@ 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 + evpa + ##---------------------------------------------------------------------------------------------------------------------- #Polarization stuff ##---------------------------------------------------------------------------------------------------------------------- -""" - Returns the momentum form in the Boyer-Lindquist basis. -""" -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, η, λ, θ)), - λ, - ] -end - -""" -Returns the Jacobian which converts a Boyer-Lindquist covector to ZAMO basis. -""" -function jac_zamo_d_bl_u(metric::Kerr{T}, r, θ) where {T} - a = metric.spin - # coords = {t, r, θ, ϕ} - Σt = Σ(metric, r, θ) - Δt = Δ(metric, r) - At = A(metric, r, θ) - - return @SMatrix [# Eq 3.1 1972ApJ...178..347B - √(At / (Σt * Δt)) zero(T) zero(T) 2*a*r/√(At * Δt * Σt) - zero(T) √(Δt / Σt) zero(T) zero(T) - zero(T) zero(T) zero(T) √(Σt / At)*csc(θ) - zero(T) zero(T) -inv(√Σt) zero(T) - ] -end - -""" -Jacobian which converts ZAMO covector to a Boyer-Lindquist basis -""" -function jac_bl_d_zamo_u(metric::Kerr{T}, r, θ) where {T} - a = metric.spin - Σt = Σ(metric, r, θ) - Δt = Δ(metric, r) - At = A(metric, r, θ) - - return @SMatrix [ - # coords = {t, r, ϕ, θ} - √((Σt * Δt) / At) zero(T) -2*a*r*sin(θ)/√(At * Σt) zero(T) - zero(T) √(Σt / Δt) zero(T) zero(T) - zero(T) zero(T) zero(T) -√Σt - zero(T) zero(T) √(At / Σt)*sin(θ) zero(T) - ] -end - -""" -Jacobian which converts Boyer-Lindquist vector to a ZAMO basis -""" -function jac_zamo_u_bl_d(metric::Kerr{T}, r, θ) where {T} - a = metric.spin - Σt = Σ(metric, r, θ) - Δt = Δ(metric, r) - At = A(metric, r, θ) - - return @SMatrix [# Eq 3.2 1972ApJ...178..347B - # coords = {t, r, ϕ, θ} - √((Σt * Δt) / At) zero(T) zero(T) zero(T) - zero(T) √(Σt / Δt) zero(T) zero(T) - -(2a * r * sin(θ))/√(At * Σt) zero(T) zero(T) √(At / Σt)*sin(θ) - zero(T) zero(T) -√Σt zero(T) - ] -end - -""" -Jacobian which converts ZAMO vector to a Boyer-Lindquist basis -""" -function jac_bl_u_zamo_d(metric::Kerr{T}, r, θ) where {T} - a = metric.spin - Σt = Σ(metric, r, θ) - Δt = Δ(metric, r) - At = A(metric, r, θ) - - return @SMatrix [ - # coords = {t, r, ϕ, θ} - √(At / (Σt * Δt)) zero(T) zero(T) zero(T) - zero(T) √(Δt / Σt) zero(T) zero(T) - zero(T) zero(T) zero(T) -inv(√Σt) - 2a*r/√(Σt * Δt * At) zero(T) √(Σt / At)*csc(θ) zero(T) - ] -end - -""" -Jacobian which expreases ZAMO vector in the fluid frame -""" -function jac_fluid_u_zamo_d(::Kerr{T}, β, θ, φ) where {T} - γ = inv(√(1 - β^2)) - sinφ = sin(φ) - cosφ = cos(φ) - sinθ = sin(θ) - cosθ = cos(θ) - - return @SMatrix [ - γ -β*γ*cosφ*sinθ -β*γ*sinφ*sinθ -β*γ*cosθ - -β*γ*cosφ*sinθ cosθ^2*cosφ^2+γ*cosφ^2*sinθ^2+sinφ^2 (γ-1)*cosφ*sinθ^2*sinφ (γ-1)*cosθ*cosφ*sinθ - -β*γ*sinθ*sinφ (γ-1)*cosφ*sinθ^2*sinφ cosφ^2+(cosθ^2+γ*sinθ^2)*sinφ^2 (γ-1)*cosθ*sinθ*sinφ - -β*γ*cosθ (γ-1)*cosθ*cosφ*sinθ (γ-1)*cosθ*sinθ*sinφ γ*cosθ^2+sinθ^2 - ] -end """ Returns the Penrose walker constant for a photon with momentum p_u emitted from a fluid particle with momentum f_u. diff --git a/src/metrics/Kerr/emission_coordinates.jl b/src/metrics/Kerr/emission_coordinates.jl index d5477631..054f4d2b 100644 --- a/src/metrics/Kerr/emission_coordinates.jl +++ b/src/metrics/Kerr/emission_coordinates.jl @@ -139,7 +139,7 @@ Emission azimuth for point at Mino time τ whose image appears at screen coordin Gϕtemp, _, _, _ = @inline Gϕ(pix, θs, isindir, n) (isnan(Gϕtemp) || !isfinite(Gϕtemp)) && return Iϕ - return -(Iϕ + λtemp * Gϕtemp + T(π / 2)) + return (Iϕ + λtemp * Gϕtemp) end """ @@ -226,7 +226,7 @@ Ray trace a point that appears at the screen coordinate (`α`, `β`) for an obse Gϕtemp, _, _, _, _ = @inline Gϕ(pix, θs, isindir, n) - emission_azimuth = -(Iϕ + λtemp * Gϕtemp + T(π / 2)) + emission_azimuth = (Iϕ + λtemp * Gϕtemp) νθ = abs(cos(θs)) < abs(cos(θo)) ? (n % 2 == 1) ⊻ (θo > θs) : !isindir ⊻ (θs > T(π / 2)) return rs, θs, emission_azimuth, νr, νθ, issuccess @@ -292,7 +292,7 @@ coordinate (`α`, `β`) for an observer located at inclination θo. Gϕtemp, _, _, _, _ = @inline Gϕ(pix, θs, isindir, n) Gttemp, _, _, _, _ = @inline Gt(pix, θs, isindir, n) - emission_azimuth = -(Iϕ + λtemp * Gϕtemp + T(π / 2)) + emission_azimuth = (Iϕ + λtemp * Gϕtemp) emission_time_regularized = (zero(T) + It + a^2 * Gttemp) # is θ̇s increasing or decreasing? @@ -345,7 +345,7 @@ Ray trace a point that appears at the screen coordinate (`α`, `β`) for an obse Gϕtemp, _, _, _, _ = @inline Gϕ(pix, θs, isindir, n) Gttemp, _, _, _, _ = @inline Gt(pix, θs, isindir, n) - emission_azimuth = -(Iϕ + λtemp * Gϕtemp + T(π / 2)) + emission_azimuth = (Iϕ + λtemp * Gϕtemp) emission_time_regularized = (It + a^2 * Gttemp) νθ = abs(cos(θs)) < abs(cos(θo)) ? (n % 2 == 1) ⊻ (θo > θs) : !isindir ⊻ (θs > T(π / 2)) diff --git a/src/metrics/Kerr/misc.jl b/src/metrics/Kerr/misc.jl index 7769792c..509dc0d7 100644 --- a/src/metrics/Kerr/misc.jl +++ b/src/metrics/Kerr/misc.jl @@ -11,7 +11,12 @@ export λ, get_radial_roots, Ir, Gθ, - total_mino_time + total_mino_time, + 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 """ Checks if a complex number is real to some tolerance @@ -122,6 +127,110 @@ function S2(α, φ, j) ) + (inv(1 + α^2) + (1 - j) / (1 - j + α^2)) * S1(α, φ, j) end +""" + Returns the momentum form in the Boyer-Lindquist basis. +""" +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, η, λ, θ)), + λ, + ] +end + +""" +Returns the Jacobian which converts a Boyer-Lindquist covector to ZAMO basis. +""" +function jac_zamo_d_bl_u(metric::Kerr{T}, r, θ) where {T} + a = metric.spin + # coords = {t, r, θ, ϕ} + Σt = Σ(metric, r, θ) + Δt = Δ(metric, r) + At = A(metric, r, θ) + + return @SMatrix [# Eq 3.1 1972ApJ...178..347B + √(At / (Σt * Δt)) zero(T) zero(T) 2*a*r/√(At * Δt * Σt) + zero(T) √(Δt / Σt) zero(T) zero(T) + zero(T) zero(T) zero(T) √(Σt / At)*csc(θ) + zero(T) zero(T) -inv(√Σt) zero(T) + ] +end + +""" +Jacobian which converts ZAMO covector to a Boyer-Lindquist basis +""" +function jac_bl_d_zamo_u(metric::Kerr{T}, r, θ) where {T} + a = metric.spin + Σt = Σ(metric, r, θ) + Δt = Δ(metric, r) + At = A(metric, r, θ) + + return @SMatrix [ + # coords = {t, r, ϕ, θ} + √((Σt * Δt) / At) zero(T) -2*a*r*sin(θ)/√(At * Σt) zero(T) + zero(T) √(Σt / Δt) zero(T) zero(T) + zero(T) zero(T) zero(T) -√Σt + zero(T) zero(T) √(At / Σt)*sin(θ) zero(T) + ] +end + +""" +Jacobian which converts Boyer-Lindquist vector to a ZAMO basis +""" +function jac_zamo_u_bl_d(metric::Kerr{T}, r, θ) where {T} + a = metric.spin + Σt = Σ(metric, r, θ) + Δt = Δ(metric, r) + At = A(metric, r, θ) + + return @SMatrix [# Eq 3.2 1972ApJ...178..347B + # coords = {t, r, ϕ, θ} + √((Σt * Δt) / At) zero(T) zero(T) zero(T) + zero(T) √(Σt / Δt) zero(T) zero(T) + -(2a * r * sin(θ))/√(At * Σt) zero(T) zero(T) √(At / Σt)*sin(θ) + zero(T) zero(T) -√Σt zero(T) + ] +end + +""" +Jacobian which converts ZAMO vector to a Boyer-Lindquist basis +""" +function jac_bl_u_zamo_d(metric::Kerr{T}, r, θ) where {T} + a = metric.spin + Σt = Σ(metric, r, θ) + Δt = Δ(metric, r) + At = A(metric, r, θ) + + return @SMatrix [ + # coords = {t, r, ϕ, θ} + √(At / (Σt * Δt)) zero(T) zero(T) zero(T) + zero(T) √(Δt / Σt) zero(T) zero(T) + zero(T) zero(T) zero(T) -inv(√Σt) + 2a*r/√(Σt * Δt * At) zero(T) √(Σt / At)*csc(θ) zero(T) + ] +end + +""" +Jacobian which expreases ZAMO vector in the fluid frame +""" +function jac_fluid_u_zamo_d(::Kerr{T}, β, θ, φ) where {T} + γ = inv(√(1 - β^2)) + sinφ = sin(φ) + cosφ = cos(φ) + sinθ = sin(θ) + cosθ = cos(θ) + + return @SMatrix [ + γ -β*γ*cosφ*sinθ -β*γ*sinφ*sinθ -β*γ*cosθ + -β*γ*cosφ*sinθ cosθ^2*cosφ^2+γ*cosφ^2*sinθ^2+sinφ^2 (γ-1)*cosφ*sinθ^2*sinφ (γ-1)*cosθ*cosφ*sinθ + -β*γ*sinθ*sinφ (γ-1)*cosφ*sinθ^2*sinφ cosφ^2+(cosθ^2+γ*sinθ^2)*sinφ^2 (γ-1)*cosθ*sinθ*sinφ + -β*γ*cosθ (γ-1)*cosθ*cosφ*sinθ (γ-1)*cosθ*sinθ*sinφ γ*cosθ^2+sinθ^2 + ] +end + + """ Energy reduced azimuthal angular momentum @@ -184,7 +293,7 @@ Defines a horizontal boundary on the assymptotic observer's screen that emission - `metric`: Kerr metric - `θs` : Emission Inclination """ -@inline function αboundary(metric::Kerr, θs) +function αboundary(metric::Kerr, θs) return metric.spin * sin(θs) end @@ -198,7 +307,7 @@ Defines a vertical boundary on the assymptotic observer's screen that emission t - `θo` : Observer inclination - `θs` : Emission Inclination """ -@inline function βboundary(metric::Kerr{T}, α, θo, θs) where {T} +function βboundary(metric::Kerr{T}, α, θo, θs) where {T} a = metric.spin cosθs2 = cos(θs)^2 return √max( @@ -251,7 +360,7 @@ Returns roots of \$r^4 + (a^2-η-λ^2)r^2 + 2(η+(a-λ)^2)r - a^2η\$ - `η` : Reduced Carter constant - `λ` : Reduced azimuthal angular momentum """ -@inline function get_radial_roots(metric::Kerr{T}, η, λ) where {T} +function get_radial_roots(metric::Kerr{T}, η, λ) where {T} a = metric.spin a2 = a * a @@ -293,7 +402,7 @@ Returns roots of \$r^4 + (a^2-η-λ^2)r^2 + 2(η+(a-λ)^2)r - a^2η\$ return roots end -@inline 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 @@ -321,10 +430,10 @@ See [`r_potential(x)`](@ref) for an implementation of \$\\mathcal{R}(r)\$. """ function Ir_inf(metric::Kerr{T}, roots) where {T} #root_diffs = _get_root_diffs(roots...) - numreals = sum(map(_isreal2, roots)) + numreals = sum(_isreal2.(roots)) if numreals == 4 #case 2 - return Ir_inf_case1_and_2(metric, map(real, roots)) + return Ir_inf_case1_and_2(metric, real.(roots)) elseif numreals == 2 #case3 return Ir_inf_case3(metric, roots) else #case 4 @@ -394,6 +503,7 @@ See [`r_potential(x)`](@ref) for an implementation of \$\\mathcal{R}(r)\$. - `νr` : Radial emission direction (Only necessary for case 1&2 geodesics) """ function Ir_s(metric::Kerr{T}, rs, roots, νr) where {T} + isinf(rs) && return Ir_inf(metric, roots) numreals = sum(_isreal2.(roots)) if numreals == 4 #case 2 @@ -748,6 +858,11 @@ function Iϕ_w_I0_terms_case4(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) whe return -2a / (rp - rm) * ((rp - a * λ / 2) * Ip - (rm - a * λ / 2) * Im) end +function Iϕ_o_m_I0_terms(metric::Kerr{T}, rs, roots::NTuple{4}, λ, νr) where {T} + isinf(rs) && return Iϕ_inf(metric, roots, λ) + return Iϕ_w_I0_terms(metric, rs, 0, roots, νr, λ) +end + """ Returns the antiderivative \$I_t=\\int\\frac{r^2\\Delta+2Mr(r^2+a^2-a\\lambda)}{\\sqrt{\\Delta\\mathcal{R(r)}}}dr\$ evaluated at infinity without I0 terms. @@ -954,7 +1069,7 @@ See [`r_potential(x)`](@ref) for an implementation of \$\\mathcal{R}(r)\$. - `λ` : Reduced azimuthal angular momentum - `νr` : Radial emission direction (Only necessary for case 1&2 geodesics) """ -@inline function It_w_I0_terms(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ, νr) where {T} +function It_w_I0_terms(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ, νr) where {T} numreals = sum(_isreal2.(roots)) if numreals == 4 #case 2 @@ -1010,10 +1125,10 @@ end ) #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 #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 coef_p = 2 / √(r31 * r42) * r43 / (rp3 * rp4) coef_m = 2 / √(r31 * r42) * r43 / (rm3 * rm4) @@ -1037,7 +1152,7 @@ end ) end -@inline function It_w_I0_terms_case3(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) where {T} +function It_w_I0_terms_case3(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) where {T} r1, r2, _, _ = roots r21, r31, r32, r41, r42, _ = _get_root_diffs(roots...) r21 = real(r21) @@ -1097,7 +1212,7 @@ end ) end -@inline function It_w_I0_terms_case4(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) where {T} +function It_w_I0_terms_case4(metric::Kerr{T}, rs, τ, roots::NTuple{4}, λ) where {T} a = metric.spin r1, _, _, r4 = roots _, r31, r32, r41, r42, _ = _get_root_diffs(roots...) @@ -1153,7 +1268,17 @@ end )# + (logdiv + lineardiv) end -@inline function radial_inf_integrals(met::Kerr{T}, roots::NTuple{4}) where {T} +function It_o_m_I0_terms(metric::Kerr{T}, ro, roots::NTuple{4}, λ, νr) where {T} + isinf(ro) && return It_inf(metric, roots, λ) + return It_w_I0_terms(metric, ro, 0, roots, λ, νr) +end + +function radial_o_m_I0_terms_integrals(met::Kerr{T}, ro::T, roots::NTuple{4}, νr) where {T} + isinf(ro) && return radial_inf_integrals(met, roots) + return radial_w_I0_terms_integrals(met, ro, roots, zero(T), νr) +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) @@ -1168,7 +1293,7 @@ end """ Returns the radial integrals for the case where there are four real roots in the radial potential, with roots outside the horizon. """ -@inline function radial_inf_integrals_case2(metric::Kerr{T}, roots::NTuple{4}) where {T} +function radial_inf_integrals_case2(metric::Kerr{T}, roots::NTuple{4}) where {T} roots = real.(roots) _, _, r3, r4 = roots _, r31, r32, r41, r42, _ = _get_root_diffs(roots...) @@ -1217,7 +1342,7 @@ end """ Returns the radial integrals for the case where there are two real roots in the radial potential """ -@inline function radial_inf_integrals_case3(metric::Kerr{T}, roots::NTuple{4}) where {T} +function radial_inf_integrals_case3(metric::Kerr{T}, roots::NTuple{4}) where {T} r1, r2, _, _ = roots r21, r31, r32, r41, r42, _ = _get_root_diffs(roots...) @@ -1260,7 +1385,7 @@ end """ Returns the radial integrals for the case where there are no real roots in the radial potential """ -@inline function radial_inf_integrals_case4(metric::Kerr{T}, roots::NTuple{4}) where {T} +function radial_inf_integrals_case4(metric::Kerr{T}, roots::NTuple{4}) where {T} r1, _, _, r4 = roots _, r31, r32, r41, r42, _ = _get_root_diffs(roots...) a2 = metric.spin^2 @@ -1335,7 +1460,7 @@ end """ Returns the radial integrals for the case where there are four real roots in the radial potential, with roots outside the horizon. """ -@inline function radial_w_I0_terms_integrals_case2( +function radial_w_I0_terms_integrals_case2( metric::Kerr{T}, rs, roots::NTuple{4}, @@ -1378,10 +1503,10 @@ Returns the radial integrals for the case where there are four real roots in the ) #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 #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 + I2_total = (r1 * r4 + r2 * r3) / 2 * τ - (-1)^νr * I2_s coef_p = 2 / √(r31 * r42) * r43 / (rp3 * rp4) coef_m = 2 / √(r31 * r42) * r43 / (rm3 * rm4) @@ -1403,7 +1528,7 @@ end """ Returns the radial integrals for the case where there are two real roots in the radial potential """ -@inline function radial_w_I0_terms_integrals_case3( +function radial_w_I0_terms_integrals_case3( metric::Kerr{T}, rs, roots::NTuple{4}, @@ -1441,9 +1566,7 @@ Returns the radial integrals for the case where there are two real roots in the Π2_s = (coef^2) * R2(αo, φ_s, k3) I0_total = τ - # Removed logarithmic divergence 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) @@ -1461,7 +1584,7 @@ end """ Returns the radial integrals for the case where there are no real roots in the radial potential """ -@inline function radial_w_I0_terms_integrals_case4( +function radial_w_I0_terms_integrals_case4( metric::Kerr{T}, rs, roots::NTuple{4}, @@ -1491,14 +1614,11 @@ Returns the radial integrals for the case where there are no real roots in the r (isnan(x4_s) || isnan(x4_p) || isnan(x4_m)) && return zero(T), zero(T), zero(T), zero(T) S1p_s = S1(gp, atan(x4_s) + atan(go), k4) S1m_s = S1(gm, atan(x4_s) + atan(go), k4) - #Fr_s = 2 / (C + D) * JacobiElliptic.F(atan(x4_s) + atan(go), k4) Π1_s = 2 / (C + D) * (a2 / go * (1 + go^2)) * S1(go, atan(x4_s) + atan(go), k4) Π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) - # 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)) * @@ -1532,6 +1652,22 @@ function total_mino_time(metric::Kerr{T}, roots::NTuple{4}) where {T} return τf end +function total_mino_time(metric::Kerr{T}, ro, roots::NTuple{4}) where {T} + isinf(ro) && return total_mino_time(metric, roots) + + numreals = unsafe_trunc(Int, sum((Krang._isreal2.(roots)))) + Io = Krang.Ir_s(metric, ro, roots, true) + + if numreals == 4 + τf = Io + Ir_inf(metric, roots) + else + rh = Krang.horizon(metric) + τf = Io - Krang.Ir_s(metric, rh, roots, true) + end + return τf +end + + ##---------------------------------------------------------------------------------------------------------------------- # Pixel utility functions ##---------------------------------------------------------------------------------------------------------------------- @@ -1547,7 +1683,7 @@ See [`r_potential(x)`](@ref) for an implementation of \$\\mathcal{R}(r)\$. - `rs` : Emission radius """ function Ir(pix::AbstractPixel, νr::Bool, rs) - return I0_inf(pix) - Ir_s(metric(pix), rs, roots(pix), νr) + return I0_o(pix) - Ir_s(metric(pix), rs, roots(pix), νr) end """ @@ -1563,7 +1699,7 @@ See [`r_potential(x)`](@ref) for an implementation of \$\\mathcal{R}(r)\$. """ function Iϕ(pix::AbstractPixel, rs, τ, νr) metric = pix.metric - λtemp = λ(metric, pix.screen_coordinate[1], pix.θo) + λtemp = λ(pix) tempIϕ_inf = Iϕ_inf(pix) tempIϕ_s = Iϕ_w_I0_terms(metric, rs, τ, pix.roots, νr, λtemp) @@ -1583,7 +1719,7 @@ See [`r_potential(x)`](@ref) for an implementation of \$\\mathcal{R}(r)\$. """ function It(pix::AbstractPixel, rs, τ, νr) metric = pix.metric - λtemp = λ(metric, pix.screen_coordinate[1], pix.θo) + λtemp = λ(pix) tempIt_inf = It_inf(pix) return tempIt_inf - It_w_I0_terms(metric, rs, τ, pix.roots, λtemp, νr) @@ -1597,14 +1733,23 @@ Return the radial integrals - `τ` : Mino time - `νr` : Sign of radial velocity direction at emission. This is always positive for case 3 and case 4 geodesics. """ -@inline function radial_integrals(pix::AbstractPixel, rs, τ, νr) +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_o, I2_o, Ip_o, Im_o = radial_o_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 end -@inline function _rs_case1_and_2(pix::AbstractPixel, rh, τ::T)::Tuple{T,Bool,Bool} where {T} +function radial_o_integrals_m_I0_terms(pix::AbstractPixel) + if isinf(pix.ro) + return radial_inf_integrals(metric(pix), pix.roots) + else + return radial_w_I0_terms_integrals(metric(pix), pix.ro, pix.roots, 0, true) + end +end + + +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...) @@ -1617,26 +1762,26 @@ end coef = 2 / √real(r31 * r42) Ir_s = !(x2_s < one(T)) ? zero(T) : coef * JacobiElliptic.F(asin(x2_s), k) - fo = I0_inf(pix) + fo = I0_o(pix) r4 < rh && τ > (fo - Ir_s) && return zero(T), false, false# invalid case2 X2 = √(r31 * r42) * (fo - τ) / 2 - if τ > 2fo + if τ > (fo + I0_inf(pix)) return zero(T), false, false end sn = r41 * JacobiElliptic.sn(X2, k)^2 return (r31 * r4 - r3 * sn) / (r31 - sn), X2 > zero(T), true end -@inline 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...) r1, r2, r21 = real.((r1, r2, r21)) - fo = I0_inf(pix) + fo = I0_o(pix) A = √abs(r32 * r42) B = √abs(r31 * r41) k = (((A + B)^2 - r21^2) / (4 * A * B)) @@ -1657,7 +1802,7 @@ end return real(num / den), X3 > zero(T), true end -@inline 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...) @@ -1676,7 +1821,7 @@ end coef = 2 / (C + D) Ir_s = coef * JacobiElliptic.F(atan(x4_s) + atan(go), k4) - fo = I0_inf(pix) + fo = I0_o(pix) τ > (fo - Ir_s) && return zero(T), false, false X4 = (C + D) / T(2) * (fo - τ) @@ -1786,11 +1931,11 @@ See [`θ_potential(x)`](@ref) for an implementation of \$\\Theta(\theta)\$. end function Gs(pix::AbstractPixel, τ::T) where {T} - α, β = screen_coordinate(pix) + _, β = screen_coordinate(pix) θo = inclination(pix) met = metric(pix) - ηtemp = η(met, α, β, θo) - λtemp = λ(met, α, θo) + ηtemp = η(pix) + λtemp = λ(pix) signβ = sign(β) τ == T(NaN) && return T(NaN) @@ -1837,14 +1982,14 @@ function Gs(pix::AbstractPixel, τ::T) where {T} return Δτ end -@inline function Gϕ(pix::AbstractPixel, θs::T, isindir, n) where {T} - α, β = screen_coordinate(pix) +function Gϕ(pix::AbstractPixel, θs::T, isindir, n) where {T} + _, β = screen_coordinate(pix) met = metric(pix) θo = inclination(pix) signβ = sign(β) - ηtemp = η(met, α, β, θo) - λtemp = λ(met, α, θo) + ηtemp = η(pix) + λtemp = λ(pix) a = met.spin Go, Ghat = absGϕo_Gϕhat(pix) @@ -1900,14 +2045,14 @@ end return ans, Gs, Go, Ghat, isvortical, true end -@inline function Gt(pix::AbstractPixel, θs::T, isindir, n) where {T} - α, β = screen_coordinate(pix) +function Gt(pix::AbstractPixel, θs::T, isindir, n) where {T} + _, β = screen_coordinate(pix) met = metric(pix) θo = inclination(pix) signβ = sign(β) - ηtemp = η(met, α, β, θo) - λtemp = λ(met, α, θo) + ηtemp = η(pix) + λtemp = λ(pix) a = met.spin Go, Ghat = absGto_Gthat(pix) Gs, ans, isvortical = zero(T), zero(T), ηtemp < zero(T) @@ -1969,7 +2114,7 @@ end # SlowLightIntensityCachedPixel utility functions ##---------------------------------------------------------------------------------------------------------------------- -@inline function _absGθo_Gθhat(metric::Kerr{T}, θo, η, λ)::NTuple{2,T} where {T} +function _absGθo_Gθhat(metric::Kerr{T}, θo, η, λ)::NTuple{2,T} where {T} a = metric.spin a2 = a^2 Go, Ghat, isvortical = zero(T), zero(T), η < zero(T) @@ -2009,7 +2154,7 @@ end return Go, Ghat end -@inline function _absGϕo_Gϕhat(metric::Kerr{T}, θo, η, λ)::NTuple{2,T} where {T} +function _absGϕo_Gϕhat(metric::Kerr{T}, θo, η, λ)::NTuple{2,T} where {T} a = metric.spin Go, Ghat, isvortical = zero(T), zero(T), η < zero(T) @@ -2049,7 +2194,7 @@ end return Go, Ghat end -@inline function _absGto_Gthat(metric::Kerr{T}, θo, η, λ)::NTuple{2,T} where {T} +function _absGto_Gthat(metric::Kerr{T}, θo, η, λ)::NTuple{2,T} where {T} a = metric.spin Go, Ghat, isvortical = zero(T), zero(T), η < zero(T) diff --git a/test/emission_coordinates_tests.jl b/test/emission_coordinates_tests.jl index 9196de7a..e3a8f86f 100644 --- a/test/emission_coordinates_tests.jl +++ b/test/emission_coordinates_tests.jl @@ -256,7 +256,7 @@ end end - @test -(Iϕ + λtemp * Gϕ + π / 2) / ϕs ≈ 1.0 atol = 1e-3 + @test (Iϕ + λtemp * Gϕ) / ϕs ≈ 1.0 atol = 1e-3 ft(r, p) = ( (r^2 * (r^2 - 2r + a^2) + 2r * (r^2 + a^2 - a * λtemp)) * diff --git a/test/finite_observer_tests.jl b/test/finite_observer_tests.jl new file mode 100644 index 00000000..bb8eb92e --- /dev/null +++ b/test/finite_observer_tests.jl @@ -0,0 +1,106 @@ +@testset "Finite Observer" begin + metric = Krang.Kerr(0.01); # Kerr metric with a spin of 0.99 + θo = 45 * π / 180; # inclination angle of the observer. θo ∈ (0, π) + θs = π/2 + ro = 1e8 + sze = 40; # resolution of the screen is sze x sze + rmin = Krang.horizon(metric); # minimum radius to be ray traced + rmax = 6.0; # maximum radius to be ray traced + ρmax = 10; # horizontal and vertical limits of the screen + + infinite_coordinates = [zeros(sze, sze) for _ = 1:3] + finite_coordinates = [zeros(sze, sze) for _ = 1:3] + infinite_camera = + Krang.SlowLightIntensityCamera(metric, θo, -ρmax, ρmax, -ρmax, ρmax, sze); + finite_camera = Krang.IntensityCamera( + metric, + θo, + ro, + -ρmax/ro, + ρmax/ro, + -ρmax/ro, + ρmax/ro, + sze, + ); + + @test all( + [i.screen_coordinate[1] .* ro for i in finite_camera.screen.pixels] .≈ [i.screen_coordinate[1] for i in infinite_camera.screen.pixels], + ) + + @test maximum([ + abs.((finite_camera.screen.pixels[i].η / infinite_camera.screen.pixels[i].η))-1 for + i in range(1, length(finite_camera.screen.pixels)) + ]) ≈ 0.0 atol = 1e-5 + @test maximum([ + abs.((finite_camera.screen.pixels[i].λ / infinite_camera.screen.pixels[i].λ))-1 for + i in range(1, length(finite_camera.screen.pixels)) + ]) ≈ 0.0 atol = 1e-5 + + 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 _ = 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) + if rs ≤ rmin || rs ≥ rmax + continue + end + coords = isnan(rs) ? observation : (ts, rs, θs, ϕs) + end + return coords + end + + + function infinite_draw!(camera, coordinates, rmin, rmax, θs, n) + times, radii, azimuths = coordinates + + geometry = Krang.ConeGeometry(θs, (n, rmin, rmax)) + + rendered_scene = coordinate_point.(camera.screen.pixels, Ref(geometry)) + for I in CartesianIndices(rendered_scene) + temp = rendered_scene[I][1] + times[I] = temp + radii[I] = rendered_scene[I][2] + azimuths[I] = mod2pi(rendered_scene[I][4]) # azimuths are in radians + end + coordinates[1] .= times + coordinates[2] .= radii + coordinates[3] .= azimuths + + end + + function finite_draw!(camera, coordinates, rmin, rmax, θs, n) + times, radii, azimuths = coordinates + + geometry = Krang.ConeGeometry(θs, (n, rmin, rmax)) + rendered_scene = coordinate_point.(camera.screen.pixels, Ref(geometry)) + for I in CartesianIndices(rendered_scene) + temp = rendered_scene[I][1] + times[I] = temp == 0 ? 0 : temp - 2log(ro) - ro + radii[I] = rendered_scene[I][2] + azimuths[I] = mod2pi(rendered_scene[I][4]) # azimuths are in radians + end + coordinates[1] .= times + coordinates[2] .= radii + coordinates[3] .= azimuths + + end + + @testset "n=$n" for n = 0:2 + finite_draw!(finite_camera, finite_coordinates, rmin, rmax, θs, n) + infinite_draw!(infinite_camera, infinite_coordinates, rmin, rmax, θs, n) + + @testset "$coord" for (i, coord) in enumerate(["radius", "azimuth"]) + residual = maximum(abs.(finite_coordinates[i+1] .- infinite_coordinates[i+1])) + @test residual ≈ 0.0 atol = 1e-4 + end + end + +end diff --git a/test/runtests.jl b/test/runtests.jl index 747094e7..6198cf5a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,7 @@ using LinearAlgebra using Integrals using StaticArrays using FileIO +import Downloads using Rotations import Downloads using GeometryBasics @@ -15,3 +16,4 @@ include("emission_coordinates_tests.jl") include("camera_tests.jl") include("raytracer_tests.jl") include("mesh_geometry_tests.jl") +include("finite_observer_tests.jl")