Documentation for debugging circle

xcoll/geometry/SimoneNotes.py

@@ -0,0 +1,10 @@
+This will be deleted later
+
+Todo: 
+* Need to fix the slicing box. Deadlock + some bug
+* Connect with Everest & Nuclear interactions. Figure out how to connect path length to scattering through probability
+* test test test
+
+Small todo:
+* Possibly move the crossings folder to scattering part of Xcoll
+* Rename FindRoot to something thats more generic and explains both crossings AND path length 

xcoll/geometry/c_init/__init__.py

@@ -3,5 +3,5 @@
 # Copyright (c) CERN, 2025.                 #
 # ######################################### #
 
-from .c_init import XC_GEOM_EPSILON, XC_GEOM_S_MAX, xo_to_ctypes, xo_to_cnames, PyMethod
+from .c_init import XC_GEOM_EPSILON, XC_GEOM_S_MAX, xo_to_ctypes, xo_to_cnames, PyMethod, define_src
 from .bounding_box import BoundingBox

xcoll/geometry/c_init/bounding_box.py

@@ -42,7 +42,7 @@ class BoundingBox(xo.Struct):
     _extra_c_sources = [
         define_src,
         _pkg_root / 'geometry' / 'c_init' / 'sort.h',
-        _pkg_root / 'geometry' / 'c_init' / 'methods.h',
+       # _pkg_root / 'geometry' / 'c_init' / 'methods.h',
         _pkg_root / 'geometry' / 'c_init' / 'bounding_box.h',
     ]
 

xcoll/geometry/c_init/c_init.py

@@ -13,6 +13,7 @@
 XC_GEOM_ROOT_NEWTON_DERIVATIVE_TOL = 1e-10 # Threshold for small derivative
 XC_GEOM_ROOT_GRID_MAX_INTER = 10           # Maximum number of intervals for grid search
 XC_GEOM_ROOT_GRID_POINTS = 1000            # Number of points to search in grid
+XC_GEOM_SIMPSON_SUBINTERVALS = 500
 
 define_src = f"""
 #ifndef XCOLL_GEOM_DEFINES_H
@@ -50,6 +51,10 @@
 #define XC_GEOM_ROOT_GRID_POINTS {XC_GEOM_ROOT_GRID_POINTS}
 #endif
 
+#ifndef XC_GEOM_SIMPSON_SUBINTERVALS
+#define XC_GEOM_SIMPSON_SUBINTERVALS {XC_GEOM_SIMPSON_SUBINTERVALS}
+#endif
+
 
 #endif /* XCOLL_GEOM_DEFINES_H */
 """

xcoll/geometry/crossings/find_root.py

@@ -2,18 +2,22 @@
 
 import xobjects as xo
 
-from ..segments import LocalSegment
-from ..trajectories import LocalTrajectory
-from ..c_init import define_src, PyMethod
+from ..segments.segment import LocalSegment
+from ..trajectories.trajectory import LocalTrajectory
+from ..c_init import PyMethod, define_src
+
 from ...general import _pkg_root
 
 # Kernel for root finder
 class FindRoot(xo.Struct):
     solution_t    = xo.Float64[:]
-    solution_l    = xo.Float64[:]
+    solution_l    = xo.Float64[:] # REMEMBER TO EDIT SO THAT SOLUTION T IS UPDATED AFTER NUCLEAR
+    guess_t       = xo.Float64[:]
+    guess_l       = xo.Float64[:] 
     converged     = xo.Int8[:]
     num_solutions = xo.Int16
     max_solutions = xo.Int8
