Fix wall_hit_cart/xend_cart inconsistency when from_cart fails

CMakeLists.txt

@@ -209,16 +209,16 @@ if(SIMPLE_ENABLE_CGAL)
     include(FetchContent)
     FetchContent_Declare(
         cgal
-        URL https://github.com/CGAL/cgal/releases/download/v5.6.1/CGAL-5.6.1.tar.xz
+        URL https://github.com/CGAL/cgal/releases/download/v6.0.1/CGAL-6.0.1.tar.xz
         DOWNLOAD_EXTRACT_TIMESTAMP TRUE
     )
     FetchContent_GetProperties(cgal)
     if(NOT cgal_POPULATED)
         FetchContent_Populate(cgal)
     endif()
-    # CGAL 5.x is header-only, just add include path
+    # CGAL 6.x is header-only, just add include path
     set(CGAL_INCLUDE_DIR "${cgal_SOURCE_DIR}/include")
-    message(STATUS "CGAL enabled via FetchContent: 5.6.1 (header-only)")
+    message(STATUS "CGAL enabled via FetchContent: 6.0.1 (header-only)")
     message(STATUS "CGAL include dir: ${CGAL_INCLUDE_DIR}")
     add_compile_definitions(SIMPLE_ENABLE_CGAL)
 else()

src/netcdf_results_output.f90

@@ -366,13 +366,14 @@ end subroutine write_results_netcdf
     subroutine compute_cartesian_positions(xstart_cart, xend_cart)
         !> Convert reference coordinates to Cartesian for all particles.
         !>
-        !> For particles that were not traced (zend = 0), xend_cart is set
-        !> to xstart_cart since they effectively remained at the start.
+        !> For particles that were not traced (zend = 0), xend_cart is set to zero.
         !>
-        !> TODO: Investigate batch coordinate transformation for performance.
-        !> Current implementation uses O(n) individual calls.
+        !> For STL wall hits (wall_hit == 1), xend_cart is set directly from
+        !> wall_hit_cart since that is the authoritative CGAL intersection point.
+        !> This avoids round-trip coordinate conversion errors when from_cart
+        !> fails in ill-conditioned chartmap regions.
 
-        use params, only: ntestpart, zstart, zend
+        use params, only: ntestpart, zstart, zend, wall_hit, wall_hit_cart
         use reference_coordinates, only: ref_coords
 
         real(dp), intent(out) :: xstart_cart(:, :)
@@ -389,14 +390,19 @@ subroutine compute_cartesian_positions(xstart_cart, xend_cart)
         do i = 1, ntestpart
             call ref_coords%evaluate_cart(zstart(1:3, i), xstart_cart(:, i))
 
-            ! Check if particle was traced: zend is set to exactly 0 for untraced
-            ! particles (see simple_main.f90:493 and classification.f90:94)
-            zend_is_zero = all(zend(1:3, i) == 0.0_dp)
-            if (zend_is_zero) then
-                ! Particle not traced - use start position
-                xend_cart(:, i) = xstart_cart(:, i)
+            ! For STL wall hits, use the authoritative CGAL intersection point
+            if (wall_hit(i) == 1) then
+                xend_cart(:, i) = wall_hit_cart(:, i)
             else
-                call ref_coords%evaluate_cart(zend(1:3, i), xend_cart(:, i))
+                ! Check if particle was traced: zend is set to exactly 0 for untraced
+                ! particles (see simple_main.f90:493 and classification.f90:94)
+                zend_is_zero = all(zend(1:3, i) == 0.0_dp)
+                if (zend_is_zero) then
+                    ! Particle not traced - keep xend_cart as zero
+                    xend_cart(:, i) = 0.0_dp
+                else
+                    call ref_coords%evaluate_cart(zend(1:3, i), xend_cart(:, i))
+                end if
             end if
         end do
 

src/simple_main.f90

@@ -542,6 +542,7 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times)
         real(dp) :: x_prev(3), x_cur(3)
         real(dp) :: x_prev_m(3), x_cur_m(3), x_hit_m(3), x_hit(3)
         real(dp) :: normal_m(3), vhat(3), vnorm, cos_inc
+        real(dp) :: segment_length, hit_distance, t_frac
         integer :: it, ierr_orbit, it_final
         integer(8) :: kt
         logical :: passing
@@ -626,6 +627,16 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times)
 
                         if (ierr_from_cart == 0) then
                             call ref_to_integ(u_ref_hit, z(1:3))
+                        else
+                            ! Fallback: linear interpolation of reference coordinates
+                            ! when from_cart fails (ill-conditioned regions of chartmap)
+                            segment_length = sqrt(sum((x_cur_m - x_prev_m)**2))
+                            if (segment_length > 0.0_dp) then
+                                hit_distance = sqrt(sum((x_hit_m - x_prev_m)**2))
+                                t_frac = hit_distance / segment_length
+                                u_ref_hit = u_ref_prev + t_frac * (u_ref_cur - u_ref_prev)
+                                call ref_to_integ(u_ref_hit, z(1:3))
+                            end if
                         end if
 
                         ierr_orbit = 77