Add 3d ray intersect by AreWeDreaming · Pull Request #269 · ProjectTorreyPines/IMAS.jl

AreWeDreaming · 2025-06-17T18:37:39Z

This adds to geometry routines that compute intersections between a ray and a torus with a rectangular cross section.
Still work in progress because I need to remove ray_rect_intersection since that is untested LLM code that I don't need.

AreWeDreaming · 2025-06-17T18:39:13Z

Should be ready to go. This has extensive tests that should cover almost every line.

bclyons12

I made a bunch of suggestions to Julia-fy the code, but does the existing toroidal_intersection() do what you want?

bclyons12 · 2025-06-17T21:45:21Z

src/geometry.jl

+            t_interval_z = [0.0 Inf]
+        end
+    else
+        t_interval_z = sort![(z_min - z0) / dz, (z_max - z0) / dz]


sort![...] can't be right?

Suggested change

t_interval_z = sort![(z_min - z0) / dz, (z_max - z0) / dz]

t_interval_z = sort([(z_min - z0) / dz, (z_max - z0) / dz])

sort!([(z_min - z0) / dz, (z_max - z0) / dz]) at the every least, but sort(@SVector[(z_min - z0) / dz, (z_max - z0) / dz]) would be better.

bclyons12 · 2025-06-17T21:47:16Z

src/geometry.jl

+    crossing = zeros(Float64, 4) # +1 -> into domain
+                                 # -1 -> out of 
+                                 # 0 -> no crossing
+    into_domain = [1.0, -1.0]


All these small, fixed-sized arrays should either be Tuples or SVector/MVector from StaticArrays.jl. This will avoid unnecessary allocations.

Suggested change

into_domain = [1.0, -1.0]

into_domain = Tuple([1.0, 1.0])

into_domain = (1.0, 1.0)

bclyons12 · 2025-06-17T21:49:05Z

src/geometry.jl

+        if (r0 > r_min && r0 < r_max)
+            return [[0.0 Inf];]
+        else
+            return zeros(Float64, 0)


You want an empty array? Float64[] is more idiomatic.

Suggested change

return zeros(Float64, 0)

return Float64[]

bclyons12 · 2025-06-17T21:53:39Z

src/geometry.jl

+    if dx == 0 && dy == 0
+        # Handle vertical ray
+        if (r0 > r_min && r0 < r_max)
+            return [[0.0 Inf];]


[0.0 Inf] will give you a 1x2 Matrix

Suggested change

return [[0.0 Inf];]

return [0.0 Inf]

bclyons12 · 2025-06-17T21:54:54Z

src/geometry.jl

+        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


If this ever becomes computationally intensive, there are plenty of redundant FLOPs here that could be avoided.

This is run once per beam so it should not come up too often and cost negligible time compared to the actual raytracing.

bclyons12 · 2025-06-17T21:55:25Z

src/geometry.jl

+        append!(crossing, 1.0)
+    end
+    i_sort = sortperm(t_crossings)
+    intervals = zeros(Float64, 0)


Float64[]

bclyons12 · 2025-06-17T21:57:00Z

src/geometry.jl

+        end
+    end
+    if length(intervals) % 2 != 0
+        throw(ErrorException("Somehow only found odd number of intersections which is impossible."))


Make sure grazing intersections are properly handled.

For my specific use-case grazing intersection can raise an error. A beam that hits the rectangular flux matrix parallel to one of its boundaries indicates an error in how the beam is set up.

AreWeDreaming · 2025-06-17T22:50:43Z

src/geometry.jl

+
+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)


Suggested change

t_crossings = zeros(Float64, 4)

t_crossings = Float64[]

This changed the size.

AreWeDreaming · 2025-06-17T22:52:44Z

src/geometry.jl

+        append!(crossing, 1.0)
+    end
+    i_sort = sortperm(t_crossings)
+    intervals = zeros(Float64, 0)


Suggested change

intervals = zeros(Float64, 0)

intervals = Float64[]

AreWeDreaming · 2025-06-17T22:57:49Z

As far as I can tell toroidal_intersection in particles.jl is iterative. This is analytic and therefore should be more accurate and computationally more efficient.

bclyons12 · 2025-06-17T23:05:09Z

toroidal_intersection() is analytic. It iterates if you give it multiple segments to try to intersect with. It was originally made for a ballistic trajectories of particles hitting a wall, but it's should work here.

IMAS.jl/src/physics/particles.jl

