Add 3d ray intersect

src/geometry.jl

@@ -3,6 +3,7 @@ document[Symbol("Geometry")] = Symbol[]
 import StaticArrays
 import MillerExtendedHarmonic: MXH
 import LinearAlgebra
+import Roots
 
 """
     centroid(x::AbstractVector{<:T}, y::AbstractVector{<:T}) where {T<:Real}
@@ -81,6 +82,130 @@ end
 
 @compat public revolution_volume
 push!(document[Symbol("Geometry")], :revolution_volume)
+"""
+    r_intersect_interval(x0::T, y0::T, dx::T, dy::T, r_min::T, r_max::T) where {T<:Real}
+
+Finds intersection of ray starting at x0 and y0 with direction (dx,dy) 
+    with reference major radius value r_ref
+
+- `x0``, `y0``: Origin
+- `dx``, `dy``: Normalized direction in x-y plane
+- `r_min`, `r_max`: Edges of r_range
+
+Returns interval [t1, t2] for which r_min <= r(t) <= r_max with r(t)^2 = (x+dx*t)^2 + (y+dy*t)^2
+
+"""
+
+function solve_r_intersect(x0::T, y0::T, dx::T, dy::T, r_min::T, r_max::T) where {T<:Real}
+    r0 = sqrt(x0^2 + y0^2)
+    t_crossings = zeros(Float64, 4)
+    crossing = zeros(Float64, 4) # +1 -> into domain
+                                 # -1 -> out of 
+                                 # 0 -> no crossing
+    into_domain = [1.0, -1.0]
+    if dx == 0 && dy == 0
+        # Handle vertical ray
+        if (r0 > r_min && r0 < r_max)
+            return [[0.0 Inf];]
+        else
+            return zeros(Float64, 0)
+        end
+    end
+    for (i_ref, r_ref) in enumerate([r_min, r_max])
+        Δ = r_ref^2 * (dx^2 + dy^2) - x0^2 * dy^2 - y0^2 * dx^2 + 2 * x0 * y0* dx*dy
+        if Δ < 0
+            # println("No intersection for ", r_ref)
+            t_crossings[2*(i_ref-1) + 1] = Inf
+            t_crossings[2*(i_ref-1) + 2] = Inf
+            crossing[2*(i_ref-1) + 1] = 0.0
+            crossing[2*(i_ref-1) + 2] = 0.0
+            continue
+        end
+        t1 = -(x0*dx + y0 * dy) - sqrt(Δ)
+        t2 = -(x0*dx + y0 * dy) + sqrt(Δ)
+        t_crossings[2*(i_ref-1) + 1] = t1 >= 0 ? t1/ (dx^2 + dy^2) : Inf
+        t_crossings[2*(i_ref-1) + 2] = t2 >= 0 ? t2/ (dx^2 + dy^2) : Inf
+        crossing[2*(i_ref-1) + 1] = t1 >= 0 ? into_domain[i_ref]*sign(x0*dx + dx^2*t1 + y0 * dy + dy^2*t1) : 0.0
+        crossing[2*(i_ref-1) + 2] = t2 >= 0 ? into_domain[i_ref]*sign(x0*dx + dx^2*t2 + y0 * dy + dy^2*t2) : 0.0
+    end
+    last = 0 # 1 in, 0 out
+    # If the first point is in we add it to the crossings
+    if (r0 > r_min && r0 < r_max)
+        append!(t_crossings, 0.0)
+        append!(crossing, 1.0)
+    end
+    i_sort = sortperm(t_crossings)
+    intervals = zeros(Float64, 0)
+    for i in i_sort
+        if t_crossings[i] == Inf
+            continue
+        end
+        if last == 0 && crossing[i] > 0.0 # looking for a crossing into the rectangle
+            append!(intervals, t_crossings[i])
+            last = 1
+        elseif last == 1 && crossing[i] < 0.0 # looking for a crossing leaving the rectangle
+            append!(intervals, t_crossings[i])
+            last = 0
+        end
+    end
+    if length(intervals) % 2 != 0
+        throw(ErrorException("Somehow only found odd number of intersections which is impossible."))
+    end
+    if length(intervals) > 2  # Return as shape (2,2). Need to transpose because column-major...
+        return permutedims(reshape(intervals, :, 2))
+    elseif length(intervals) > 0
+        return reshape(intervals, :, 2)  # Return as shape (1,2)
+    else
+        return intervals # Return empty
+    end
+end
+
+
+"""
+    ray_torus_intersect(origin, direction, ρ_bounds, z_bounds)
+
+Finds intersection of a ray with a rectangular torus defined in cylindrical coordinates.
+
+- `origin`: (x, y, z) vector
+- `direction`: normalized (dx, dy, dz)
+- `r_bounds`: (r_min, r_max)
+- `z_bounds`: (z_min, z_max)
+
+Returns intersection point between ray closest to the origin or `nothing` if no intersection.
+"""
+function ray_torus_intersect(origin, direction, r_bounds, z_bounds)
+    x0, y0, z0 = origin
+    if norm(direction) == 0
+        throw(ArgumentError("Direction most not have norm zero"))
+    end
+    dx, dy, dz = direction/norm(direction)
+    r_min, r_max = r_bounds
+    z_min, z_max = z_bounds
+    # Intervals in between the ray passes throught the r limits of the box
+    t_intervals_R = solve_r_intersect(x0, y0, dx, dy, r_min, r_max)
+    if dz == 0
+        if z0 < z_min || z0 > z_max
+            return nothing, nothing
+        else
+            t_interval_z = [0.0 Inf]
+        end
+    else
+        t_interval_z = sort![(z_min - z0) / dz, (z_max - z0) / dz]
+    end
+    t_interval_z[1] = t_interval_z[1] > 0.0 ? t_interval_z[1] : 0.0 # Remove intersection at negative values
+
+
+    # Parametrize ρ(t) and z(t)
+    for t_interval_R in eachrow(t_intervals_R)
+        t_min = max(t_interval_R[1], t_interval_z[1])
+        t_max = min(t_interval_R[2], t_interval_z[2])
+        if t_min < t_max
+            return [x0 y0 z0] + [dx dy dz] * t_min
+        end
+    end
+    return nothing
+end
+
 
 """
     intersection_angles(

src/physics/sources.jl

@@ -656,7 +656,6 @@ function new_source(
         setproperty!(electrons, :energy, value; error_on_missing_coordinates = false)
         setproperty!(electrons, :power_inside, cumtrapz(volume, value); error_on_missing_coordinates = false)
     elseif electrons_power_inside !== missing
-        value = electrons_power_inside
         setproperty!(electrons, :power_inside, value; error_on_missing_coordinates = false)
         setproperty!(electrons, :energy, gradient(volume, value); error_on_missing_coordinates = false)
     else

test/geometry/test_ray_torus_intersect.jl

@@ -0,0 +1,16 @@
+using IMAS
+using Test
+
+@testset "Test ray_torus_intersect" begin
+    test_cases = [
+        ([6.0,0.0,0.0], [-1.0,0.0,0.0], [1.0, 4.0], [-1.0, 1.0], [4.0 0.0 0.0], )
+    ]
+    for (origin, direction, r_bounds, z_bounds, expected) in test_cases
+        @testset "origin=$origin, direction=$direction, r_bounds=$r_bounds, z_bounds=$z_bounds" begin
+            # println("--- x0=$x0, y0=$y0, dx=$dx, dy=$dy, r_min=$r_min, r_max=$r_max ---")
+            test_result = IMAS.ray_torus_intersect(origin, direction, r_bounds, z_bounds)
+            println("--- Result | expected: $test_result | $expected  ---")
+            @test test_result ≈ expected
+        end
+    end
+end

test/geometry/test_solve_r_intersect.jl

@@ -0,0 +1,24 @@
+using IMAS
+using Test
+
+@testset "Test solve_r_intersect" begin
+    test_cases = [
+        (3.0, 0.0, 0.0, 0.0, 1.0, 4.0, [[0.0 Inf];]), # ray is not moving, but starts on the grid
+        (6.0, 0.0, 0.0, 0.0, 1.0, 4.0, zeros(Float64, 0)), # ray is not moving, and starts off the grid
+        (6.0, 0.0, -1.0, 0.0, 1.0, 4.0, [2. 5.; 7. 10.]), # ray in x-direction 
+        (0.0, 6.0, 0.0, -1.0, 2.0, 4.0, [2. 4.; 8. 10.]), # ray in y-direction
+        (3.0, 0.0, -1.0, 0.0, 1.0, 4.0, [0.0 2.; 4.0 7.0]), # ray in x-direction with start in torus
+        (6.0, 6.0, 0.0, -1.0, 2.0, 4.0, zeros(Float64, 0)), # ray in x-direction with offset causing no intersection
+        (3.0, 6.0, -3.0/5.0, -4.0/5.0, 1.0, 2.0, [[5.0 41.0/5.0];]), # Single crossing with x-y ray
+        (3.0, 8.0, -3.0/5.0, -4.0/5.0, 3.0, 4.0, [5.0 32.0/5.0; 10.0 57.0/5.0]) # double crossing with x-y ray
+    ]
+
+    for (x0, y0, dx, dy, r_min, r_max, expected) in test_cases
+        @testset "x0=$x0, y0=$y0, dx=$dx, dy=$dy, r_min=$r_min, r_max=$r_max" begin
+            # println("--- x0=$x0, y0=$y0, dx=$dx, dy=$dy, r_min=$r_min, r_max=$r_max ---")
+            test_result = IMAS.solve_r_intersect(x0, y0, dx, dy, r_min, r_max)
+            println("--- Result | expected: $test_result | $expected  ---")
+            @test test_result ≈ expected
+        end
+    end
+end

test/runtests.jl

@@ -11,4 +11,8 @@ else
     include("runtests_extract.jl")
 
     include("runtests_plot_recipes.jl")
+
+    include("geometry/test_ray_torus_intersect.jl")
+
+    include("geometry/test_solve_r_intersect.jl")
 end