+    path_length   = xo.Float64
 
     _depends_on = [LocalTrajectory, LocalSegment]
     _kernels = {'newton': xo.Kernel(
@@ -22,11 +26,19 @@ class FindRoot(xo.Struct):
                                 xo.Arg(xo.ThisClass, name="finder"),
                                 xo.Arg(LocalSegment, name="seg"),
                                 xo.Arg(LocalTrajectory, name="traj"),
-                                xo.Arg(xo.Float64, pointer=True, name="guess_t"),
-                                xo.Arg(xo.Float64, pointer=True, name="guess_l"),
+                                xo.Arg(xo.Float64, pointer=False, name="guess_t"),
+                                xo.Arg(xo.Float64, pointer=False, name="guess_l"),
+                                xo.Arg(xo.Int8, name="num"),
                             ], ret=None),
                 'find_crossing': xo.Kernel(
                             c_name="FindRoot_find_crossing",
+                            args=[
+                                xo.Arg(xo.ThisClass, name="finder"),
+                                xo.Arg(LocalSegment, name="seg"),
+                                xo.Arg(LocalTrajectory, name="traj"),
+                            ], ret=None),
+                'find_path_length': xo.Kernel(
+                            c_name="FindRoot_find_path_length",
                             args=[
                                 xo.Arg(xo.ThisClass, name="finder"),
                                 xo.Arg(LocalSegment, name="seg"),
@@ -35,14 +47,18 @@ class FindRoot(xo.Struct):
     _needs_compilation = True
     _extra_c_sources = [
         define_src,
-        _pkg_root / 'geometry' / 'crossings' / 'find_root.h']
+        _pkg_root / 'geometry' / 'crossings' / 'find_root.h',
+        _pkg_root / 'geometry' / 'crossings' / 'path_length.h']
 
     def __init__(self, **kwargs):
         kwargs.setdefault('max_solutions', 100)
         kwargs.setdefault('num_solutions', 0)
-        kwargs.setdefault('solution_t', np.zeros(kwargs['max_solutions'], dtype=np.float64))
-        kwargs.setdefault('solution_l', np.zeros(kwargs['max_solutions'], dtype=np.float64))
-        kwargs.setdefault('converged',  np.zeros(kwargs['max_solutions'], dtype=np.int8))
+        kwargs.setdefault('solution_t', np.full(kwargs['max_solutions'], 1e21, dtype=np.float64))
+        kwargs.setdefault('solution_l', np.full(kwargs['max_solutions'], 1e21, dtype=np.float64))
+        kwargs.setdefault('guess_t',    np.full(kwargs['max_solutions'], 1e21, dtype=np.float64))
+        kwargs.setdefault('guess_l',    np.full(kwargs['max_solutions'], 1e21, dtype=np.float64))
+        kwargs.setdefault('converged',  np.full(kwargs['max_solutions'], 1e21, dtype=np.int8))
+        kwargs.setdefault('path_length', 0)
         super().__init__(**kwargs)
 
     def __getattr__(self, attr):

xcoll/geometry/crossings/path_length.h

@@ -0,0 +1,68 @@
+// copyright ############################### #
+// This file is part of the Xcoll package.   #
+// Copyright (c) CERN, 2025.                 #
+// ######################################### #
+
+#ifndef XCOLL_GEOM_SIMPSON_H
+#define XCOLL_GEOM_SIMPSON_H
+#include <stdio.h>
+#include <math.h>
+// Simpson's Rule function for numerical integration
+// (b - a) / (3 * subintervals) * (f(l1) + 4f(subinterval1) + 2f(even subinterval) + 4f(odd subinterval)) + ... + f(l2))
+/*gpufun*/
+void simpson(FindRoot finder, LocalTrajectory traj, int32_t subintervals) {
+    double l1 = 0;
+    double l2 = FindRoot_get_solution_l(finder, 0);
+    if (subintervals % 2 != 1) {
+        subintervals++;
+    }
+    double step_size = (l2 - l1) / subintervals;
+    double sum = sqrt(1 + LocalTrajectory_deriv_x(traj, l1)*LocalTrajectory_deriv_x(traj, l1)) + 
+                 sqrt(1 + LocalTrajectory_deriv_x(traj, l2)*LocalTrajectory_deriv_x(traj, l2));  // f(l1) + f(l2)
+    // Add subintervals to the sum
+    for (int i = 1; i < subintervals; i++) {
+        double l = l1 + i * step_size;
+        // Even indices (except the endpoints) are multiplied by 2. Odd indices are multiplied by 4
+        if (i % 2 == 0) {
+            sum += 2.0 * sqrt(1 + LocalTrajectory_deriv_x(traj, l)*LocalTrajectory_deriv_x(traj, l));
+        } else {
+            sum += 4.0 * sqrt(1 + LocalTrajectory_deriv_x(traj, l)*LocalTrajectory_deriv_x(traj, l));
+        }
+    }
+    sum *= step_size / 3.;
+    if (sum < 0) {
+        printf("Warning: Computed path length is negative. Setting to 1e21.\n");
+        fflush(stdout);
+        FindRoot_set_path_length(finder, 1e21);
+        return;
+    }
+    FindRoot_set_path_length(finder, sum);
+    return;
+}
+
+/*gpufun*/
+void find_path_length_analytic(FindRoot finder, LocalTrajectory traj){
+    double l2  = FindRoot_get_solution_l(finder, 0); // We assume we only care for the first solution
+    double s0  = LocalTrajectory_func_s(traj, 0);  // l1 = 0
+    double x0  = LocalTrajectory_func_x(traj, 0);  
+    double s1  = LocalTrajectory_func_s(traj, l2); // l2 = solution for l
+    double x1  = LocalTrajectory_func_x(traj, l2);
+    double path_length = sqrt((s1 - s0)*(s1 - s0) + (x1 - x0)*(x1 - x0));
+    FindRoot_set_path_length(finder, path_length);
+    return;
+    // only for drift for now
+}
+/*gpufun*/
+void FindRoot_find_path_length(FindRoot finder, LocalSegment seg, LocalTrajectory traj){
+    if (LocalSegment_typeid(seg) == LocalSegment_BezierSegment_t){
+        return simpson(finder, traj, XC_GEOM_SIMPSON_SUBINTERVALS);
+    }
+    // TODO: Check with Frederik. I assume we always use the first solution as the second solution might not
+    // even be relevant (say out - in again, then we actually never go in again. IN-OUT -> we change to mcs before out)
+    if (LocalTrajectory_typeid(traj) == LocalTrajectory_DriftTrajectory_t){
+        return find_path_length_analytic(finder, traj);
+    } else {
+        return simpson(finder, traj, XC_GEOM_SIMPSON_SUBINTERVALS);
+    }
+}
+#endif /* XCOLL_GEOM_SIMPSON_H */

xcoll/geometry/segments/bezier.h

@@ -54,7 +54,7 @@ double BezierSegment_deriv_x(BezierSegment seg, double t){
 }
 
 /*gpufun*/
-void BezierSegment_init_bounding_box(BezierSegment seg, BoundingBox box, double t1, double t2){
+void BezierSegment_update_box(BezierSegment seg, double t1, double t2){
     double ts1 = BezierSegment_get__ts1(seg);
     double ts2 = BezierSegment_get__ts2(seg);
     double s1  = BezierSegment_func_s(seg, t1);
@@ -121,7 +121,7 @@ void BezierSegment_init_bounding_box(BezierSegment seg, BoundingBox box, double
     double w = xmax - xmin;  // width of the box
     double sin_tb = 0;       // orientation of the box (angle of length wrt horizontal)
     double cos_tb = 1;
-    BoundingBox_set_params(box, rC, sin_tC, cos_tC, l, w, sin_tb, cos_tb);
+    BoundingBox_set_params(BezierSegment_getp_box(seg), rC, sin_tC, cos_tC, l, w, sin_tb, cos_tb);
 }
 
 
@@ -188,37 +188,6 @@ void BezierSegment_calculate_extrema(BezierSegment seg){
     BezierSegment_set__ex1(seg, BezierSegment_func_x(seg, t1));
     BezierSegment_set__ex2(seg, BezierSegment_func_x(seg, t2));
 }
-
-
-double BezierSegment_prepare_newton(BezierSegment seg, BoundingBox MCSbox, double tol){
-    // Prepare initial guess for Newton-Raphson root finding
-    double org_t1 = BezierSegment_get__t1(seg);
-    double org_t2 = BezierSegment_get__t2(seg);
-    while ((BezierSegment_get__t2(seg) -  BezierSegment_get__t1(seg)) > tol){
-        double t1_old = BezierSegment_get__t1(seg);
-        double t2_old = BezierSegment_get__t2(seg);
-        double t_middle = 0.5 * (BezierSegment_get__t2(seg) + BezierSegment_get__t1(seg));
-
-        // first half
-        BezierSegment_set__t2(seg, t_middle);
-        double overlap_lower = BoundingBox_overlaps(MCSbox, 
-                                                    BezierSegment_getp_box(seg));
-        // second half
-        BezierSegment_set__t2(seg, t2_old);
-        BezierSegment_set__t1(seg, t_middle);
-        double overlap_upper = BoundingBox_overlaps(MCSbox, 
-                                                    BezierSegment_getp_box(seg));
-        if (overlap_lower && !overlap_upper){
-            BezierSegment_set__t1(seg, t1_old);
-            BezierSegment_set__t2(seg, t_middle);
-        }
-    }
-    double guess_t = 0.5 * (BezierSegment_get__t2(seg) + BezierSegment_get__t1(seg));
-    // Reset to original values
-    BezierSegment_set__t1(seg, org_t1);
-    BezierSegment_set__t2(seg, org_t2);
-    return guess_t;
-}
 // /*gpufun*/
 // void _hit_s_bezier(BezierSegment seg, double t, double multiplicity, int8_t* n_hit, double* s){
 //     double s1  = BezierSegment_get_s1(seg);
@@ -232,108 +201,6 @@ double BezierSegment_prepare_newton(BezierSegment seg, BoundingBox MCSbox, doubl
 //     }
 // }
 
-// /*gpufun*/
-// void BezierSegment_crossing_drift(BezierSegment seg, int8_t* n_hit, double* s, double s0, double x0, double xm){
-//     // Get segment data
-//     double s1  = BezierSegment_get_s1(seg);
-//     double x1  = BezierSegment_get_x1(seg);
-//     double s2  = BezierSegment_get_s2(seg);
-//     double x2  = BezierSegment_get_x2(seg);
-//     double cs1 = BezierSegment_get_cs1(seg);
-//     double cx1 = BezierSegment_get_cx1(seg);
-//     double cs2 = BezierSegment_get_cs2(seg);
-//     double cx2 = BezierSegment_get_cx2(seg);
-//     // The Bézier curve is defined by the parametric equations (with t in [0, 1]):
-//     // s(t) = (1-t)^3*s1 + 3(1-t)^2*t*cs1 + 3(1-t)*t^2*cs2 + t^3*s2
-//     // x(t) = (1-t)^3*x1 + 3(1-t)^2*t*cx1 + 3(1-t)*t^2*cx2 + t^3*x2
-//     // Plug the parametric eqs into the drift trajectory x(t) = m*(s(t) - s0) + x0 and solve for t
-//     // The solutions for t (which we get by Cardano's method) are valid if in [0, 1]
-//     double a = (xm*s1 - x1) - (xm*s2 - x2) - 3*(xm*cs1 - cx1) + 3*(xm*cs2 - cx2);
-//     double b = 6*(xm*cs1 - cx1) - 3*(xm*cs2 - cx2) - 3*(xm*s1 - x1);
-//     double c = 3*(xm*s1 - x1) - 3*(xm*cs1 - cx1);
-//     double d = (xm*s0 - x0) - (xm*s1 - x1);
-//     double t;
-//     // Edge cases
-//     if (fabs(a) < XC_GEOM_EPSILON){
-//         if (fabs(b) < XC_GEOM_EPSILON){
-//             if (fabs(c) < XC_GEOM_EPSILON){
-//                 if (fabs(d) < XC_GEOM_EPSILON){
-//                     // The trajectory is on the Bézier curve
-//                     // TODO: This cannot happen because we don't have a cubic trajectory.
-//                     //       Adapt if these ever would be implemented.
-//                     return;
-//                 } else {
-//                     // No solutions
-//                     return;
-//                 }
-//             } else {
-//                 // This is a linear equation
-//                 t = -d/c;
-//                 if (0 <= t && t <= 1){
-//                     _hit_s_bezier(seg, t, 1, n_hit, s);
-//                 }
-//             }
-//         } else {
-//             // This is a quadratic equation
-//             double disc = c*c - 4*b*d;
-//             if (disc < 0){
-//                 // No solutions
-//                 return;
-//             }
-//             for (int8_t i = 0; i < 2; i++) {
-//                 double sgnD = i*2-1; // negative and positive solutions; if multiplicity 2, we add the same solution twice
-//                 t = (-c + sgnD*sqrt(fabs(disc)))/(2*b);
-//                 if (0 <= t && t <= 1){
-//                     _hit_s_bezier(seg, t, 1, n_hit, s);
-//                 }
-//             }
-//         }
-//     } else {
-//         // Full cubic equation. Coefficients for the depressed cubic t^3 + p*t + q = 0:
-//         double p = (3*a*c - b*b)/(3*a*a);
-//         double q = (2*b*b*b - 9*a*b*c + 27*a*a*d)/(27*a*a*a);
-//         double disc = -p*p*p/27 - q*q/4;  // This is the discriminant of the depressed cubic but divided by (4*27)
-//         if (fabs(disc) < XC_GEOM_EPSILON){
-//             if (fabs(p) < XC_GEOM_EPSILON){
-//                 // One real root with multiplicity 3
-//                 t = -b/(3*a);
-//                 if (0 <= t && t <= 1){
-//                     _hit_s_bezier(seg, t, 3, n_hit, s);
-//                 }
-//             } else {
-//                 // Two real roots (one simple and one with multiplicity 2)
-//                 t = 3*q/p - b/(3*a);
-//                 if (0 <= t && t <= 1){
-//                     _hit_s_bezier(seg, t, 1, n_hit, s);
-//                 }
-//                 t = -3*q/(2*p) - b/(3*a);
-//                 if (0 <= t && t <= 1){
-//                     _hit_s_bezier(seg, t, 2, n_hit, s);
-//                 }
-//             }
-//         } else if (disc < 0){
-//             // One real root
-//             t = cbrt(-q/2 + sqrt(fabs(disc))) + cbrt(-q/2 - sqrt(fabs(disc))) - b/(3*a);
-//             if (0 <= t && t <= 1){
-//                 _hit_s_bezier(seg, t, 1, n_hit, s);
-//             }
-//         } else {
-//             // Three real roots
-//             double phi = acos(3*q/(2*p)*sqrt(fabs(3/p)));
-//             t = 2*sqrt(fabs(p/3))*cos(phi/3) - b/(3*a);
-//             if (0 <= t && t <= 1){
-//                 _hit_s_bezier(seg, t, 1, n_hit, s);
-//             }
-//             t = 2*sqrt(fabs(p/3))*cos((phi + 2*M_PI)/3) - b/(3*a);
-//             if (0 <= t && t <= 1){
-//                 _hit_s_bezier(seg, t, 1, n_hit, s);
-//             }
-//             t = 2*sqrt(fabs(p/3))*cos((phi + 4*M_PI)/3) - b/(3*a);
-//             if (0 <= t && t <= 1){
-//                 _hit_s_bezier(seg, t, 1, n_hit, s);
-//             }
-//         }
-//     }
-// }
+
 
 #endif /* XCOLL_GEOM_TRAJ_BEZIER_H */