Lines 296 to 379 in ad841d4

    
           """ 
        
               toroidal_intersection(r1::Real, z1::Real, r2::Real, z2::Real, px::Real, py::Real, pz::Real, vx::Real, vy::Real, vz::Real, v2::Real, vz2::Real) 
        
           Compute the time of intersection between a moving particle and a toroidal surface defined by two points `(r1, z1)` and `(r2, z2)`. 
        
           Returns the smallest positive time at which the particle intersects the surface segment. Returns `NaN` if no valid intersection occurs. 
        
             - `r1`, `z1`: Coordinates of the first endpoint of the toroidal surface segment. 
        
             - `r2`, `z2`: Coordinates of the second endpoint of the toroidal surface segment. 
        
             - `px`, `py`, `pz`: Current position of the particle in Cartesian coordinates. 
        
             - `vx`, `vy`, `vz`: Velocity components of the particle. 
        
             - `v2`: Squared radial velocity component, `vx^2 + vy^2`. 
        
             - `vz2`: Squared z-velocity component, `vz^2`. 
        
           """ 
        
           function toroidal_intersection(r1::Real, z1::Real, r2::Real, z2::Real, px::Real, py::Real, pz::Real, vx::Real, vy::Real, vz::Real, v2::Real, vz2::Real) 
        
               t = Inf 
        
               if isapprox(z1, z2) 
        
                   # a horizontal segment: time for how long it takes to get to this z, then check if r is between segment bounds 
        
                   ti = (vz == 0.0) ? -Inf : (z1 - pz) / vz 
        
                   if ti >= 0 
        
                       ri = sqrt((px + vx * ti)^2 + (py + vy * ti)^2) 
        
                       if ri >= r1 && ri <= r2 
        
                           t = ti 
        
                       end 
        
                   end 
        
               else 
        
                   # parameterize parabola intersection with a line 
        
                   m = (z2 - z1) / (r2 - r1) 
        
                   z0 = z1 - m * r1 
        
                   m2 = m^2 
        
                   a = 0.0 
        
                   b = 0.0 
        
                   c = 0.0 
        
                   if isfinite(m) 
        
                       z_z0 = pz - z0 
        
                       a = m2 * v2 - vz2 
        
                       b = 2 * (m2 * (vx * px + vy * py) - z_z0 * vz) 
        
                       c = m2 * (px^2 + py^2) - z_z0^2 
        
                   else 
        
                       a = v2 
        
                       b = 2 * (vx * px + vy * py) 
        
                       c = px^2 + py^2 - r1^2 
        
                   end 
        
                   t1 = -Inf 
        
                   t2 = -Inf 
        
                   if a == 0 
        
                       b != 0 && (t2 = -c / b) 
        
                   else 
        
                       nb_2a = -b / (2a) 
        
                       b2a_ca = (b^2 - 4 * a * c) / (4 * a^2) 
        
                       if b2a_ca > 0 
        
                           t_sq = sqrt(b2a_ca) 
        
                           t1 = nb_2a - t_sq 
        
                           t2 = nb_2a + t_sq 
        
                       elseif b2a_ca ≈ 0.0 
        
                           t2 = nb_2a 
        
                       end 
        
                   end 
        
                   # check if intersection points are at positive time and fall between segment 
        
                   t = Inf 
        
                   if t1 >= 0 
        
                       zi = vz * t1 + pz 
        
                       if zi >= z1 && zi <= z2 
        
                           t = t1 
        
                       end 
        
                   end 
        
                   if !isfinite(t) && t2 >= 0 
        
                       zi = vz * t2 + pz 
        
                       if zi >= z1 && zi <= z2 
        
                           t = t2 
        
                       end 
        
                   end 
        
               end 
        
               if t == Inf 
        
                   t = NaN 
        
               end 
        
               return t 
        
           end

AreWeDreaming · 2025-06-18T00:06:04Z

@bclyons12 you were right, the existing routines do what I need to do. I am not sure why I didn't want to use them initially. I almost everything. Now all this PR does is add another dispatch to the toroidal_intersection method such that it accepts p and v as vectors. I also rewrote the test I contribute here so that it tests the toroidal_intersection method instead.

AreWeDreaming · 2025-06-18T00:06:53Z

(I also tested the regression test in TorJ which now uses toroidal_intersection and that has passed as well.)

src/physics/particles.jl

Add assertion

AreWeDreaming added 5 commits May 13, 2025 13:32

This adds a routine for ray intersection with a torus

99251b1

Merge remote-tracking branch 'origin/master' into add_3d_ray_intersect

5f45534

Sneak in little addition for new source

2115990

Add tests and fix bugs

4fd08d1

Merge branch 'master' into add_3d_ray_intersect

4d6516f

AreWeDreaming requested review from bclyons12 and orso82 June 17, 2025 18:37

AreWeDreaming and others added 2 commits June 17, 2025 11:38

Remove uneeded, untested routine

da464d5

Merge branch 'master' into add_3d_ray_intersect

756f013

AreWeDreaming marked this pull request as ready for review June 17, 2025 18:38

Make sure geometry tests are run

eb2edce

bclyons12 requested changes Jun 17, 2025

View reviewed changes

AreWeDreaming commented Jun 17, 2025

View reviewed changes

AreWeDreaming added 4 commits June 17, 2025 16:54

Welp theere goes my contribution

544e4ac

Fix accidental deletion when merging

ef97206

Wrong spot

91a5b95

Remove cruft

14f7d46

bclyons12 reviewed Jun 18, 2025

View reviewed changes

src/physics/particles.jl Show resolved Hide resolved

bclyons12 approved these changes Jun 18, 2025

View reviewed changes

AreWeDreaming added 2 commits August 12, 2025 10:59

Update src/physics/particles.jl

381cc89

Add assertion

Merge branch 'master' into add_3d_ray_intersect

6195ab4

	t_interval_z = sort![(z_min - z0) / dz, (z_max - z0) / dz]
	t_interval_z = sort([(z_min - z0) / dz, (z_max - z0) / dz])

Conversation

AreWeDreaming commented Jun 17, 2025

Uh oh!

AreWeDreaming commented Jun 17, 2025

Uh oh!

bclyons12 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

AreWeDreaming commented Jun 17, 2025

Uh oh!

bclyons12 commented Jun 17, 2025

Uh oh!

AreWeDreaming commented Jun 18, 2025

Uh oh!

AreWeDreaming commented Jun 18, 2025

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants