From 8d48ed6b2549ad4f8879c8682c2f6b14c64e8168 Mon Sep 17 00:00:00 2001
From: NetchAna <vaibhavi.gawas@cern.ch>
Date: Mon, 16 Feb 2026 16:07:22 +0100
Subject: [PATCH 1/4] Add coupled transverse motion to LineSegmentMap with Cij
 matrix

---
 xtrack/beam_elements/elements.py              |  34 ++++
 .../elements_src/linesegmentmap.h             | 149 ++++++++++++++++--
 2 files changed, 174 insertions(+), 9 deletions(-)

diff --git a/xtrack/beam_elements/elements.py b/xtrack/beam_elements/elements.py
index 97158c560..cc34f3a93 100644
--- a/xtrack/beam_elements/elements.py
+++ b/xtrack/beam_elements/elements.py
@@ -3619,6 +3619,11 @@ class LineSegmentMap(BeamElement):
         'alfx': xo.Float64[2],
         'alfy': xo.Float64[2],
 
+        'c11': xo.Float64[2],
+        'c12': xo.Float64[2],
+        'c21': xo.Float64[2],
+        'c22': xo.Float64[2],
+
         'dx': xo.Float64[2],
         'dpx': xo.Float64[2],
         'dy': xo.Float64[2],
@@ -3665,6 +3670,7 @@ class LineSegmentMap(BeamElement):
 
     def __init__(self, length=0., qx=0, qy=0,
             betx=1., bety=1., alfx=0., alfy=0.,
+            c11=0., c12=0., c21=0., c22=0.,
             dx=0., dpx=0., dy=0., dpy=0.,
             x_ref=0.0, px_ref=0.0, y_ref=0.0, py_ref=0.0,
             longitudinal_mode=None,
@@ -3707,6 +3713,18 @@ def __init__(self, length=0., qx=0, qy=0,
         alfy : tuple of length 2 or float
             Vertical alpha function at the entrance and exit of the segment.
             If a float is given, the same value is used for both entrance and exit.
+        c11 : tuple of length 2 or float
+            Coupling matrix element
+            If a float is given, the same value is used for both entrance and exit
+        c12 : tuple of length 2 or float
+            Coupling matrix element
+            If a float is given, the same value is used for both entrance and exit
+        c21 : tuple of length 2 or float
+            Coupling matrix element
+            If a float is given, the same value is used for both entrance and exit
+        c22 : tuple of length 2 or float
+            Coupling matrix element
+            If a float is given, the same value is used for both entrance and exit
         dx : tuple of length 2 or float
             Horizontal dispersion at the entrance and exit of the segment.
             If a float is given, the same value is used for both entrance and exit.
@@ -3963,6 +3981,18 @@ def __init__(self, length=0., qx=0, qy=0,
         if np.isscalar(alfy): alfy = [alfy, alfy]
         else: assert len(alfy) == 2
 
+        if np.isscalar(c11): c11 = [c11, c11]
+        else: assert len(c11) == 2    
+
+        if np.isscalar(c12): c12 = [c12, c12]
+        else: assert len(c12) == 2
+
+        if np.isscalar(c21): c21 = [c21, c21]
+        else: assert len(c21) == 2
+
+        if np.isscalar(c22): c22 = [c22, c22]
+        else: assert len(c22) == 2
+
         if np.isscalar(dx): dx = [dx, dx]
         else: assert len(dx) == 2
 
@@ -3991,6 +4021,10 @@ def __init__(self, length=0., qx=0, qy=0,
         nargs['bety'] = bety
         nargs['alfx'] = alfx
         nargs['alfy'] = alfy
+        nargs['c11'] = c11
+        nargs['c12'] = c12
+        nargs['c21'] = c21
+        nargs['c22'] = c22
         nargs['dx'] = dx
         nargs['dpx'] = dpx
         nargs['dy'] = dy
diff --git a/xtrack/beam_elements/elements_src/linesegmentmap.h b/xtrack/beam_elements/elements_src/linesegmentmap.h
index db6d3b5f3..9abbb421e 100644
--- a/xtrack/beam_elements/elements_src/linesegmentmap.h
+++ b/xtrack/beam_elements/elements_src/linesegmentmap.h
@@ -101,6 +101,19 @@ void transverse_motion(LocalParticle *part0, LineSegmentMapData el){
     int64_t const ndqx = LineSegmentMapData_len_coeffs_dqx(el);
     int64_t const ndqy = LineSegmentMapData_len_coeffs_dqy(el);
 
+    double const c11_0= LineSegmentMapData_get_c11(el,0);
+    double const c12_0= LineSegmentMapData_get_c12(el,0);
+    double const c21_0= LineSegmentMapData_get_c21(el,0);
+    double const c22_0= LineSegmentMapData_get_c22(el,0);
+
+    double const c11_1= LineSegmentMapData_get_c11(el,1);;
+    double const c12_1= LineSegmentMapData_get_c12(el,1);;
+    double const c21_1= LineSegmentMapData_get_c21(el,1);;
+    double const c22_1= LineSegmentMapData_get_c22(el,1);
+
+    double const coupling_0 = c11_0 + c22_0 + c12_0 + c21_0;
+    double const coupling_1 = c11_1 + c22_1 + c12_1 + c21_1;
+
     int64_t detuning;
     double sin_x = 0;
     double cos_x = 0;
@@ -180,18 +193,136 @@ void transverse_motion(LocalParticle *part0, LineSegmentMapData el){
                     )/sqrt_betprod_y;
         double const M11_y = (cos_y-alfy_1*sin_y)/sqrt_betratio_y;
 
-        double const x_out = M00_x*LocalParticle_get_x(part) + M01_x * LocalParticle_get_px(part);
-        double const px_out = M10_x*LocalParticle_get_x(part) + M11_x * LocalParticle_get_px(part);
-        double const y_out = M00_y*LocalParticle_get_y(part) + M01_y * LocalParticle_get_py(part);
-        double const py_out = M10_y*LocalParticle_get_y(part) + M11_y * LocalParticle_get_py(part);
+        if (coupling_0 == 0.0 && coupling_1 == 0.0) {
+            // No coupling matrices, use simple transfer matrix
+            double const x_out = M00_x*LocalParticle_get_x(part) + M01_x * LocalParticle_get_px(part);
+            double const px_out = M10_x*LocalParticle_get_x(part) + M11_x * LocalParticle_get_px(part);
+            double const y_out = M00_y*LocalParticle_get_y(part) + M01_y * LocalParticle_get_py(part);
+            double const py_out = M10_y*LocalParticle_get_y(part) + M11_y * LocalParticle_get_py(part);
+
+            LocalParticle_set_x(part, x_out);
+            LocalParticle_set_px(part, px_out);
+            LocalParticle_set_y(part, y_out);
+            LocalParticle_set_py(part, py_out);
+        }
 
-        LocalParticle_set_x(part, x_out);
-        LocalParticle_set_px(part, px_out);
-        LocalParticle_set_y(part, y_out);
-        LocalParticle_set_py(part, py_out);
-    END_PER_PARTICLE_BLOCK;
+        /* Construct 4x4 symplectic matrix R and its inverse from a 2x2 matrix C */
+        void compute_matrices_V(double const c11, double const c12,
+                          double const c21, double const c22,
+                          double const gamma,
+                          double V[16]) 
+        {
+        /* Initialize R and Rinv to zero */
+            for(int i=0;i<16;i++) {
+            V[i] = 0.0;
+            }
+            /* Identity 2x2 blocks (top-left and bottom-right) */
+            V[0] = gamma;  V[5] = gamma;
+            V[10] = gamma; V[15] = gamma;
+            /* C block (top-right) */
+            V[2] = c11;  V[3] = -c12;
+            V[6] = -c21;  V[7] = c22;
+            /* -C^inv block (bottom-left) */
+            V[8]  = -c22; V[9]  =  -c12;
+            V[12] =  -c21; V[13] = -c11;
+        }
+
+        /* Construct 4x4 symplectic matrix R and its inverse from a 2x2 matrix C */
+        void compute_matrices_Vinv(double const c11, double const c12,
+                          double const c21, double const c22,
+                          double const gamma,
+                          double V_inv[16]) 
+        {
+            /* Initialize R and Rinv to zero */
+            for(int i=0;i<16;i++) {
+                V_inv[i] = 0.0;
+            }
+            /* Identity 2x2 blocks (top-left and bottom-right) */
+            V_inv[0] = gamma;  V_inv[5] = gamma;
+            V_inv[10] = gamma; V_inv[15] = gamma;
+            /* -C block (top-right) */
+            V_inv[2] = -c11;  V_inv[3] =  c12;
+            V_inv[6] =  c21;  V_inv[7] = -c22;
+            /*C^inv block (bottom-left) */
+            V_inv[8]  =  c22; V_inv[9]  =  c12;
+            V_inv[12] =  c21; V_inv[13] =  c11;
+        }
+
+
+        if (coupling_0 != 0.0 || coupling_1 != 0.0) {
+        
+            double const det_0 = c11_0 * c22_0 - c12_0 * c21_0;
+            double const det_1 = c11_1 * c22_1 - c12_1 * c21_1;
+
+            double const gamma_0 = sqrt(1 - det_0);
+            double const gamma_1 = sqrt(1 - det_1);
+
+            double V_1[16], V_0inv[16];
+            double M[16];
+            double T[16];
+            double temp[16];
+
+            compute_matrices_V(c11_1, c12_1, c21_1, c22_1, gamma_1, V_1);  
+            compute_matrices_Vinv(c11_0, c12_0, c21_0, c22_0, gamma_0, V_0inv);
+            
+            //the final transformation is V_1 * M * V_0inv
+
+            for(int i=0;i<16;i++) {
+                M[i] = 0.0;
+            }
+            M[0] = M00_x; M[1] = M01_x; M[2]=0.0; M[3]=0.0;
+            M[4] = M10_x; M[5] = M11_x; M[6]=0.0; M[7]=0.0;
+            M[8] = 0.0; M[9] = 0.0; M[10]=M00_y; M[11]=M01_y;
+            M[12]= 0.0; M[13]= 0.0; M[14]=M10_y; M[15]=M11_y;       
+
+            for(int i=0;i<16;i++) {
+                T[i] = 0.0;
+            }
+
+            for(int i=0;i<16;i++) {
+                temp[i] = 0.0;
+            }
+
+            // temp = M * V_0inv
+            for (int i=0;i<4;i++) {
+                for (int j=0;j<4;j++) {
+                    temp[i*4+j] = 0.0;
+                    for (int k=0;k<4;k++)
+                        temp[i*4+j] += M[i*4+k]*V_0inv[k*4+j];
+                }
+            }
+            // T = V_1 * temp
+            for (int i=0;i<4;i++) {
+                for (int j=0;j<4;j++) {
+                    T[i*4+j] = 0.0;
+                    for (int k=0;k<4;k++)
+                        T[i*4+j] += V_1[i*4+k]*temp[k*4+j];
+                }
+            }
+
+                double const old_x = LocalParticle_get_x(part);
+                double const old_px = LocalParticle_get_px(part);
+                double const old_y = LocalParticle_get_y(part);
+                double const old_py = LocalParticle_get_py(part);
+
+                double const x_out = (
+                    T[0]*old_x + T[1]*old_px + T[2]*old_y + T[3]*old_py);
+                double const px_out = (
+                    T[4]*old_x + T[5]*old_px + T[6]*old_y + T[7]*old_py);
+                double const y_out = (
+                    T[8]*old_x + T[9]*old_px + T[10]*old_y + T[11]*old_py);
+                double const py_out = (
+                    T[12]*old_x + T[13]*old_px + T[14]*old_y + T[15]*old_py);
+
+                LocalParticle_set_x(part, x_out);
+                LocalParticle_set_px(part, px_out);
+                LocalParticle_set_y(part, y_out);
+                LocalParticle_set_py(part, py_out);
+            END_PER_PARTICLE_BLOCK;
+    }
 }
 
+
 GPUFUN
 void longitudinal_motion(LocalParticle *part0,
                          LineSegmentMapData el){

From 9c665168eb5633b92ceabbb6b23b6e5a9acc3f03 Mon Sep 17 00:00:00 2001
From: NetchAna <vaibhavi.gawas@cern.ch>
Date: Wed, 25 Mar 2026 14:47:16 +0100
Subject: [PATCH 2/4] Fix signs of coupling matrix + add tests for coupled
 LineSegmentMap

---
 xtrack/beam_elements/elements_src/.__afs0C09  | 160 +++++
 xtrack/beam_elements/elements_src/.__afs1C1F  | 323 +++++++++
 xtrack/beam_elements/elements_src/.__afs2C4D  | 157 +++++
 xtrack/beam_elements/elements_src/.__afs5693  | 624 +++++++++++++++++
 xtrack/beam_elements/elements_src/.__afs7D11  | 131 ++++
 xtrack/beam_elements/elements_src/.__afsB12B  | 447 ++++++++++++
 xtrack/beam_elements/elements_src/.__afsEE4E  | 422 ++++++++++++
 xtrack/beam_elements/elements_src/.__afsF1AA  | 635 ++++++++++++++++++
 .../elements_src/linesegmentmap.h             |  16 +-
 xtrack/particles/rng_src/.__afsCCC4           |  64 ++
 10 files changed, 2971 insertions(+), 8 deletions(-)
 create mode 100644 xtrack/beam_elements/elements_src/.__afs0C09
 create mode 100644 xtrack/beam_elements/elements_src/.__afs1C1F
 create mode 100644 xtrack/beam_elements/elements_src/.__afs2C4D
 create mode 100644 xtrack/beam_elements/elements_src/.__afs5693
 create mode 100644 xtrack/beam_elements/elements_src/.__afs7D11
 create mode 100644 xtrack/beam_elements/elements_src/.__afsB12B
 create mode 100644 xtrack/beam_elements/elements_src/.__afsEE4E
 create mode 100644 xtrack/beam_elements/elements_src/.__afsF1AA
 create mode 100644 xtrack/particles/rng_src/.__afsCCC4

diff --git a/xtrack/beam_elements/elements_src/.__afs0C09 b/xtrack/beam_elements/elements_src/.__afs0C09
new file mode 100644
index 000000000..26edc6fa1
--- /dev/null
+++ b/xtrack/beam_elements/elements_src/.__afs0C09
@@ -0,0 +1,160 @@
+// copyright ############################### //
+// This file is part of the Xtrack Package.  //
+// Copyright (c) CERN, 2025.                 //
+// ######################################### //
+#ifndef XTRACK_TRACK_MAGNET_EDGE_H
+#define XTRACK_TRACK_MAGNET_EDGE_H
+
+#include <headers/track.h>
+#include <beam_elements/elements_src/track_dipole_edge_linear.h>
+#include <beam_elements/elements_src/track_yrotation.h>
+#include <beam_elements/elements_src/track_wedge.h>
+#include <beam_elements/elements_src/track_mult_fringe.h>
+#include <beam_elements/elements_src/track_dipole_fringe.h>
+
+
+GPUFUN
+void track_magnet_edge_particles(
+    LocalParticle* part0,
+    const int8_t model,  // 0: linear, 1: full, 2: dipole-only, 3: ax ay cancellation
+    const uint8_t is_exit,
+    const double half_gap,
+    const double* knorm,
+    const double* kskew,
+    const int64_t k_order,
+    const double* knl,
+    const double* ksl,
+    const double factor_knl_ksl,
+    const int64_t kl_order,
+    const double ksol,
+    const double x0_solenoid,
+    const double y0_solenoid,
+    const double length,
+    const double face_angle,
+    const double face_angle_feed_down,
+    const double fringe_integral,
+    const double factor_for_backtrack // -1 for backtracking, 1 for forward tracking
+) {
+    double k0 = 0;
+    if (k_order > -1) k0 += knorm[0];
+    if (fabs(length) > 1e-10 && kl_order > -1) k0 += factor_knl_ksl * knl[0] / length;
+
+    // Assume we are coming from or going to a drift
+    if (is_exit) {
+        START_PER_PARTICLE_BLOCK(part0, part);
+            LocalParticle_set_ax(part, 0.);
+            LocalParticle_set_ay(part, 0.);
+        END_PER_PARTICLE_BLOCK;
+    }
+    else {
+        START_PER_PARTICLE_BLOCK(part0, part);
+            LocalParticle_set_ax(part, -0.5 * ksol * (LocalParticle_get_y(part) - y0_solenoid));
+            LocalParticle_set_ay(part, 0.5 * ksol * (LocalParticle_get_x(part) - x0_solenoid));
+        END_PER_PARTICLE_BLOCK;
+    }
+
+    if (model == 0) {  // Linear model
+        // Calculate coefficients for x and y to compute the px and py kicks
+        double r21, r43;
+        compute_dipole_edge_linear_coefficients(
+            k0, face_angle, face_angle_feed_down, half_gap, fringe_integral,
+            &r21, &r43
+        );
+
+        r21 = r21 * factor_for_backtrack;
+        r43 = r43 * factor_for_backtrack;
+
+        START_PER_PARTICLE_BLOCK(part0, part);
+            DipoleEdgeLinear_single_particle(part, r21, r43);
+        END_PER_PARTICLE_BLOCK;
+        return;
+    }
+    else if (model == 1 || model == 2) { // Full model
+
+        if (factor_for_backtrack < 0) {
+            START_PER_PARTICLE_BLOCK(part0, part);
+                LocalParticle_kill_particle(part, -32);
+            END_PER_PARTICLE_BLOCK;
+        }
+
+        uint8_t should_rotate = 0;
+        double sin_ = 0, cos_ = 1, tan_ = 0;
+        if (fabs(face_angle) > 10e-10) {
+            should_rotate = 1;
+            sin_ = sin(face_angle);
+            cos_ = cos(face_angle);
+            tan_ = tan(face_angle);
+        }
+
+        if (is_exit) k0 = -k0;
+
+        #define MAGNET_Y_ROTATE(PART) \
+            if (should_rotate) YRotation_single_particle((PART), -sin_, cos_, -tan_)
+
+        #define MAGNET_DIPOLE_FRINGE(PART) \
+            DipoleFringe_single_particle((PART), fringe_integral, half_gap, k0)
+
+        #define MAGNET_MULTIPOLE_FRINGE(PART) \
+            MultFringe_track_single_particle( \
+                (PART), \
+                knorm, \
+                kskew, \
+                k_order, \
+                knl, \
+                ksl, \
+                kl_order, \
+                length / factor_knl_ksl, \
+                is_exit, \
+                /* min_order */ 1 \
+            );
+        // Above, I use the length to rescale knl and ksl. Here I am relying on
+        // the fact that the length is only used to obtain kn and ks in
+        // MultFringe_track_single_particle. To be remembered if the fringe
+        // model changes!
+
+        #define MAGNET_WEDGE(PART) \
+            if (should_rotate & (k_order >= 0)) Wedge_single_particle((PART), -face_angle, knorm[0])
+
+        #define MAGNET_QUAD_WEDGE(PART) \
+            if (should_rotate & (k_order >= 1)) Quad_wedge_single_particle((PART), -face_angle, knorm[1])
+
+        if (is_exit == 0){ // entry
+            START_PER_PARTICLE_BLOCK(part0, part);
+                MAGNET_Y_ROTATE(part);
+                MAGNET_DIPOLE_FRINGE(part);
+                if (model == 1){
+                    MAGNET_MULTIPOLE_FRINGE(part);
+                }
+                if (model == 1){
+                    MAGNET_QUAD_WEDGE(part);
+                }
+                MAGNET_WEDGE(part);
+            END_PER_PARTICLE_BLOCK;
+        }
+        else { // exit
+            START_PER_PARTICLE_BLOCK(part0, part);
+                MAGNET_WEDGE(part);
+                if (model == 1){
+                    MAGNET_QUAD_WEDGE(part);
+                }
+                if (model == 1){
+                    MAGNET_MULTIPOLE_FRINGE(part);
+                }
+                MAGNET_DIPOLE_FRINGE(part);
+                MAGNET_Y_ROTATE(part);
+            END_PER_PARTICLE_BLOCK;
+        }
+
+        #undef MAGNET_Y_ROTATE
+        #undef MAGNET_DIPOLE_FRINGE
+        #undef MAGNET_MULTIPOLE_FRINGE
+        #undef MAGNET_WEDGE
+        #undef MAGNET_QUAD_WEDGE
+    }
+    else if (model == 3) { // only ax ay cancellation (already done above)
+        // do nothing
+    }
+    // If model is not 0 or 1, do nothing
+}
+
+#endif // XTRACK_TRACK_MAGNET_EDGE_H
diff --git a/xtrack/beam_elements/elements_src/.__afs1C1F b/xtrack/beam_elements/elements_src/.__afs1C1F
new file mode 100644
index 000000000..c400d0797
--- /dev/null
+++ b/xtrack/beam_elements/elements_src/.__afs1C1F
@@ -0,0 +1,323 @@
+// copyright ############################### //
+// This file is part of the Xtrack Package.  //
+// Copyright (c) CERN, 2023.                 //
+// ######################################### //
+#ifndef XTRACK_TRACK_MAGNET_KICK_H
+#define XTRACK_TRACK_MAGNET_KICK_H
+
+#include <headers/track.h>
+
+
+GPUFUN
+void kick_simple_single_particle(
+    LocalParticle* part,
+    int64_t order,
+    double inv_factorial,
+    const double* knl,
+    const double* ksl,
+    double factor,
+    double kick_weight
+);
+
+
+GPUFUN
+void track_magnet_kick_single_particle(
+    LocalParticle* part,
+    double length,
+    int64_t order,
+    double inv_factorial_order,
+    GPUGLMEM const double* knl,
+    GPUGLMEM const double* ksl,
+    double const factor_knl_ksl,
+    double kick_weight,
+    double k0,
+    double k1,
+    double k2,
+    double k3,
+    double k0s,
+    double k1s,
+    double k2s,
+    double k3s,
+    double h,
+    double hxl,
+    double k0_h_correction,
+    double k1_h_correction,
+    uint8_t rot_frame
+){
+
+    double const chi = LocalParticle_get_chi(part);
+    double const x = LocalParticle_get_x(part);
+    double const y = LocalParticle_get_y(part);
+
+    double knl_main[4] = {k0, k1, k2, k3};
+    double ksl_main[4] = {k0s, k1s, k2s, k3s};
+
+    for (int index = 0; index < 4; index++) {
+        knl_main[index] = knl_main[index] * length;
+        ksl_main[index] = ksl_main[index] * length;
+    }
+
+    // multipolar kick
+    kick_simple_single_particle(
+        part,
+        order,
+        inv_factorial_order,
+        knl,
+        ksl,
+        factor_knl_ksl,
+        kick_weight
+    );
+
+    kick_simple_single_particle(
+        part,
+        /* order */ 3,
+        /* inv_factorial_order */ 1. / (3 * 2),
+        knl_main,
+        ksl_main,
+        /* factor_knl_ksl */ 1,
+        kick_weight
+    );
+
+    // Correct for the curvature
+    double dpx = 0;
+    double dpy = 0;
+    double dzeta = 0;
+
+    if (rot_frame) {
+        double const hl = h * length * kick_weight + hxl * kick_weight;
+        dpx += hl * (1. + LocalParticle_get_delta(part));
+        double const rv0v = 1./LocalParticle_get_rvv(part);
+        dzeta += -rv0v * hl * x;
+    }
+
+    double htot = h;
+    if (length != 0) {
+        htot += hxl / length;
+    }
+
+    // Correct for the curvature
+    // k0h correction can be computed from this term in the hamiltonian
+    // H = 1/2 h k0 x^2
+    // (see MAD 8 physics manual, eq. 5.15, and apply Hamilton's eq. dp/ds = -dH/dx)
+    double k0l_mult = 0;
+    if (order >= 0) {
+        k0l_mult = knl[0] * factor_knl_ksl;
+    }
+    dpx += -chi * (k0_h_correction  *length + k0l_mult) * kick_weight * htot * x;
+
+    // k1h correction can be computed from this term in the hamiltonian
+    // H = 1/3 hk1 x^3 - 1/2 hk1 xy^2
+    // (see MAD 8 physics manual, eq. 5.15, and apply Hamilton's eq. dp/ds = -dH/dx)
+    double k1l_mult = 0;
+    if (order >= 1) {
+        k1l_mult = knl[1] * factor_knl_ksl;
+    }
+    dpx += htot * chi * (k1_h_correction * length + k1l_mult) * kick_weight * (-x * x + 0.5 * y * y);
+    dpy += htot * chi * (k1_h_correction * length  + k1l_mult) * kick_weight * x * y;
+
+    LocalParticle_add_to_px(part, dpx);
+    LocalParticle_add_to_py(part, dpy);
+    LocalParticle_add_to_zeta(part, dzeta);
+
+}
+
+
+
+GPUFUN
+uint8_t kick_is_inactive(
+    int64_t order,
+    GPUGLMEM const double* knl,
+    GPUGLMEM const double* ksl,
+    double k0,
+    double k1,
+    double k2,
+    double k3,
+    double k0s,
+    double k1s,
+    double k2s,
+    double k3s,
+    double h
+){
+    if (h != 0) return 0;
+    if (k0 != 0) return 0;
+    if (k1 != 0) return 0;
+    if (k2 != 0) return 0;
+    if (k3 != 0) return 0;
+    if (k0s != 0) return 0;
+    if (k1s != 0) return 0;
+    if (k2s != 0) return 0;
+    if (k3s != 0) return 0;
+
+    for (int index = order; index >= 0; index--) {
+        if (knl[index] != 0) return 0;
+        if (ksl[index] != 0) return 0;
+    }
+
+    return 1;
+
+}
+
+GPUFUN
+void kick_simple_single_coordinates(
+    double const x,
+    double const y,
+    double const chi,
+    int64_t order,
+    double inv_factorial,
+    const double* knl,
+    const double* ksl,
+    double factor,
+    double kick_weight,
+    double *dpx,
+    double *dpy
+) {
+
+    int64_t index = order;
+
+    double dpx_mul = chi * knl[index] * factor * inv_factorial;
+    double dpy_mul = chi * ksl[index] * factor * inv_factorial;
+
+    while( index > 0 )
+    {
+        double const zre = dpx_mul * x - dpy_mul * y;
+        double const zim = dpx_mul * y + dpy_mul * x;
+
+        inv_factorial *= index;
+        index -= 1;
+
+        double this_knl = chi * knl[index] * factor;
+        double this_ksl = chi * ksl[index] * factor;
+
+        dpx_mul = this_knl * inv_factorial + zre;
+        dpy_mul = this_ksl * inv_factorial + zim;
+    }
+
+    dpx_mul = -dpx_mul; // rad
+
+    *dpx = kick_weight * dpx_mul;
+    *dpy = kick_weight * dpy_mul;
+}
+
+
+GPUFUN
+void kick_simple_single_particle(
+    LocalParticle* part,
+    int64_t order,
+    double inv_factorial,
+    const double* knl,
+    const double* ksl,
+    double factor,
+    double kick_weight
+) {
+    double const chi = LocalParticle_get_chi(part);
+    double const x = LocalParticle_get_x(part);
+    double const y = LocalParticle_get_y(part);
+
+    double dpx, dpy;
+
+    kick_simple_single_coordinates(
+        x,
+        y,
+        chi,
+        order,
+        inv_factorial,
+        knl,
+        ksl,
+        factor,
+        kick_weight,
+        &dpx,
+        &dpy);
+
+    LocalParticle_add_to_px(part, dpx);
+    LocalParticle_add_to_py(part, dpy);
+}
+
+GPUFUN
+void evaluate_field_from_strengths(
+    double const p0c,
+    double const q0,
+    double const x,
+    double const y,
+    double length,
+    int64_t order,
+    double inv_factorial_order,
+    GPUGLMEM const double* knl,
+    GPUGLMEM const double* ksl,
+    double const factor_knl_ksl,
+    double k0,
+    double k1,
+    double k2,
+    double k3,
+    double k0s,
+    double k1s,
+    double k2s,
+    double k3s,
+    double ks,
+    double dks_ds,
+    double x0_solenoid,
+    double y0_solenoid,
+    double *Bx_T,
+    double *By_T,
+    double *Bz_T
+){
+    if (length == 0.0) {
+        *Bx_T = 0.0;
+        *By_T = 0.0;
+        *Bz_T = 0.0;
+        return;
+    }
+
+    double knl_main[4] = {k0, k1, k2, k3};
+    double ksl_main[4] = {k0s, k1s, k2s, k3s};
+
+    for (int index = 0; index < 4; index++) {
+        knl_main[index] = knl_main[index] * length;
+        ksl_main[index] = ksl_main[index] * length;
+    }
+
+    // multipolar kick
+    double dpx_mul = 0.;
+    double dpy_mul = 0.;
+    kick_simple_single_coordinates(
+        x,
+        y,
+        1., // chi
+        order,
+        inv_factorial_order,
+        knl,
+        ksl,
+        factor_knl_ksl,
+        1., // kick_weight
+        &dpx_mul,
+        &dpy_mul);
+
+
+    // main kick
+    double dpx_main=0.;
+    double dpy_main=0.;
+    kick_simple_single_coordinates(
+        x,
+        y,
+        1., // chi
+        3, // order
+        1. / (3 * 2), //inv_factorial_order
+        knl_main,
+        ksl_main,
+        1, // factor_knl_ksl,
+        1., // kick_weight
+        &dpx_main,
+        &dpy_main);
+
+    double const dpx = dpx_mul + dpx_main;
+    double const dpy = dpy_mul + dpy_main;
+
+    double const brho_0 = p0c / C_LIGHT / q0; // [T m]
+
+
+    *Bx_T = dpy * brho_0 / length - 0.5 * dks_ds * brho_0 * (x - x0_solenoid); // [T]
+    *By_T = -dpx * brho_0 / length - 0.5 * dks_ds * brho_0 * (y - y0_solenoid); // [T]
+    *Bz_T = ks * brho_0; // [T]
+
+}
+
+#endif
\ No newline at end of file
diff --git a/xtrack/beam_elements/elements_src/.__afs2C4D b/xtrack/beam_elements/elements_src/.__afs2C4D
new file mode 100644
index 000000000..7a6c2624f
--- /dev/null
+++ b/xtrack/beam_elements/elements_src/.__afs2C4D
@@ -0,0 +1,157 @@
+#ifndef XTRACK_TRACK_MAGNET_CONFIGURE_H
+#define XTRACK_TRACK_MAGNET_CONFIGURE_H
+
+#define H_TOLERANCE (1e-8)
+
+GPUFUN
+void configure_tracking_model(
+    int64_t model,
+    double k0,
+    double k1,
+    double h,
+    double ks,
+    double* k0_drift,
+    double* k1_drift,
+    double* h_drift,
+    double* ks_drift,
+    double* k0_kick,
+    double* k1_kick,
+    double* h_kick,
+    double* k0_h_correction,
+    double* k1_h_correction,
+    int8_t* kick_rot_frame,
+    int8_t* out_drift_model
+){
+
+    // model = 0 or 1 : adaptive
+    // model = 2: bend-kick-bend
+    // model = 3: rot-kick-rot
+    // model = 4: mat-kick-mat (previously called `expanded`)
+    // model = 5: drift-kick-drift-exact
+    // model = 6: drift-kick-drift-expanded
+    // model = -1: kick only (not exposed in python)
+    // model = -2: sol-kick-sol (not exposed in python)
+
+    if (model==1){
+        model = 3; // backward compatibility
+    }
+
+    int8_t h_is_zero = (fabs(h) < H_TOLERANCE);
+    int64_t drift_model = 0;
+
+    if (model == 2){ // bend-kick-bend
+        if (h_is_zero){
+            drift_model = 5; // bend without h
+        }
+        else{
+            drift_model = 4; // bend with h
+        }
+    }
+    else if(model == 3){ // rot-kick-rot
+        if (h_is_zero){
+            drift_model = 1; // drift exact
+        }
+        else{
+            drift_model = 2; // polar drift
+        }
+    }
+    else if(model == 4){ // mat-kick-mat
+        drift_model = 3; // expanded
+    }
+    else if(model == 5){ // drift-kick-drift-exact
+        drift_model = 1; // drift exact
+    }
+    else if(model == 6){ // drift-kick-drift-expanded
+        drift_model = 0; // drift expanded
+    }
+    else if(model == -1){ // kick only
+        drift_model = -1;
+    }
+    else if(model == -2){ // sol-kick-sol
+        drift_model = 6; // solenoid
+    }
+    else{
+        // This should never happen, but just in case
+        drift_model = 99999999;
+    }
+
+    if (drift_model == -1 || drift_model == 0 || drift_model == 1){ // drift expanded, drift exact, kick only
+        *k0_drift = 0.0;
+        *k1_drift = 0.0;
+        *h_drift = 0.0;
+        *ks_drift = 0.0;
+        *k0_kick = k0;
+        *k1_kick = k1;
+        *h_kick = h;
+        *k0_h_correction = k0;
+        *k1_h_correction = k1;
+        *kick_rot_frame = 1;
+    }
+    else if (drift_model == 2){ // polar drift
+        *k0_drift = 0.0;
+        *k1_drift = 0.0;
+        *h_drift = h;
+        *ks_drift = 0.0;
+        *k0_kick = k0;
+        *k1_kick = k1;
+        *h_kick = h;
+        *k0_h_correction = k0;
+        *k1_h_correction = k1;
+        *kick_rot_frame = 0;
+    }
+    else if (drift_model == 3){ // expanded dipole-quadrupole
+        *k0_drift = k0;
+        *k1_drift = k1;
+        *h_drift = h;
+        *ks_drift = 0.0;
+        *k0_kick = 0.0;
+        *k1_kick = 0.0;
+        *h_kick = h;
+        *k0_h_correction = 0.;
+        *k1_h_correction = k1;
+        *kick_rot_frame = 0;
+    }
+    else if (drift_model == 4){ // bend with h
+        *k0_drift = k0;
+        *k1_drift = 0.0;
+        *h_drift = h;
+        *ks_drift = 0.0;
+        *k0_kick = 0.0;
+        *k1_kick = k1;
+        *h_kick = h;
+        *k0_h_correction = 0.;
+        *k1_h_correction = k1;
+        *kick_rot_frame = 0;
+    }
+    else if (drift_model == 5){ // bend without h
+        *k0_drift = k0;
+        *k1_drift = 0.0;
+        *h_drift = 0.0;
+        *ks_drift = 0.0;
+        *k0_kick = 0.0;
+        *k1_kick = k1;
+        *h_kick = 0.0;
+        *k0_h_correction = 0.;
+        *k1_h_correction = 0.;
+        *kick_rot_frame = 0;
+    }
+    else if (drift_model == 6){ // solenoid
+        *k0_drift = 0.0;
+        *k1_drift = 0.0;
+        *h_drift = 0.0;
+        *ks_drift = ks;
+        *k0_kick = k0;
+        *k1_kick = k1;
+        *h_kick = h;
+        *k0_h_correction = k0;
+        *k1_h_correction = k1;
+        *kick_rot_frame = 1;
+    }
+
+
+    *out_drift_model = drift_model;
+}
+
+#undef H_TOLERANCE
+
+#endif // XTRACK_TRACK_MAGNET_CONFIGURE_H
\ No newline at end of file
diff --git a/xtrack/beam_elements/elements_src/.__afs5693 b/xtrack/beam_elements/elements_src/.__afs5693
new file mode 100644
index 000000000..299a4799a
--- /dev/null
+++ b/xtrack/beam_elements/elements_src/.__afs5693
@@ -0,0 +1,624 @@
+// copyright ############################### //
+// This file is part of the Xtrack Package.  //
+// Copyright (c) CERN, 2023.                 //
+// ######################################### //
+#ifndef XTRACK_TRACK_MAGNET_H
+#define XTRACK_TRACK_MAGNET_H
+
+#include <headers/track.h>
+#include <beam_elements/elements_src/track_magnet_kick.h>
+#include <beam_elements/elements_src/track_magnet_drift.h>
+#include <beam_elements/elements_src/track_magnet_edge.h>
+#include <beam_elements/elements_src/track_magnet_configure.h>
+
+#ifndef XTRACK_MULTIPOLE_NO_SYNRAD
+#include <beam_elements/elements_src/track_magnet_radiation.h>
+#endif
+
+#define H_TOLERANCE (1e-8)
+
+#ifndef VSWAP
+    #define VSWAP(a, b) { double tmp = a; a = b; b = tmp; }
+#endif
+
+
+GPUFUN
+void track_magnet_body_single_particle(
+    LocalParticle* part,
+    const double length,
+    const int64_t order,
+    const double inv_factorial_order,
+    GPUGLMEM const double* knl,
+    GPUGLMEM const double* ksl,
+    const double factor_knl_ksl,
+    const int64_t num_multipole_kicks,
+    const int8_t kick_rot_frame,
+    const int8_t drift_model,
+    const int8_t integrator,
+    const double k0_drift,
+    const double k1_drift,
+    const double ks_drift,
+    const double h_drift,
+    const double k0_kick,
+    const double k1_kick,
+    const double h_kick,
+    const double hxl,
+    const double k0_h_correction,
+    const double k1_h_correction,
+    const double k2,
+    const double k3,
+    const double k0s,
+    const double k1s,
+    const double k2s,
+    const double k3s,
+    const double dks_ds,
+    const double x0_solenoid,
+    const double y0_solenoid,
+    const int64_t radiation_flag,
+    const int64_t spin_flag,
+    SynchrotronRadiationRecordData radiation_record,
+    double* dp_record_exit,
+    double* dpx_record_exit,
+    double* dpy_record_exit
+) {
+
+    #define MAGNET_KICK(part, weight) \
+        track_magnet_kick_single_particle(\
+            part, length, order, inv_factorial_order, \
+            knl, ksl, factor_knl_ksl, (weight), \
+            k0_kick, k1_kick, k2, k3, k0s, k1s, k2s, k3s, h_kick,\
+            hxl, k0_h_correction, k1_h_correction, kick_rot_frame\
+        )
+
+    #define MAGNET_DRIFT(part, dlength) \
+        track_magnet_drift_single_particle(\
+            part, (dlength), k0_drift, k1_drift, ks_drift, h_drift,\
+            x0_solenoid, y0_solenoid, drift_model\
+        )
+
+    #ifdef XTRACK_MULTIPOLE_NO_SYNRAD
+        #define WITH_RADIATION(ll, code)\
+        {\
+            code;\
+        }
+    #else
+        #define WITH_RADIATION(ll, code) \
+        { \
+            double const old_x = LocalParticle_get_x(part); \
+            double const old_y = LocalParticle_get_y(part); \
+            double const old_zeta = LocalParticle_get_zeta(part); \
+            double const old_kin_px = LocalParticle_get_px(part) - LocalParticle_get_ax(part);\
+            double const old_kin_py = LocalParticle_get_py(part) - LocalParticle_get_ay(part);\
+            code; \
+            if ((radiation_flag || spin_flag) && length > 0){ \
+                double h_for_rad = h_kick + hxl / length; \
+                if (fabs(h_drift) > 0){ h_for_rad = h_drift; } \
+                double Bx_T, By_T, Bz_T; \
+                double const p0c = LocalParticle_get_p0c(part); \
+                double const q0 = LocalParticle_get_q0(part); \
+                double const new_x = LocalParticle_get_x(part); \
+                double const new_y = LocalParticle_get_y(part); \
+                double const new_kin_px = LocalParticle_get_px(part) - LocalParticle_get_ax(part);\
+                double const new_kin_py = LocalParticle_get_py(part) - LocalParticle_get_ay(part);\
+                double const mean_x = 0.5 * (old_x + new_x); \
+                double const mean_y = 0.5 * (old_y + new_y); \
+                double const mean_kin_px = 0.5 * (old_kin_px + new_kin_px); \
+                double const mean_kin_py = 0.5 * (old_kin_py + new_kin_py); \
+                evaluate_field_from_strengths( \
+                    p0c, \
+                    q0, \
+                    mean_x, \
+                    mean_y, \
+                    length, \
+                    order, \
+                    inv_factorial_order, \
+                    knl, \
+                    ksl, \
+                    factor_knl_ksl, \
+                    k0_drift + k0_kick, \
+                    k1_drift + k1_kick, \
+                    k2, \
+                    k3, \
+                    k0s, \
+                    k1s, \
+                    k2s, \
+                    k3s, \
+                    ks_drift, \
+                    dks_ds, \
+                    x0_solenoid, \
+                    y0_solenoid, \
+                    &Bx_T, \
+                    &By_T, \
+                    &Bz_T \
+                ); \
+                double const dzeta = LocalParticle_get_zeta(part) - old_zeta; \
+                double const rvv = LocalParticle_get_rvv(part); \
+                double l_path = rvv * (ll - dzeta); \
+                if (spin_flag){ \
+                    magnet_spin( \
+                        part, \
+                        Bx_T, \
+                        By_T, \
+                        Bz_T, \
+                        h_for_rad, \
+                        ll, \
+                        l_path); \
+                } \
+                if (radiation_flag){ \
+                    double const B_perp_T = compute_b_perp_mod( \
+                        mean_kin_px, \
+                        mean_kin_py, \
+                        LocalParticle_get_delta(part), \
+                        Bx_T, \
+                        By_T, \
+                        Bz_T); \
+                    magnet_radiation( \
+                        part, \
+                        B_perp_T, \
+                        ll, \
+                        l_path, \
+                        radiation_flag, \
+                        radiation_record, \
+                        dp_record_exit, dpx_record_exit, dpy_record_exit \
+                    ); \
+                } \
+            }\
+        }
+    #endif
+
+    if (num_multipole_kicks == 0 && k0_kick == 0 && k1_kick == 0 && h_kick == 0) { //only drift
+        WITH_RADIATION(length,
+            MAGNET_DRIFT(part, length);
+        )
+    }
+    else{
+
+        
+    // START GENERATED INTEGRATION CODE
+
+    if (integrator == 1){ // TEAPOT
+
+        WITH_RADIATION(length,
+            const double kick_weight = 1. / num_multipole_kicks;
+            double edge_drift_weight = 0.5;
+            double inside_drift_weight = 0;
+            if (num_multipole_kicks > 1) {
+                edge_drift_weight = 1. / (2 * (1 + num_multipole_kicks));
+                inside_drift_weight = (
+                    ((double) num_multipole_kicks)
+                        / ((double)(num_multipole_kicks*num_multipole_kicks) - 1));
+            }
+
+            MAGNET_DRIFT(part, edge_drift_weight*length);
+            for (int i_kick=0; i_kick<num_multipole_kicks - 1; i_kick++) {
+                MAGNET_KICK(part, kick_weight);
+                MAGNET_DRIFT(part, inside_drift_weight*length);
+            }
+            MAGNET_KICK(part, kick_weight);
+            MAGNET_DRIFT(part, edge_drift_weight*length);
+        )
+
+    }
+    else if (integrator==3){ // uniform
+
+        const double kick_weight = 1. / num_multipole_kicks;
+        const double drift_weight = kick_weight;
+
+        for (int i_kick=0; i_kick<num_multipole_kicks; i_kick++) {
+            WITH_RADIATION(drift_weight*length,
+                MAGNET_DRIFT(part, 0.5*drift_weight*length);
+                MAGNET_KICK(part, kick_weight);
+                MAGNET_DRIFT(part, 0.5*drift_weight*length);
+            )
+        }
+
+    }
+    else if (integrator==2){ // YOSHIDA 4
+
+        const int64_t n_kicks_yoshida = 7;
+        const int64_t num_slices = (num_multipole_kicks / n_kicks_yoshida
+                                + (num_multipole_kicks % n_kicks_yoshida != 0));
+
+        const double slice_length = length / (num_slices);
+        const double kick_weight = 1. / num_slices;
+        const double d_yoshida[] =
+                     // From MAD-NG
+                     {3.922568052387799819591407413100e-01,
+                      5.100434119184584780271052295575e-01,
+                      -4.710533854097565531482416645304e-01,
+                      6.875316825251809316199569366290e-02};
+                    //  {0x1.91abc4988937bp-2, 0x1.052468fb75c74p-1, // same in hex
+                    //  -0x1.e25bd194051b9p-2, 0x1.199cec1241558p-4 };
+                    //  {1/8.0, 1/8.0, 1/8.0, 1/8.0}; // Uniform, for debugging
+        const double k_yoshida[] =
+                     // From MAD-NG
+                     {7.845136104775599639182814826199e-01,
+                      2.355732133593569921359289764951e-01,
+                      -1.177679984178870098432412305556e+00,
+                      1.315186320683906284756403692882e+00};
+                    //  {0x1.91abc4988937bp-1, 0x1.e2743579895b4p-3, // same in hex
+                    //  -0x1.2d7c6f7933b93p+0, 0x1.50b00cfb7be3ep+0 };
+                    //  {1/7.0, 1/7.0, 1/7.0, 1/7.0}; // Uniform, for debugging
+
+            for (int ii = 0; ii < num_slices; ii++) {
+                WITH_RADIATION(slice_length,
+                    MAGNET_DRIFT(part, slice_length * d_yoshida[0]);
+                    MAGNET_KICK(part, kick_weight * k_yoshida[0]);
+                    MAGNET_DRIFT(part, slice_length * d_yoshida[1]);
+                    MAGNET_KICK(part, kick_weight * k_yoshida[1]);
+                    MAGNET_DRIFT(part, slice_length * d_yoshida[2]);
+                    MAGNET_KICK(part, kick_weight * k_yoshida[2]);
+                    MAGNET_DRIFT(part, slice_length * d_yoshida[3]);
+                    MAGNET_KICK(part, kick_weight * k_yoshida[3]);
+                    MAGNET_DRIFT(part, slice_length * d_yoshida[3]);
+                    MAGNET_KICK(part, kick_weight * k_yoshida[2]);
+                    MAGNET_DRIFT(part, slice_length * d_yoshida[2]);
+                    MAGNET_KICK(part, kick_weight * k_yoshida[1]);
+                    MAGNET_DRIFT(part, slice_length * d_yoshida[1]);
+                    MAGNET_KICK(part, kick_weight * k_yoshida[0]);
+                    MAGNET_DRIFT(part, slice_length * d_yoshida[0]);
+                ) // WITH_RADIATION
+            }
+    } // integrator if
+
+    // END GENERATED INTEGRATION CODE
+
+
+    }
+
+    #undef MAGNET_KICK
+    #undef MAGNET_DRIFT
+    #undef WITH_RADIATION
+
+}
+
+
+GPUFUN
+void track_magnet_particles(
+    double const weight,
+    LocalParticle* part0,
+    double length,
+    int64_t order,
+    double inv_factorial_order,
+    GPUGLMEM const double* knl,
+    GPUGLMEM const double* ksl,
+    int64_t num_multipole_kicks,
+    int8_t model,
+    int8_t default_model,
+    int8_t integrator,
+    int8_t default_integrator,
+    int64_t radiation_flag,
+    int64_t radiation_flag_parent,
+    SynchrotronRadiationRecordData radiation_record,
+    double delta_taper,
+    double h,
+    double hxl,
+    double k0,
+    double k1,
+    double k2,
+    double k3,
+    double k0s,
+    double k1s,
+    double k2s,
+    double k3s,
+    double ks,
+    double dks_ds,
+    double x0_solenoid,
+    double y0_solenoid,
+    int64_t rbend_model, // -1: not used, 0: auto, 1: curved body, 2: straight body
+    int64_t rbend_compensate_sagitta,
+    double rbend_shift,
+    double rbend_angle_diff,
+    double length_straight,
+    int64_t body_active,
+    int64_t edge_entry_active,
+    int64_t edge_exit_active,
+    int64_t edge_entry_model,
+    int64_t edge_exit_model,
+    double edge_entry_angle,
+    double edge_exit_angle,
+    double edge_entry_angle_fdown,
+    double edge_exit_angle_fdown,
+    double edge_entry_fint,
+    double edge_exit_fint,
+    double edge_entry_hgap,
+    double edge_exit_hgap
+) {
+
+    double factor_knl_ksl = 1.0;
+
+    // Used only for rbend-straight-body
+    double theta_in = 0.;
+    double theta_out = 0.;
+    double cos_theta_in = 1.;
+    double sin_theta_in = 0.;
+    double cos_theta_out = 1.;
+    double sin_theta_out = 0.;
+    double length_curved = 0.;
+    double x0_mid = 0.;
+    double x0_in = 0.;
+    double x0_out = 0.;
+
+    if (rbend_model == 0){
+        // auto mode, curved body
+        rbend_model = 1;
+    }
+
+    double const angle = h * length;
+    if (rbend_model == 1){
+        // curved body
+        edge_entry_angle += (angle - rbend_angle_diff) / 2.0;
+        edge_exit_angle += (angle + rbend_angle_diff) / 2.0;
+    }
+    else if (rbend_model == 2){
+        // straight body
+        theta_in = (angle - rbend_angle_diff) / 2.0;
+        if (fabs(theta_in) > 1e-10){
+            sin_theta_in = sin(theta_in);
+            cos_theta_in = cos(theta_in);
+        }
+        theta_out = (angle + rbend_angle_diff) / 2.0;
+        if (fabs(theta_out) > 1e-10){
+            sin_theta_out = sin(theta_out);
+            cos_theta_out = cos(theta_out);
+        }
+
+        length_curved = length;
+        length = length_straight;
+
+        x0_mid -= rbend_shift;
+
+        if (rbend_compensate_sagitta && fabs(angle) > 1e-10){
+            // shift by half the sagitta
+            double cos_rbha = cos(angle / 2.);
+            x0_mid += 0.5 / h * (1 - cos_rbha);
+        }
+
+        x0_in = x0_mid;
+        x0_out = x0_mid;
+        if (fabs(angle) > 1e-10){
+            double const px0_in = sin(theta_in);
+            double const px0_mid = px0_in - h * length_straight / 2;
+            double const sqrt_mid = sqrt(1 - px0_mid * px0_mid);
+            x0_in -= 1/h *(sqrt_mid - cos_theta_in);
+            x0_out += 1/h * (cos_theta_out - sqrt_mid);
+        }
+        ;
+        h = 0; // treat magnet as straight
+        // We are entering the fringe with an angle, the linear fringe
+        // needs to come from the  expansion around the right angle
+        edge_entry_angle_fdown += theta_in;
+        edge_exit_angle_fdown += theta_out;
+    }
+
+    double core_length, core_length_curved,factor_knl_ksl_body,
+           factor_knl_ksl_edge,factor_backtrack_edge;
+    // Backtracking
+    if (LocalParticle_check_track_flag(part0, XS_FLAG_BACKTRACK)) {
+        core_length = -length * weight;
+        core_length_curved = -length_curved * weight;
+        factor_knl_ksl_body = -factor_knl_ksl * weight;
+        factor_knl_ksl_edge = factor_knl_ksl; // Edge has a specific factor for backtracking
+        factor_backtrack_edge = -1.;
+        hxl = -hxl;
+        VSWAP(edge_entry_active, edge_exit_active);
+        VSWAP(edge_entry_model, edge_exit_model);
+        VSWAP(edge_entry_angle, edge_exit_angle);
+        VSWAP(edge_entry_angle_fdown, edge_exit_angle_fdown);
+        VSWAP(edge_entry_fint, edge_exit_fint);
+        VSWAP(edge_entry_hgap, edge_exit_hgap);
+        VSWAP(theta_in, theta_out);
+        VSWAP(cos_theta_in, cos_theta_out);
+        VSWAP(sin_theta_in, sin_theta_out);
+        theta_in = -theta_in;
+        theta_out = -theta_out;
+        sin_theta_in = -sin_theta_in;
+        sin_theta_out = -sin_theta_out;
+        VSWAP(x0_in, x0_out);
+    }
+    else {
+        core_length = length * weight;
+        core_length_curved = length_curved * weight;
+        factor_knl_ksl_body = factor_knl_ksl * weight;
+        factor_knl_ksl_edge = factor_knl_ksl;
+        factor_backtrack_edge = 1.;
+    }
+
+    #ifndef XTRACK_MULTIPOLE_NO_SYNRAD
+        if (radiation_flag == 10){ // from parent
+            radiation_flag = radiation_flag_parent;
+        }
+    #endif
+
+    // Tapering
+    if (LocalParticle_check_track_flag(part0, XS_FLAG_SR_TAPER)){
+        part0->ipart = 0;
+        delta_taper = LocalParticle_get_delta(part0); // I can use part0 because
+                                                      // there is only one particle
+                                                      // when doing the tapering
+    }
+
+    #ifndef XTRACK_MULTIPOLE_NO_SYNRAD
+        if (radiation_flag){
+            // knl and ksl are scaled by the called functions below using factor_knl_ksl
+            factor_knl_ksl_body *= (1. + delta_taper);
+            factor_knl_ksl_edge *= (1. + delta_taper);
+            // k0, k1, k2, k3, k0s, k1s, k2s, k3s are scaled directly here
+            k0 *= (1 + delta_taper);
+            k1 *= (1 + delta_taper);
+            k2 *= (1 + delta_taper);
+            k3 *= (1 + delta_taper);
+            k0s *= (1 + delta_taper);
+            k1s *= (1 + delta_taper);
+            k2s *= (1 + delta_taper);
+            k3s *= (1 + delta_taper);
+        }
+    #endif
+
+    if (edge_entry_active){
+
+        if (rbend_model == 2){
+            // straight body --> curvature in the edges
+            START_PER_PARTICLE_BLOCK(part0, part);
+                YRotation_single_particle(part, -sin_theta_in, cos_theta_in,
+                                          -sin_theta_in/cos_theta_in);
+                LocalParticle_add_to_x(part, x0_in);
+            END_PER_PARTICLE_BLOCK;
+        }
+
+        double knorm[] = {k0, k1, k2, k3};
+        double kskew[] = {k0s, k1s, k2s, k3s};
+
+        track_magnet_edge_particles(
+            part0,
+            edge_entry_model,
+            0, // is_exit
+            edge_entry_hgap,
+            knorm,
+            kskew,
+            3, // k_order,
+            knl,
+            ksl,
+            factor_knl_ksl_edge,
+            order,
+            ks,
+            x0_solenoid,
+            y0_solenoid,
+            length,
+            edge_entry_angle,
+            edge_entry_angle_fdown,
+            edge_entry_fint,
+            factor_backtrack_edge
+        );
+    }
+
+    if (body_active){
+
+        if (integrator == 0){
+            integrator = default_integrator;
+        }
+        if (model == 0){
+            model = default_model;
+        }
+        if (model==-1){ // kick only
+            integrator = 3; // uniform
+            num_multipole_kicks = 1;
+        }
+
+        // Adjust the number of multipole kicks based on the weight
+        if (weight != 1.0 && num_multipole_kicks > 0){
+            num_multipole_kicks = (int64_t) ceil(num_multipole_kicks * weight);
+        }
+
+        // Compute the number of kicks for auto mode
+        if (num_multipole_kicks == 0) { // num_multipole_kicks = 0 means auto mode
+            // If there are active kicks the number of kicks is guessed. Otherwise,
+            // only the drift is performed.
+            if (!kick_is_inactive(order, knl, ksl, k0, k1, k2, k3, k0s, k1s, k2s, k3s, h)){
+                if (fabs(h) < 1e-8){
+                    num_multipole_kicks = 1; // straight magnet, one multipole kick in the middle
+                }
+                else{
+                    double b_circum = 2 * 3.14159 / fabs(h);
+                    num_multipole_kicks = fabs(core_length) / b_circum / 0.5e-3; // 0.5 mrad per kick (on average)
+                    if (num_multipole_kicks < 1){
+                        num_multipole_kicks = 1;
+                    }
+                }
+            }
+        }
+
+        double k0_drift=0, k1_drift=0, h_drift=0, ks_drift=0;
+        double k0_kick=0, k1_kick=0, h_kick=0;
+        double k0_h_correction=0, k1_h_correction=0;
+        int8_t kick_rot_frame=0;
+        int8_t drift_model=0;
+        configure_tracking_model(
+            model,
+            k0,
+            k1,
+            h,
+            ks,
+            &k0_drift,
+            &k1_drift,
+            &h_drift,
+            &ks_drift,
+            &k0_kick,
+            &k1_kick,
+            &h_kick,
+            &k0_h_correction,
+            &k1_h_correction,
+            &kick_rot_frame,
+            &drift_model
+        );
+
+        double dp_record_exit, dpx_record_exit, dpy_record_exit;
+
+        START_PER_PARTICLE_BLOCK(part0, part);
+            track_magnet_body_single_particle(
+                part, core_length, order, inv_factorial_order,
+                knl, ksl,
+                factor_knl_ksl_body,
+                num_multipole_kicks, kick_rot_frame, drift_model, integrator,
+                k0_drift, k1_drift, ks_drift, h_drift,
+                k0_kick, k1_kick, h_kick, hxl,
+                k0_h_correction, k1_h_correction,
+                k2, k3, k0s, k1s, k2s, k3s,
+                dks_ds,
+                x0_solenoid, y0_solenoid,
+                radiation_flag,
+                1, // spin_flag
+                radiation_record,
+                &dp_record_exit, &dpx_record_exit, &dpy_record_exit
+            );
+        END_PER_PARTICLE_BLOCK;
+
+        if (rbend_model == 2){
+            // straight body --> correct s to match the curved frame
+            double const ds = (core_length_curved - core_length);
+            START_PER_PARTICLE_BLOCK(part0, part);
+                LocalParticle_add_to_s(part, ds);
+                LocalParticle_add_to_zeta(part, ds);
+            END_PER_PARTICLE_BLOCK;
+        }
+    }
+
+    if (edge_exit_active){
+        double knorm[] = {k0, k1, k2, k3};
+        double kskew[] = {k0s, k1s, k2s, k3s};
+
+        track_magnet_edge_particles(
+            part0,
+            edge_exit_model,
+            1, // is_exit
+            edge_exit_hgap,
+            knorm,
+            kskew,
+            3, // k_order,
+            knl,
+            ksl,
+            factor_knl_ksl_edge,
+            order,
+            ks,
+            x0_solenoid,
+            y0_solenoid,
+            length,
+            edge_exit_angle,
+            edge_exit_angle_fdown,
+            edge_exit_fint,
+            factor_backtrack_edge
+        );
+
+        if (rbend_model == 2){
+            // straight body --> curvature in the edges
+            START_PER_PARTICLE_BLOCK(part0, part);
+                LocalParticle_add_to_x(part, -x0_out); // shift by half sagitta
+                YRotation_single_particle(part, -sin_theta_out, cos_theta_out,
+                    -sin_theta_out/cos_theta_out);
+            END_PER_PARTICLE_BLOCK;
+        }
+    }
+
+}
+
+#endif
\ No newline at end of file
diff --git a/xtrack/beam_elements/elements_src/.__afs7D11 b/xtrack/beam_elements/elements_src/.__afs7D11
new file mode 100644
index 000000000..78ec3403a
--- /dev/null
+++ b/xtrack/beam_elements/elements_src/.__afs7D11
@@ -0,0 +1,131 @@
+// copyright ############################### //
+// This file is part of the Xtrack Package.  //
+// Copyright (c) CERN, 2024.                 //
+// ######################################### //
+#ifndef XTRACK_TRACK_MULT_FRINGE_H
+#define XTRACK_TRACK_MULT_FRINGE_H
+
+#include <headers/track.h>
+
+// This functionality is ported from MAD-NG
+
+GPUFUN
+void MultFringe_track_single_particle(
+    LocalParticle* part,  // Particle to be tracked
+    const double* kn,  // Normal components; array of length `order`
+    const double* ks,  // Skew components; array of length `order`
+    int64_t k_order,  // Order components
+    const double* knl,  // Second set of normal components; array of length kl_order
+    const double* ksl,  // Second set of skey components; array of length kl_order
+    int64_t kl_order,  // Order of the fringe
+    const double length, // Effective length of the magnet corresponding to knl, ksl
+    const uint8_t is_exit,  // If truthy it's the exit fringe, otherwise the entry
+    uint64_t min_order  // Minimum order of the fringe, ignore the lower components
+) {
+    if (k_order == -1 && kl_order == -1) return;
+
+    if (LocalParticle_check_track_flag(part, XS_FLAG_BACKTRACK)) {
+        LocalParticle_kill_particle(part, -32);
+        return;
+    }
+
+    const double beta0 = LocalParticle_get_beta0(part);
+    const double q = LocalParticle_get_q0(part) * LocalParticle_get_charge_ratio(part);
+    const double direction = is_exit ? -1 : 1;
+
+    // Particle coordinates
+    const double x = LocalParticle_get_x(part);
+    const double px = LocalParticle_get_px(part);
+    const double y = LocalParticle_get_y(part);
+    const double py = LocalParticle_get_py(part);
+    const double t = LocalParticle_get_zeta(part) / beta0;
+    const double pt = LocalParticle_get_ptau(part);
+
+    const double rpp = LocalParticle_get_rpp(part);
+
+    double rx = 1;
+    double ix = 0;
+    double fx = 0;
+    double fxx = 0;
+    double fxy = 0;
+    double fy = 0;
+    double fyx = 0;
+    double fyy = 0;
+
+    uint32_t order = (k_order > kl_order) ? k_order : kl_order;
+    double inv_factorial = 1;
+
+    for (uint32_t ii = 0; ii <= order; ii++)
+    {
+        if (ii > 1) inv_factorial /= ii;
+        double component = ii + 1;
+        double drx = rx;
+        double dix = ix;
+        rx = drx * x - dix * y;
+        ix = drx * y + dix * x;
+
+        double kn_total = 0;
+        double ks_total = 0;
+
+        if (ii >= min_order) {
+            if (ii <= k_order) {
+                kn_total += kn[ii] * inv_factorial;
+                ks_total += ks[ii] * inv_factorial;
+            }
+            if (ii <= kl_order && length != 0.) {
+                kn_total += knl[ii] / length * inv_factorial;
+                ks_total += ksl[ii] / length * inv_factorial;
+            }
+        }
+
+        double nj = -q * direction / (4 * (component + 1));
+        double nf = (component + 2) / component;
+        double kj = kn_total;
+        double ksj = ks_total;
+        double u, v, du, dv;
+
+
+        if (ii == 0) {
+            u = nj * (-ksj * ix);
+            v = nj * (ksj * rx);
+            du = nj * (-ksj * dix);
+            dv = nj * (ksj * drx);
+        } else {
+            u = nj * (kj * rx - ksj * ix);
+            v = nj * (kj * ix + ksj * rx);
+            du = nj * (kj * drx - ksj * dix);
+            dv = nj * (kj * dix + ksj * drx);
+        }
+
+        double dux = component * du;
+        double dvx = component * dv;
+        double duy = -component * dv;
+        double dvy = component * du;
+
+        fx = fx + u * x + nf * v * y;
+        fy = fy + u * y - nf * v * x;
+        fxx = fxx + dux * x + nf * dvx * y + u;
+        fyy = fyy + duy * y - nf * dvy * x + u;
+        fxy = fxy + duy * x + nf * (dvy * y + v);
+        fyx = fyx + dux * y - nf * (dvx * x + v);
+
+    }
+
+    double a = 1 - fxx * rpp;
+    double b = -fyx * rpp;
+    double c = -fxy * rpp;
+    double d = 1 - fyy * rpp;
+    double det = (a * d - b * c);
+
+    double new_px = (d * px - b * py) / det;
+    double new_py = (a * py - c * px) / det;
+    double delta_t = (1 / beta0 + pt) * (new_px * fx + new_py * fy) * POW3(rpp);
+
+    LocalParticle_add_to_x(part, -fx * rpp);
+    LocalParticle_add_to_y(part, -fy * rpp);
+    LocalParticle_set_px(part, new_px);
+    LocalParticle_set_py(part, new_py);
+    LocalParticle_set_zeta(part, (t + delta_t) * beta0);
+}
+
+#endif // XTRACK_TRACK_MULT_FRINGE_H
diff --git a/xtrack/beam_elements/elements_src/.__afsB12B b/xtrack/beam_elements/elements_src/.__afsB12B
new file mode 100644
index 000000000..fd745eabc
--- /dev/null
+++ b/xtrack/beam_elements/elements_src/.__afsB12B
@@ -0,0 +1,447 @@
+// copyright ############################### //
+// This file is part of the Xtrack Package.  //
+// Copyright (c) CERN, 2023.                 //
+// ######################################### //
+#ifndef XTRACK_TRACK_MAGNET_DRIFT_H
+#define XTRACK_TRACK_MAGNET_DRIFT_H
+
+#include <headers/track.h>
+
+#define IS_ZERO(X) (fabs(X) < 1e-9)
+
+GPUFUN
+void track_expanded_drift_single_particle(LocalParticle* part, double length){
+    double const rpp    = LocalParticle_get_rpp(part);
+    double const rv0v    = 1./LocalParticle_get_rvv(part);
+    double const xp     = LocalParticle_get_px(part) * rpp;
+    double const yp     = LocalParticle_get_py(part) * rpp;
+    double const dzeta  = 1 - rv0v * ( 1. + ( xp*xp + yp*yp ) / 2. );
+
+    LocalParticle_add_to_x(part, xp * length );
+    LocalParticle_add_to_y(part, yp * length );
+    LocalParticle_add_to_s(part, length);
+    LocalParticle_add_to_zeta(part, length * dzeta );
+}
+
+
+GPUFUN
+void track_exact_drift_single_particle(LocalParticle* part, double length){
+    double const px = LocalParticle_get_px(part);
+    double const py = LocalParticle_get_py(part);
+    double const rv0v    = 1./LocalParticle_get_rvv(part);
+    double const one_plus_delta = 1. + LocalParticle_get_delta(part);
+
+    double const one_over_pz = 1./sqrt(one_plus_delta*one_plus_delta
+                                       - px * px - py * py);
+    double const dzeta = 1 - rv0v * one_plus_delta * one_over_pz;
+
+    LocalParticle_add_to_x(part, px * one_over_pz * length);
+    LocalParticle_add_to_y(part, py * one_over_pz * length);
+    LocalParticle_add_to_zeta(part, dzeta * length);
+    LocalParticle_add_to_s(part, length);
+}
+
+GPUFUN
+void track_polar_drift_single_particle(
+    LocalParticle* part,  // LocalParticle to track
+    const double length,  // length of the element
+    const double h        // curvature
+) {
+
+    // Based on SUBROUTINE Sprotr in PTC and curex_drift in MAD-NG
+
+    const double rvv = LocalParticle_get_rvv(part);
+    // Particle coordinates
+    const double x = LocalParticle_get_x(part);
+    const double y = LocalParticle_get_y(part);
+    const double px = LocalParticle_get_px(part);
+    const double py = LocalParticle_get_py(part);
+    const double s = length;
+
+    const double one_plus_delta = LocalParticle_get_delta(part) + 1.0;
+    const double pz = sqrt(POW2(one_plus_delta) - POW2(px) - POW2(py));
+
+    double new_x, new_px, new_y, delta_ell;
+
+    // Polar drift
+    double const rho = 1 / h;
+    const double ca = cos(h * s);
+    const double sa = sin(h * s);
+    const double sa2 = sin(0.5 * h * s);
+    const double _pz = 1 / pz;
+    const double pxt = px * _pz;
+    const double _ptt = 1 / (ca - sa * pxt);
+    const double pst = (x + rho) * sa * _pz * _ptt;
+
+    new_x = (x + rho * (2 * sa2 * sa2 + sa * pxt)) * _ptt;
+    new_px = ca * px + sa * pz;
+    new_y = y + pst * py;
+    delta_ell = one_plus_delta * (x + rho) * sa / ca / pz / (1 - px * sa / ca / pz);
+
+    // Update Particles object
+    LocalParticle_set_x(part, new_x);
+    LocalParticle_set_px(part, new_px);
+    LocalParticle_set_y(part, new_y);
+    LocalParticle_add_to_zeta(part, length - delta_ell / rvv);
+    LocalParticle_add_to_s(part, s);
+}
+
+
+GPUFUN
+void track_expanded_combined_dipole_quad_single_particle(
+    LocalParticle* part,  // LocalParticle to track
+    const double length,  // length of the element
+    const double k0_,     // normal dipole strength
+    const double k1_,     // normal quadrupole strength
+    const double h        // curvature
+) {
+    // From madx: https://github.com/MethodicalAcceleratorDesign/MAD-X/blob/8695bd422dc403a01aa185e9fea16603bbd5b3e1/src/trrun.f90#L4320
+    // Particle coordinates
+    const double x = LocalParticle_get_x(part);
+    const double y = LocalParticle_get_y(part);
+    const double px = LocalParticle_get_px(part);
+    const double py = LocalParticle_get_py(part);
+    const double rvv = LocalParticle_get_rvv(part);
+
+    // In MAD-X (delta + 1) is computed:
+    // const double delta_plus_1 = sqrt(pt*pt + 2.0*pt*beti + 1.0);
+    const double delta_plus_1 = LocalParticle_get_delta(part) + 1;
+    const double chi = LocalParticle_get_chi(part);
+
+    const double k0 = chi * k0_ / delta_plus_1;
+    const double k1 = chi * k1_ / delta_plus_1;
+
+    const double Kx = k0 * h + k1;
+    const double Ky = -k1;
+
+    double x_, px_, y_, py_, Sx, Sy, Cx, Cy;
+
+    if (Kx > 0.0) {
+        double sqrt_Kx = sqrt(Kx);
+        Sx = sin(sqrt_Kx * length) / sqrt_Kx;
+        Cx = cos(sqrt_Kx * length);
+    }
+    else if (Kx < 0.0) {
+        double sqrt_Kx = sqrt(-Kx); // the imaginary part
+        Sx = sinh(sqrt_Kx * length) / sqrt_Kx; // sin(ix) = i sinh(x)
+        Cx = cosh(sqrt_Kx * length); // cos(ix) = cosh(x)
+    }
+    else { // Kx == 0.0
+        Sx = length;
+        Cx = 1.0;
+    }
+
+    if (Ky > 0.0) {
+        double sqrt_Ky = sqrt(Ky);
+        Sy = sin(sqrt_Ky * length) / sqrt_Ky;
+        Cy = cos(sqrt_Ky * length);
+    }
+    else if (Ky < 0.0) {
+        double sqrt_Ky = sqrt(-Ky); // the imaginary part
+        Sy = sinh(sqrt_Ky * length) / sqrt_Ky; // sin(ix) = i sinh(x)
+        Cy = cosh(sqrt_Ky * length);  // cos(ix) = cosh(x)
+    }
+    else { // Ky == 0.0
+        Sy = length;
+        Cy = 1.0;
+    }
+
+    // useful quantities
+    const double xp = px / delta_plus_1;
+    const double yp = py / delta_plus_1;
+    const double A = -Kx * x - k0 + h;
+    const double B = xp;
+    const double C = -Ky * y;
+    const double D = yp;
+
+    // transverse map
+    x_ = x * Cx + xp * Sx;
+    y_ = y * Cy + yp * Sy;
+    px_ = (A * Sx + B * Cx) * delta_plus_1;
+    py_ = (C * Sy + D * Cy) * delta_plus_1;
+
+    if (NONZERO(Kx))
+        x_ = x_ + (k0 - h) * (Cx - 1.0) / Kx;
+    else
+        x_ = x_ - (k0 - h) * 0.5 * POW2(length);
+
+    // longitudinal map
+    double length_ = length; // will be the total path length traveled by the particle
+    if (NONZERO(Kx)) {
+        length_ -= (h * ((Cx - 1.0) * xp + Sx * A + length * (k0 - h))) / Kx;
+        length_ += 0.5 * (
+            - (POW2(A) * Cx * Sx) / (2.0 * Kx) \
+            + (POW2(B) * Cx * Sx) / 2.0 \
+            + (POW2(A) * length) / (2.0 * Kx) \
+            + (POW2(B) * length) / 2.0 \
+            - (A * B * POW2(Cx)) / Kx \
+            + (A * B) / Kx
+        );
+    }
+    else {
+        length_ += h * length * (
+            3.0 * length * xp \
+            + 6.0 * x \
+            - (k0 - h) * POW2(length)
+        ) / 6.0;
+        length_ += 0.5 * (POW2(B)) * length;
+    }
+
+    if (NONZERO(Ky)) {
+        length_ += 0.5 * (
+            - (POW2(C) * Cy * Sy) / (2.0 * Ky) \
+            + (POW2(D) * Cy * Sy) / 2.0 \
+            + (POW2(C) * length) / (2.0 * Ky) \
+            + (POW2(D) * length) / 2.0 \
+            - (C * D * POW2(Cy)) / Ky \
+            + (C * D) / Ky
+        );
+    }
+    else {
+        length_ += 0.5 * POW2(D) * length;
+    }
+
+    const double dzeta = length - length_ / rvv;
+
+    LocalParticle_set_x(part, x_);
+    LocalParticle_set_px(part, px_);
+    LocalParticle_set_y(part, y_);
+    LocalParticle_set_py(part, py_);
+    LocalParticle_add_to_zeta(part, dzeta);
+    LocalParticle_add_to_s(part, length);
+
+}
+
+GPUFUN
+void track_curved_exact_bend_single_particle(
+    LocalParticle* part,  // LocalParticle to track
+    const double length,  // length of the element
+    const double k0,      // normal dipole strength
+    const double h        // curvature
+) {
+
+    // Here we assume that the caller has ensured h != 0
+
+    double const k0_chi = k0 * LocalParticle_get_chi(part);
+
+    if (fabs(k0_chi) < 1e-8) {
+        track_polar_drift_single_particle(part, length, h);
+        return;
+    }
+
+    const double rvv = LocalParticle_get_rvv(part);
+    // Particle coordinates
+    const double x = LocalParticle_get_x(part);
+    const double y = LocalParticle_get_y(part);
+    const double px = LocalParticle_get_px(part);
+    const double py = LocalParticle_get_py(part);
+    const double s = length;
+
+    const double one_plus_delta = LocalParticle_get_delta(part) + 1.0;
+    const double A = 1.0 / sqrt(POW2(one_plus_delta) - POW2(py));
+    const double pz = sqrt(POW2(one_plus_delta) - POW2(px) - POW2(py));
+
+    double new_x, new_px, new_y, delta_ell;
+
+    // The case for non-zero curvature, s is arc length
+    // Useful constants
+    const double C = pz - k0_chi * ((1 / h) + x);
+    new_px = px * cos(s * h) + C * sin(s * h);
+    double const new_pz = sqrt(POW2(one_plus_delta) - POW2(new_px) - POW2(py));
+    // double const d_new_px_ds = new_px / new_pz;
+
+    const double d_new_px_ds = C * h * cos(h * s) - h * px * sin(h * s);
+
+    // Update particle coordinates
+    new_x = (new_pz * h - d_new_px_ds - k0_chi) / (h * k0_chi);
+    const double D = asin(A * px) - asin(A * new_px);
+    new_y = y + ((py * s) / (k0_chi / h)) + (py / k0_chi) * D;
+
+    delta_ell = ((one_plus_delta * s * h) / k0_chi) + (one_plus_delta / k0_chi) * D;
+
+    // Update Particles object
+    LocalParticle_set_x(part, new_x);
+    LocalParticle_set_px(part, new_px);
+    LocalParticle_set_y(part, new_y);
+    LocalParticle_add_to_zeta(part, length - delta_ell / rvv);
+    LocalParticle_add_to_s(part, s);
+}
+
+GPUFUN
+void track_straight_exact_bend_single_particle(
+    LocalParticle* part,  // LocalParticle to track
+    const double length,  // length of the element
+    const double k0       // normal dipole strength
+) {
+
+    // Here we assume that the caller has ensured h != 0
+
+    double const k0_chi = k0 * LocalParticle_get_chi(part);
+
+    if (fabs(k0_chi) < 1e-8) {
+        track_exact_drift_single_particle(part, length);
+        return;
+    }
+
+    const double rvv = LocalParticle_get_rvv(part);
+    // Particle coordinates
+    const double x = LocalParticle_get_x(part);
+    const double y = LocalParticle_get_y(part);
+    const double px = LocalParticle_get_px(part);
+    const double py = LocalParticle_get_py(part);
+    const double s = length;
+
+    const double one_plus_delta = LocalParticle_get_delta(part) + 1.0;
+    const double A = 1.0 / sqrt(POW2(one_plus_delta) - POW2(py));
+    const double pz = sqrt(POW2(one_plus_delta) - POW2(px) - POW2(py));
+
+    double new_x, new_px, new_y, delta_ell;
+
+    // STRAIGHT EXACT BEND
+    // The case for zero curvature -- straight bend, s is Cartesian length
+    new_px = px - k0_chi * s;
+    new_x = x + (sqrt(POW2(one_plus_delta) - POW2(new_px) - POW2(py)) - pz) / k0_chi;
+
+    const double D = asin(A * px) - asin(A * new_px);
+    new_y = y + (py / k0_chi) * D;
+
+    delta_ell = (one_plus_delta / k0_chi) * D;
+
+    // Update Particles object
+    LocalParticle_set_x(part, new_x);
+    LocalParticle_set_px(part, new_px);
+    LocalParticle_set_y(part, new_y);
+    LocalParticle_add_to_zeta(part, length - delta_ell / rvv);
+    LocalParticle_add_to_s(part, s);
+}
+
+GPUFUN
+void track_solenoid_single_particle(
+    LocalParticle* part,
+    double length,
+    double ks,
+    double x0_solenoid,
+    double y0_solenoid
+) {
+    const double sk = ks / 2;
+
+    if (IS_ZERO(sk)) {
+        track_exact_drift_single_particle(part, length);
+        LocalParticle_set_ax(part, 0);
+        LocalParticle_set_ay(part, 0);
+        return;
+    }
+
+    if (IS_ZERO(length)){
+        return;
+    }
+
+    const double skl = sk * length;
+
+    // Particle coordinates
+    const double x = LocalParticle_get_x(part) - x0_solenoid;
+    const double px = LocalParticle_get_px(part);
+    const double y = LocalParticle_get_y(part) - y0_solenoid;
+    const double py = LocalParticle_get_py(part);
+    const double delta = LocalParticle_get_delta(part);
+    const double rvv = LocalParticle_get_rvv(part);
+
+    // set up constants
+    const double pk1 = px + sk * y;
+    const double pk2 = py - sk * x;
+    const double ptr2 = pk1 * pk1 + pk2 * pk2;
+    const double one_plus_delta = 1 + delta;
+    const double one_plus_delta_sq = one_plus_delta * one_plus_delta;
+    const double pz = sqrt(one_plus_delta_sq - ptr2);
+
+    // set up constants
+    const double cosTh = cos(skl / pz);
+    const double sinTh = sin(skl / pz);
+
+    const double si = sin(skl / pz) / sk;
+    const double rps[4] = {
+        cosTh * x + sinTh * y,
+        cosTh * px + sinTh * py,
+        cosTh * y - sinTh * x,
+        cosTh * py - sinTh * px
+    };
+    double const new_x = cosTh * rps[0] + si * rps[1];
+    double const new_px = cosTh * rps[1] - sk * sinTh * rps[0];
+    double const new_y = cosTh * rps[2] + si * rps[3];
+    double const new_py = cosTh * rps[3] - sk * sinTh * rps[2];
+    double const add_to_zeta = length * (1 - one_plus_delta / (pz * rvv));
+    double const new_ax = -0.5 * ks * new_y;
+    double const new_ay = 0.5 * ks * new_x;
+
+    LocalParticle_set_x(part, new_x + x0_solenoid);
+    LocalParticle_set_px(part, new_px);
+    LocalParticle_set_y(part, new_y + y0_solenoid);
+    LocalParticle_set_py(part, new_py);
+    LocalParticle_add_to_zeta(part, add_to_zeta);
+    LocalParticle_add_to_s(part, length);
+    LocalParticle_set_ax(part, new_ax);
+    LocalParticle_set_ay(part, new_ay);
+
+}
+
+
+
+GPUFUN
+void track_magnet_drift_single_particle(
+    LocalParticle* part,  // LocalParticle to track
+    const double length,  // length of the element
+    const double k0,      // normal dipole strength
+    const double k1,      // normal quadrupole strength
+    const double ks,      // solenoid strength
+    const double h,       // curvature
+    const double x0_solenoid,
+    const double y0_solenoid,
+    const int64_t drift_model      // drift model
+) {
+
+    // drift_model = 0 : drift expanded (caller has ensured k0=0, k1=0, h=0)
+    // drift_model = 1 : drift exact (caller has ensured k0=0, k1=0, h=0)
+    // drift_model = 2 : polar drift (caller has ensured k0=0, k1=0, h!=0)
+    // drift_model = 3 : k0, k1, h expanded map (this is general for all possible values)
+    // drift_model = 4 : bend with h (caller has ensured k1=0, h!=0)
+    // drift_model = 5 : bend without h (caller has ensured k1=0, h=0)
+    // drift_model = -1 : no drift
+
+    if (drift_model == -1) {
+        return;
+    }
+
+    if (length == 0.0) {
+        return;
+    }
+    switch (drift_model) {
+        case 0:
+            track_expanded_drift_single_particle(part, length);
+            break;
+        case 1:
+            track_exact_drift_single_particle(part, length);
+            break;
+        case 2:
+            track_polar_drift_single_particle(part, length, h);
+            break;
+        case 3:
+            track_expanded_combined_dipole_quad_single_particle(part, length, k0, k1, h);
+            break;
+        case 4:
+            track_curved_exact_bend_single_particle(part, length, k0, h);
+            break;
+        case 5:
+            track_straight_exact_bend_single_particle(part, length, k0);
+            break;
+        case 6:
+            track_solenoid_single_particle(part, length, ks, x0_solenoid, y0_solenoid);
+            break;
+        default:
+            break;
+    }
+
+}
+
+#undef IS_ZERO
+
+#endif
diff --git a/xtrack/beam_elements/elements_src/.__afsEE4E b/xtrack/beam_elements/elements_src/.__afsEE4E
new file mode 100644
index 000000000..639e04581
--- /dev/null
+++ b/xtrack/beam_elements/elements_src/.__afsEE4E
@@ -0,0 +1,422 @@
+// copyright ############################### //
+// This file is part of the Xtrack Package.  //
+// Copyright (c) CERN, 2025.                 //
+// ######################################### //
+#ifndef XTRACK_TRACK_MISALIGNMENT_H
+#define XTRACK_TRACK_MISALIGNMENT_H
+
+#include <headers/track.h>
+#include <beam_elements/elements_src/track_xrotation.h>
+#include <beam_elements/elements_src/track_yrotation.h>
+#include <beam_elements/elements_src/track_srotation.h>
+#include <beam_elements/elements_src/track_drift.h>
+#include <beam_elements/elements_src/track_xyshift.h>
+
+
+#define IF_NONZERO(VALUE, EXPRESSION) { if (VALUE != 0.0) { EXPRESSION; } }
+#define LOOP_PARTICLES(PART0, CODE) {\
+    START_PER_PARTICLE_BLOCK(PART0, part) { \
+        CODE ; \
+    } END_PER_PARTICLE_BLOCK; }
+#define Y_ROTATE(PART0, THETA) IF_NONZERO(THETA, LOOP_PARTICLES(PART0, YRotation_single_particle(part, sin(THETA), cos(THETA), tan(THETA))))
+// Flip some signs so that the input matches the MAD-X survey convention
+#define X_ROTATE(PART0, PHI) IF_NONZERO(PHI, LOOP_PARTICLES(PART0, XRotation_single_particle(part, -sin(PHI), cos(PHI), -tan(PHI))))
+#define S_ROTATE(PART0, PSI) IF_NONZERO(PSI, LOOP_PARTICLES(PART0, SRotation_single_particle(part, sin(PSI), cos(PSI))))
+#define XY_SHIFT(PART0, DX, DY) LOOP_PARTICLES(PART0, XYShift_single_particle(part, DX, DY))
+#define S_SHIFT(PART0, DS) IF_NONZERO(DS, LOOP_PARTICLES(PART0, { \
+        Drift_single_particle_exact(part, DS); \
+        LocalParticle_add_to_zeta(part, -DS); \
+        LocalParticle_add_to_s(part, -DS); \
+    }))
+
+
+GPUFUN void matrix_multiply_4x4(const double[4][4], const double[4][4], double[4][4]);
+GPUFUN void matrix_rigid_affine_inverse(const double[4][4], double[4][4]);
+
+
+GPUFUN
+void track_misalignment_entry_straight(
+    LocalParticle* part0,  // LocalParticle to track
+    double dx,  // misalignment in x
+    double dy,  // misalignment in y
+    double ds,  // misalignment in s
+    double theta, // rotation around y, yaw, positive s to x
+    double phi,  // rotation around x, pitch, positive s to y
+    double psi_no_frame,  // rotation around s, roll, positive y to x
+    double anchor, // anchor of the misalignment as offset in m from entry
+    double length,  // length of the misaligned element
+    double psi_with_frame,  // psi_with_frame of the element, positive s to x
+    int8_t backtrack
+) {
+    // Silence the warning about unused variable length
+    (void)length; // kept for API consistency with track_misalignment_exit_straight
+
+    const double mis_x = dx - anchor * cos(phi) * sin(theta);
+    const double mis_y = dy - anchor * sin(phi);
+    const double mis_s = ds - anchor * (cos(phi) * cos(theta) - 1);
+
+    // Apply transformations
+    if (!backtrack){
+        XY_SHIFT(part0, mis_x, mis_y);
+        S_SHIFT(part0, mis_s);
+        Y_ROTATE(part0, theta);
+        X_ROTATE(part0, phi);
+        S_ROTATE(part0, psi_no_frame);
+        S_ROTATE(part0, psi_with_frame);
+    } else {
+        S_ROTATE(part0, -psi_with_frame);
+        S_ROTATE(part0, -psi_no_frame);
+        X_ROTATE(part0, -phi);
+        Y_ROTATE(part0, -theta);
+        S_SHIFT(part0, -mis_s);
+        XY_SHIFT(part0, -mis_x, -mis_y);
+    }
+}
+
+
+GPUFUN
+void track_misalignment_exit_straight(
+    LocalParticle* part0,  // LocalParticle to track
+    double dx,  // misalignment in x
+    double dy,  // misalignment in y
+    double ds,  // misalignment in s
+    double theta, // rotation around y, yaw, positive s to x
+    double phi,  // rotation around x, pitch, positive s to y
+    double psi_no_frame,  // rotation around s, roll, positive y to x
+    double anchor, // anchor of the misalignment as offset in m from entry
+    double length,  // length of the misaligned element
+    double psi_with_frame,  // psi_with_frame of the element, positive s to x
+    int8_t backtrack
+) {
+    const double neg_part_length = anchor - length;
+    const double mis_x = neg_part_length * cos(phi) * sin(theta) - dx;
+    const double mis_y = neg_part_length * sin(phi) - dy;
+    const double mis_s = neg_part_length * (cos(phi) * cos(theta) - 1) - ds;
+
+    // Apply transformations
+    if (!backtrack){
+        S_ROTATE(part0, -psi_with_frame);
+        S_ROTATE(part0, -psi_no_frame);
+        X_ROTATE(part0, -phi);
+        Y_ROTATE(part0, -theta);
+        S_SHIFT(part0, mis_s);
+        XY_SHIFT(part0, mis_x, mis_y);
+    } else {
+        XY_SHIFT(part0, -mis_x, -mis_y);
+        S_SHIFT(part0, -mis_s);
+        Y_ROTATE(part0, theta);
+        X_ROTATE(part0, phi);
+        S_ROTATE(part0, psi_no_frame);
+        S_ROTATE(part0, psi_with_frame);
+    }
+}
+
+
+GPUFUN
+void track_misalignment_entry_curved(
+    LocalParticle* part0,  // LocalParticle to track
+    double dx,  // misalignment in x
+    double dy,  // misalignment in y
+    double ds,  // misalignment in s
+    double theta, // rotation around y, yaw, positive s to x
+    double phi,  // rotation around x, pitch, positive s to y
+    double psi_no_frame,  // rotation around s, roll, positive y to x
+    double anchor, // anchor of the misalignment as offset in m from entry
+    double length,  // length of the misaligned element
+    double angle,  // angle by which the element bends the reference frame
+    double h,  // curvature, only used when length == 0
+    double psi_with_frame,  // psi_with_frame of the element, positive s to x
+    int8_t backtrack
+) {
+    // Handle as a straight case if a thick element has no bending.
+    // For a thin element, we still perform the curved calculation if h is given.
+    if (angle == 0.0 && (length != 0.0 || h == 0.0)) {
+        track_misalignment_entry_straight(part0, dx, dy, ds, theta, phi,
+            psi_no_frame, anchor, length, psi_with_frame, backtrack);
+        return;
+    }
+    // Precompute trigonometric functions
+    const double s_phi = sin(phi), c_phi = cos(phi);
+    const double s_theta = sin(theta), c_theta = cos(theta);
+    const double s_psi = sin(psi_no_frame), c_psi = cos(psi_no_frame);
+
+    /* We need to compute the transformation that takes us from the aligned
+       frame to the entry of the element in the misaligned frame:
+
+       misaligned_entry = matrix_first_part @ misalignment_matrix @ inv(matrix_first_part)
+       
+       where:
+       - misalignment_matrix is the matrix that applies the misalignment
+       - matrix_first_part is the matrix that takes us from the entry of the
+         element to the anchor of the misalignment
+    */
+
+    // Misalignment matrix
+    const double misalignment_matrix[4][4] = {
+        {
+            -s_phi * s_psi * s_theta + c_psi * c_theta,
+            -c_psi * s_phi * s_theta - c_theta * s_psi,
+            c_phi * s_theta,
+            dx
+        },
+        {
+            c_phi * s_psi,
+            c_phi * c_psi,
+            s_phi, dy
+        },
+        {
+            -c_theta * s_phi * s_psi - c_psi * s_theta,
+            -c_psi * c_theta * s_phi + s_psi * s_theta,
+            c_phi * c_theta,
+            ds
+        },
+        {0, 0, 0, 1}
+    };
+
+    // Compute matrix that takes us from the reference point of the misalignment
+    // to the entry of the element
+    if (length != 0.0) h = angle / length;
+    const double part_angle = anchor * h;
+    const double delta_x_first_part = (cos(part_angle) - 1) * cos(psi_with_frame) / h;
+    const double delta_y_first_part = (cos(part_angle) - 1) * sin(psi_with_frame) / h;
+    const double delta_s_first_part = sin(part_angle) / h;
+
+    const double matrix_first_part[4][4] = {
+            {
+                (cos(part_angle) - 1) * POW2(cos(psi_with_frame)) + 1,
+                (cos(part_angle) - 1) * cos(psi_with_frame) * sin(psi_with_frame),
+                -cos(psi_with_frame) * sin(part_angle),
+                delta_x_first_part
+            },
+            {
+                (cos(part_angle) - 1) * cos(psi_with_frame) * sin(psi_with_frame),
+                (cos(part_angle) - 1) * POW2(sin(psi_with_frame)) + 1,
+                -sin(part_angle) * sin(psi_with_frame),
+                delta_y_first_part
+            },
+            {
+                cos(psi_with_frame) * sin(part_angle),
+                sin(part_angle) * sin(psi_with_frame),
+                cos(part_angle),
+                delta_s_first_part
+            },
+            {0, 0, 0, 1}
+        };
+
+    double inv_matrix_first_part[4][4];
+    matrix_rigid_affine_inverse(matrix_first_part, inv_matrix_first_part);
+
+    // Compute the transformation that takes us from the aligned frame to the
+    // entry of the element in the misaligned frame
+    double misaligned_entry[4][4], temp[4][4];
+    matrix_multiply_4x4(matrix_first_part, misalignment_matrix, temp);
+    matrix_multiply_4x4(temp, inv_matrix_first_part, misaligned_entry);
+
+    // Extract the basic transformations from the misalignment matrix
+    const double mis_x = misaligned_entry[0][3];
+    const double mis_y = misaligned_entry[1][3];
+    const double mis_s = misaligned_entry[2][3];
+
+    const double rot_theta = atan2(misaligned_entry[0][2], misaligned_entry[2][2]);
+    const double rot_phi = atan2(misaligned_entry[1][2], sqrt(misaligned_entry[1][0] * misaligned_entry[1][0] + misaligned_entry[1][1] * misaligned_entry[1][1]));
+    const double rot_psi = atan2(misaligned_entry[1][0], misaligned_entry[1][1]);
+
+    // Apply transformations
+    if (!backtrack) {
+        XY_SHIFT(part0, mis_x, mis_y);
+        S_SHIFT(part0, mis_s);
+        Y_ROTATE(part0, rot_theta);
+        X_ROTATE(part0, rot_phi);
+        S_ROTATE(part0, rot_psi);
+        S_ROTATE(part0, psi_with_frame);
+    } else {
+        S_ROTATE(part0, -psi_with_frame);
+        S_ROTATE(part0, -rot_psi);
+        X_ROTATE(part0, -rot_phi);
+        Y_ROTATE(part0, -rot_theta);
+        S_SHIFT(part0, -mis_s);
+        XY_SHIFT(part0, -mis_x, -mis_y);
+    }
+}
+
+
+GPUFUN
+void track_misalignment_exit_curved(
+    LocalParticle* part0,  // LocalParticle to track
+    double dx,  // misalignment in x
+    double dy,  // misalignment in y
+    double ds,  // misalignment in s
+    double theta, // rotation around y, yaw, positive s to x
+    double phi,  // rotation around x, pitch, positive s to y
+    double psi_no_frame,  // rotation around s, roll, positive y to x
+    double anchor, // anchor of the misalignment as a fraction of the length
+    double length,  // length of the misaligned element
+    double angle,  // angle by which the element bends the reference frame
+    double h,  // curvature, only used when length == 0
+    double psi_with_frame,  // psi_with_frame of the element, positive s to x
+    int8_t backtrack  // whether to backtrack the particle
+) {
+    // Handle as a straight case if a thick element has no bending.
+    // For a thin element, we still perform the curved calculation if h is given.
+    if (angle == 0.0 && (length != 0.0 || h == 0.0)) {
+        track_misalignment_exit_straight(
+            part0, dx, dy, ds, theta, phi, psi_no_frame, anchor, length,
+            psi_with_frame, backtrack);
+        return;
+    }
+    // Precompute trigonometric functions
+    double s_phi = sin(phi), c_phi = cos(phi);
+    double s_theta = sin(theta), c_theta = cos(theta);
+    double s_psi = sin(psi_no_frame), c_psi = cos(psi_no_frame);
+
+    /* We need to compute the transformation that takes us from the misaligned
+       frame to the aligned frame:
+
+       realign = inv(matrix_second_part) * inv(misalignment_matrix) * matrix_second_part
+
+       where:
+       - misalignment_matrix is the matrix that applies the misalignment
+       - matrix_second_part is the matrix that takes us from the frame in the
+         middle of the element (anchor) to the end of the element
+    */
+
+    // Misalignment matrix
+    const double misalignment_matrix[4][4] = {
+        {
+            -s_phi * s_psi * s_theta + c_psi * c_theta,
+            -c_psi * s_phi * s_theta - c_theta * s_psi,
+            c_phi * s_theta,
+            dx
+        },
+        {
+            c_phi * s_psi,
+            c_phi * c_psi,
+            s_phi,
+            dy
+        },
+        {
+            -c_theta * s_phi * s_psi - c_psi * s_theta,
+            -c_psi * c_theta * s_phi + s_psi * s_theta,
+            c_phi * c_theta,
+            ds
+        },
+        {0, 0, 0, 1}
+    };
+
+    // Compute the inverse of the misalignment matrix
+    double inv_misalignment_matrix[4][4];
+    matrix_rigid_affine_inverse(misalignment_matrix, inv_misalignment_matrix);
+
+    // Compute the inverse of the matrix that takes us from the point of the
+    // misalignment to the exit of the element.
+    if (length != 0.0) h = angle / length;
+    const double part_angle = angle - h * anchor;
+
+    const double s_part_angle = sin(part_angle), c_part_angle = cos(part_angle);
+    const double s_tilt = sin(psi_with_frame), c_tilt = cos(psi_with_frame);
+
+    const double delta_x_second_part = (c_part_angle - 1) * c_tilt / h;
+    const double delta_y_second_part = (c_part_angle - 1) * s_tilt / h;
+    const double delta_s_second_part = s_part_angle / h;
+
+    const double matrix_second_part[4][4] = {
+        {
+            (c_part_angle - 1) * POW2(c_tilt) + 1,
+            (c_part_angle - 1) * c_tilt * s_tilt,
+            c_tilt * -s_part_angle,
+            delta_x_second_part
+        },
+        {
+            (c_part_angle - 1) * c_tilt * s_tilt,
+            (c_part_angle - 1) * POW2(s_tilt) + 1,
+            -s_part_angle * s_tilt,
+            delta_y_second_part
+        },
+        {
+            c_tilt * s_part_angle,
+            s_part_angle * s_tilt,
+            c_part_angle,
+            delta_s_second_part
+        },
+        {0, 0, 0, 1}
+    };
+
+    double inv_matrix_second_part[4][4];
+    matrix_rigid_affine_inverse(matrix_second_part, inv_matrix_second_part);
+
+    // Compute the realignment matrix
+    double realign[4][4], temp[4][4];
+    matrix_multiply_4x4(inv_matrix_second_part, inv_misalignment_matrix, temp);
+    matrix_multiply_4x4(temp, matrix_second_part, realign);
+
+    // Extract the basic transformations from the realignment matrix
+    double mis_x = realign[0][3];
+    double mis_y = realign[1][3];
+    double mis_s = realign[2][3];
+
+    double rot_theta = atan2(realign[0][2], realign[2][2]);
+    double rot_phi = atan2(realign[1][2], sqrt(realign[1][0] * realign[1][0] + realign[1][1] * realign[1][1]));
+    double rot_psi = atan2(realign[1][0], realign[1][1]);
+
+    // Apply transformations
+    if (!backtrack){
+        S_ROTATE(part0, -psi_with_frame);
+        XY_SHIFT(part0, mis_x, mis_y);
+        S_SHIFT(part0, mis_s);
+        Y_ROTATE(part0, rot_theta);
+        X_ROTATE(part0, rot_phi);
+        S_ROTATE(part0, rot_psi);
+    } else {
+        S_ROTATE(part0, -rot_psi);
+        X_ROTATE(part0, -rot_phi);
+        Y_ROTATE(part0, -rot_theta);
+        S_SHIFT(part0, -mis_s);
+        XY_SHIFT(part0, -mis_x, -mis_y);
+        S_ROTATE(part0, psi_with_frame);
+    }
+}
+
+
+GPUFUN
+void matrix_multiply_4x4(const double a[4][4], const double b[4][4], double result[4][4]) {
+    // Multiply two 4x4 matrices `a` and `b`, and store the result in `result`.
+    for (int i = 0; i < 4; i++) {
+        result[i][0] = a[i][0] * b[0][0] + a[i][1] * b[1][0] + a[i][2] * b[2][0] + a[i][3] * b[3][0];
+        result[i][1] = a[i][0] * b[0][1] + a[i][1] * b[1][1] + a[i][2] * b[2][1] + a[i][3] * b[3][1];
+        result[i][2] = a[i][0] * b[0][2] + a[i][1] * b[1][2] + a[i][2] * b[2][2] + a[i][3] * b[3][2];
+        result[i][3] = a[i][0] * b[0][3] + a[i][1] * b[1][3] + a[i][2] * b[2][3] + a[i][3] * b[3][3];
+    }
+}
+
+
+GPUFUN
+void matrix_rigid_affine_inverse(const double m[4][4], double inv_m[4][4]) {
+    // Compute the inverse `inv_m` of a rigid affine transformation matrix `m`.
+    const double m00 = m[0][0], m01 = m[0][1], m02 = m[0][2], m03 = m[0][3];
+    const double m10 = m[1][0], m11 = m[1][1], m12 = m[1][2], m13 = m[1][3];
+    const double m20 = m[2][0], m21 = m[2][1], m22 = m[2][2], m23 = m[2][3];
+
+    // Invert the rotation part of `m`: as it's rigid, it's simply a transpose
+    inv_m[0][0] = m00;
+    inv_m[0][1] = m10;
+    inv_m[0][2] = m20;
+    inv_m[1][0] = m01;
+    inv_m[1][1] = m11;
+    inv_m[1][2] = m21;
+    inv_m[2][0] = m02;
+    inv_m[2][1] = m12;
+    inv_m[2][2] = m22;
+
+    // Compute the translation part, -m[:3, :3]^{-1} @ m[:3, 3]
+    inv_m[0][3] = -m00 * m03 - m10 * m13 - m20 * m23;
+    inv_m[1][3] = -m01 * m03 - m11 * m13 - m21 * m23;
+    inv_m[2][3] = -m02 * m03 - m12 * m13 - m22 * m23;
+
+    // Fill the last row
+    inv_m[3][0] = 0;
+    inv_m[3][1] = 0;
+    inv_m[3][2] = 0;
+    inv_m[3][3] = 1;
+}
+
+#endif  // XTRACK_TRACK_MISALIGNMENT_H
diff --git a/xtrack/beam_elements/elements_src/.__afsF1AA b/xtrack/beam_elements/elements_src/.__afsF1AA
new file mode 100644
index 000000000..12e929645
--- /dev/null
+++ b/xtrack/beam_elements/elements_src/.__afsF1AA
@@ -0,0 +1,635 @@
+// copyright ############################### //
+// This file is part of the Xtrack Package.  //
+// Copyright (c) CERN, 2021.                 //
+// ######################################### //
+
+#ifndef XTRACK_LINESEGMENTMAP_H
+#define XTRACK_LINESEGMENTMAP_H
+
+#include <headers/track.h>
+
+
+GPUFUN
+void remove_closed_orbit(
+    LocalParticle* part0, double const x_ref, double const px_ref,
+    double const y_ref, double const py_ref){
+
+    START_PER_PARTICLE_BLOCK(part0, part);
+        // Remove closed orbit
+        LocalParticle_add_to_x(part, -x_ref);
+        LocalParticle_add_to_px(part, -px_ref);
+        LocalParticle_add_to_y(part, -y_ref);
+        LocalParticle_add_to_py(part, -py_ref);
+    END_PER_PARTICLE_BLOCK;
+}
+
+GPUFUN
+void add_closed_orbit(
+    LocalParticle* part0, double const x_ref, double const px_ref,
+    double const y_ref, double const py_ref){
+
+    START_PER_PARTICLE_BLOCK(part0, part);
+        // Add closed orbit
+        LocalParticle_add_to_x(part, x_ref);
+        LocalParticle_add_to_px(part, px_ref);
+        LocalParticle_add_to_y(part, y_ref);
+        LocalParticle_add_to_py(part, py_ref);
+    END_PER_PARTICLE_BLOCK;
+}
+
+GPUFUN
+void remove_dispersion(
+    LocalParticle* part0, double const dx_0, double const dpx_0,
+    double const dy_0, double const dpy_0){
+
+    START_PER_PARTICLE_BLOCK(part0, part);
+        // Remove dispersion
+        // Symplecticity correction (not working, to be investigated)
+        // LocalParticle_add_to_zeta(part, (
+        //     dpx_0 * LocalParticle_get_x(part)
+        //     - dx_0 * LocalParticle_get_px(part)
+        //     + dpy_0 * LocalParticle_get_y(part)
+        //     - dy_0 * LocalParticle_get_py(part)
+        //     )/LocalParticle_get_rvv(part));
+        double const delta = LocalParticle_get_delta(part);
+        LocalParticle_add_to_x(part, -dx_0 * delta);
+        LocalParticle_add_to_px(part, -dpx_0 * delta);
+        LocalParticle_add_to_y(part, -dy_0 * delta);
+        LocalParticle_add_to_py(part, -dpy_0 * delta);
+    END_PER_PARTICLE_BLOCK;
+}
+
+GPUFUN
+void add_dispersion(
+    LocalParticle* part0, double const dx_1, double const dpx_1,
+    double const dy_1, double const dpy_1){
+
+    START_PER_PARTICLE_BLOCK(part0, part);
+        // Add dispersion
+        // Symplecticity correction (not working, to be investigated)
+        // LocalParticle_add_to_zeta(part, (
+        //     dpx_1 * LocalParticle_get_x(part)
+        //     - dx_1 * LocalParticle_get_px(part)
+        //     + dpy_1 * LocalParticle_get_y(part)
+        //     - dy_1 * LocalParticle_get_py(part)
+        //     )/LocalParticle_get_rvv(part));
+        double const delta = LocalParticle_get_delta(part);
+        LocalParticle_add_to_x(part, dx_1 * delta);
+        LocalParticle_add_to_px(part, dpx_1 * delta);
+        LocalParticle_add_to_y(part, dy_1 * delta);
+        LocalParticle_add_to_py(part, dpy_1 * delta);
+    END_PER_PARTICLE_BLOCK;
+}
+
+GPUFUN
+void transverse_motion(LocalParticle *part0, LineSegmentMapData el){
+
+    double const qx = LineSegmentMapData_get_qx(el);
+    double const qy = LineSegmentMapData_get_qy(el);
+    double const det_xx = LineSegmentMapData_get_det_xx(el);
+    double const det_xy = LineSegmentMapData_get_det_xy(el);
+    double const det_yx = LineSegmentMapData_get_det_yx(el);
+    double const det_yy = LineSegmentMapData_get_det_yy(el);
+    double const alfx_0 = LineSegmentMapData_get_alfx(el, 0);
+    double const betx_0 = LineSegmentMapData_get_betx(el, 0);
+    double const alfy_0 = LineSegmentMapData_get_alfy(el, 0);
+    double const bety_0 = LineSegmentMapData_get_bety(el, 0);
+    double const alfx_1 = LineSegmentMapData_get_alfx(el, 1);
+    double const betx_1 = LineSegmentMapData_get_betx(el, 1);
+    double const alfy_1 = LineSegmentMapData_get_alfy(el, 1);
+    double const bety_1 = LineSegmentMapData_get_bety(el, 1);
+    int64_t const ndqx = LineSegmentMapData_len_coeffs_dqx(el);
+    int64_t const ndqy = LineSegmentMapData_len_coeffs_dqy(el);
+
+    double const c11_0= LineSegmentMapData_get_c11(el,0);
+    double const c12_0= LineSegmentMapData_get_c12(el,0);
+    double const c21_0= LineSegmentMapData_get_c21(el,0);
+    double const c22_0= LineSegmentMapData_get_c22(el,0);
+
+    double const c11_1= LineSegmentMapData_get_c11(el,1);;
+    double const c12_1= LineSegmentMapData_get_c12(el,1);;
+    double const c21_1= LineSegmentMapData_get_c21(el,1);;
+    double const c22_1= LineSegmentMapData_get_c22(el,1);
+
+    double const coupling_0 = c11_0 + c22_0 + c12_0 + c21_0;
+    double const coupling_1 = c11_1 + c22_1 + c12_1 + c21_1;
+
+    int64_t detuning;
+    double sin_x = 0;
+    double cos_x = 0;
+    double sin_y = 0;
+    double cos_y = 0;
+
+    int64_t any_chroma = 0;
+
+    for (int i_dqx=0; i_dqx<ndqx; i_dqx++){
+        any_chroma += LineSegmentMapData_get_coeffs_dqx(el, 0) != 0;
+    }
+
+    for (int i_dqy=0; i_dqy<ndqy; i_dqy++){
+        any_chroma += LineSegmentMapData_get_coeffs_dqy(el, 0) != 0;
+    }
+
+    if (any_chroma ||
+        det_xx != 0.0 || det_xy != 0.0 || det_yx != 0.0 || det_yy != 0.0){
+        detuning = 1;
+    }
+    else{
+        detuning = 0;
+        sin_x = sin(2 * PI * qx);
+        cos_x = cos(2 * PI * qx);
+        sin_y = sin(2 * PI * qy);
+        cos_y = cos(2 * PI * qy);
+    }
+
+    double const sqrt_betprod_x = sqrt(betx_1 * betx_0);
+    double const sqrt_betprod_y = sqrt(bety_1 * bety_0);
+
+    double const sqrt_betratio_x = sqrt(betx_1 / betx_0);
+    double const sqrt_betratio_y = sqrt(bety_1 / bety_0);
+
+    START_PER_PARTICLE_BLOCK(part0, part);
+        if (detuning){
+            double const J_x = 0.5 * (
+                (1.0 + alfx_0 * alfx_0) / betx_0
+                    * LocalParticle_get_x(part)*LocalParticle_get_x(part)
+                + 2 * alfx_0
+                    * LocalParticle_get_x(part)*LocalParticle_get_px(part)
+                + betx_0
+                    * LocalParticle_get_px(part)*LocalParticle_get_px(part));
+            double const J_y = 0.5 * (
+                (1.0 + alfy_0 * alfy_0) /bety_0
+                    * LocalParticle_get_y(part)*LocalParticle_get_y(part)
+                + 2*alfy_0
+                    * LocalParticle_get_y(part)*LocalParticle_get_py(part)
+                + bety_0
+                    * LocalParticle_get_py(part)*LocalParticle_get_py(part));
+            double phase = 2*PI*(qx + det_xx * J_x + det_xy * J_y);
+            for (int i_dqx=1; i_dqx<ndqx; i_dqx++){
+                phase += 2*PI*(LineSegmentMapData_get_coeffs_dqx(el, i_dqx) *
+                               pow(LocalParticle_get_delta(part), (double)i_dqx));
+            }
+            cos_x = cos(phase);
+            sin_x = sin(phase);
+            phase = 2*PI*(qy + det_yx * J_x + det_yy * J_y);
+            for (int i_dqy=1; i_dqy<ndqy; i_dqy++){
+                phase += 2*PI*(LineSegmentMapData_get_coeffs_dqy(el, i_dqy) *
+                               pow(LocalParticle_get_delta(part), (double)i_dqy));
+            }
+            cos_y = cos(phase);
+            sin_y = sin(phase);
+        }
+
+        double const M00_x = sqrt_betratio_x*(cos_x+alfx_0*sin_x);
+        double const M01_x = sqrt_betprod_x*sin_x;
+        double const M10_x = ((alfx_0-alfx_1)*cos_x
+                    -(1+alfx_0*alfx_1)*sin_x
+                    )/sqrt_betprod_x;
+        double const M11_x = (cos_x-alfx_1*sin_x)/sqrt_betratio_x;
+        double const M00_y = sqrt_betratio_y*(cos_y+alfy_0*sin_y);
+        double const M01_y = sqrt_betprod_y*sin_y;
+        double const M10_y = ((alfy_0-alfy_1)*cos_y
+                    -(1+alfy_0*alfy_1)*sin_y
+                    )/sqrt_betprod_y;
+        double const M11_y = (cos_y-alfy_1*sin_y)/sqrt_betratio_y;
+
+        if (coupling_0 == 0.0 && coupling_1 == 0.0) {
+            // No coupling matrices, use simple transfer matrix
+            double const x_out = M00_x*LocalParticle_get_x(part) + M01_x * LocalParticle_get_px(part);
+            double const px_out = M10_x*LocalParticle_get_x(part) + M11_x * LocalParticle_get_px(part);
+            double const y_out = M00_y*LocalParticle_get_y(part) + M01_y * LocalParticle_get_py(part);
+            double const py_out = M10_y*LocalParticle_get_y(part) + M11_y * LocalParticle_get_py(part);
+
+            LocalParticle_set_x(part, x_out);
+            LocalParticle_set_px(part, px_out);
+            LocalParticle_set_y(part, y_out);
+            LocalParticle_set_py(part, py_out);
+        }
+
+        /* Construct 4x4 symplectic matrix R and its inverse from a 2x2 matrix C */
+        void compute_matrices_V(double const c11, double const c12,
+                          double const c21, double const c22,
+                          double const gamma,
+                          double V[16]) 
+        {
+        /* Initialize R and Rinv to zero */
+            for(int i=0;i<16;i++) {
+            V[i] = 0.0;
+            }
+            /* Identity 2x2 blocks (top-left and bottom-right) */
+            V[0] = gamma;  V[5] = gamma;
+            V[10] = gamma; V[15] = gamma;
+            /* C block (top-right) */
+            V[2] = c11;  V[3] = -c12;
+            V[6] = -c21;  V[7] = c22;
+            /* -C^inv block (bottom-left) */
+            V[8]  = -c22; V[9]  =  -c12;
+            V[12] =  -c21; V[13] = -c11;
+        }
+
+        /* Construct 4x4 symplectic matrix R and its inverse from a 2x2 matrix C */
+        void compute_matrices_Vinv(double const c11, double const c12,
+                          double const c21, double const c22,
+                          double const gamma,
+                          double V_inv[16]) 
+        {
+            /* Initialize R and Rinv to zero */
+            for(int i=0;i<16;i++) {
+                V_inv[i] = 0.0;
+            }
+            /* Identity 2x2 blocks (top-left and bottom-right) */
+            V_inv[0] = gamma;  V_inv[5] = gamma;
+            V_inv[10] = gamma; V_inv[15] = gamma;
+            /* -C block (top-right) */
+            V_inv[2] = -c11;  V_inv[3] =  c12;
+            V_inv[6] =  c21;  V_inv[7] = -c22;
+            /*C^inv block (bottom-left) */
+            V_inv[8]  =  c22; V_inv[9]  =  c12;
+            V_inv[12] =  c21; V_inv[13] =  c11;
+        }
+
+
+        if (coupling_0 != 0.0 || coupling_1 != 0.0) {
+        
+            double const det_0 = c11_0 * c22_0 - c12_0 * c21_0;
+            double const det_1 = c11_1 * c22_1 - c12_1 * c21_1;
+
+            double const gamma_0 = sqrt(1 - det_0);
+            double const gamma_1 = sqrt(1 - det_1);
+
+            double V_1[16], V_0inv[16];
+            double M[16];
+            double T[16];
+            double temp[16];
+
+            compute_matrices_V(c11_1, c12_1, c21_1, c22_1, gamma_1, V_1);  
+            compute_matrices_Vinv(c11_0, c12_0, c21_0, c22_0, gamma_0, V_0inv);
+            
+            //the final transformation is V_1 * M * V_0inv
+
+            for(int i=0;i<16;i++) {
+                M[i] = 0.0;
+            }
+            M[0] = M00_x; M[1] = M01_x; M[2]=0.0; M[3]=0.0;
+            M[4] = M10_x; M[5] = M11_x; M[6]=0.0; M[7]=0.0;
+            M[8] = 0.0; M[9] = 0.0; M[10]=M00_y; M[11]=M01_y;
+            M[12]= 0.0; M[13]= 0.0; M[14]=M10_y; M[15]=M11_y;       
+
+            for(int i=0;i<16;i++) {
+                T[i] = 0.0;
+            }
+
+            for(int i=0;i<16;i++) {
+                temp[i] = 0.0;
+            }
+
+            // temp = M * V_0inv
+            for (int i=0;i<4;i++) {
+                for (int j=0;j<4;j++) {
+                    temp[i*4+j] = 0.0;
+                    for (int k=0;k<4;k++)
+                        temp[i*4+j] += M[i*4+k]*V_0inv[k*4+j];
+                }
+            }
+            // T = V_1 * temp
+            for (int i=0;i<4;i++) {
+                for (int j=0;j<4;j++) {
+                    T[i*4+j] = 0.0;
+                    for (int k=0;k<4;k++)
+                        T[i*4+j] += V_1[i*4+k]*temp[k*4+j];
+                }
+            }
+
+                double const old_x = LocalParticle_get_x(part);
+                double const old_px = LocalParticle_get_px(part);
+                double const old_y = LocalParticle_get_y(part);
+                double const old_py = LocalParticle_get_py(part);
+
+                double const x_out = (
+                    T[0]*old_x + T[1]*old_px + T[2]*old_y + T[3]*old_py);
+                double const px_out = (
+                    T[4]*old_x + T[5]*old_px + T[6]*old_y + T[7]*old_py);
+                double const y_out = (
+                    T[8]*old_x + T[9]*old_px + T[10]*old_y + T[11]*old_py);
+                double const py_out = (
+                    T[12]*old_x + T[13]*old_px + T[14]*old_y + T[15]*old_py);
+
+                LocalParticle_set_x(part, x_out);
+                LocalParticle_set_px(part, px_out);
+                LocalParticle_set_y(part, y_out);
+                LocalParticle_set_py(part, py_out);
+            END_PER_PARTICLE_BLOCK;
+    }
+}
+
+
+GPUFUN
+void longitudinal_motion(LocalParticle *part0,
+                         LineSegmentMapData el){
+
+    int64_t const mode_flag = LineSegmentMapData_get_longitudinal_mode_flag(el);
+
+    if (mode_flag==1){ // linear motion fixed qs
+        double const qs = LineSegmentMapData_get_qs(el);
+        double const bets = LineSegmentMapData_get_bets(el);
+        double const bucket_length = LineSegmentMapData_get_bucket_length(el)*LocalParticle_get_beta0(part0)*C_LIGHT;
+        double const sin_s = sin(2 * PI * qs);
+        double const cos_s = cos(2 * PI * qs);
+        uint8_t kill_energy_kick = LocalParticle_check_track_flag(
+                    part0, XS_FLAG_KILL_CAVITY_KICK);
+        START_PER_PARTICLE_BLOCK(part0, part);
+            // We set cos_s = 999 if long map is to be skipped
+            double shift = 0.0;
+            if (bucket_length > 0.0) {
+                shift = bucket_length*floor(LocalParticle_get_zeta(part)/bucket_length+0.5);
+                LocalParticle_add_to_zeta(part,-shift);
+            }
+            double const new_zeta = cos_s * LocalParticle_get_zeta(part) - bets * sin_s * LocalParticle_get_pzeta(part);
+            double const new_pzeta = sin_s * LocalParticle_get_zeta(part) / bets + cos_s * LocalParticle_get_pzeta(part);
+
+            LocalParticle_set_zeta(part, new_zeta+shift);
+            if (!kill_energy_kick){
+                LocalParticle_update_pzeta(part, new_pzeta);
+            }
+        END_PER_PARTICLE_BLOCK;
+    }
+    else if (mode_flag==2){ // non-linear motion
+
+        double const alfp =
+            LineSegmentMapData_get_momentum_compaction_factor(el);
+        double const slippage_length =
+            LineSegmentMapData_get_slippage_length(el);
+
+        START_PER_PARTICLE_BLOCK(part0, part);
+            double const gamma0 = LocalParticle_get_gamma0(part);
+            double const eta = alfp - 1.0 / (gamma0 * gamma0);
+            LocalParticle_add_to_zeta(part,
+                -0.5 * eta * slippage_length * LocalParticle_get_delta(part));
+        END_PER_PARTICLE_BLOCK;
+
+        int64_t const n_rf = LineSegmentMapData_len_voltage_rf(el);
+        for (int i_rf=0; i_rf<n_rf; i_rf++){
+
+            double const v_rf = LineSegmentMapData_get_voltage_rf(el,i_rf);
+            double const f_rf = LineSegmentMapData_get_frequency_rf(el,i_rf);
+            double const lag_rf = LineSegmentMapData_get_lag_rf(el,i_rf);
+
+            if (f_rf == 0) continue;
+
+            START_PER_PARTICLE_BLOCK(part0, part);
+                double const K_FACTOR = ( ( double )2.0 *PI ) / C_LIGHT;
+                double const   beta0  = LocalParticle_get_beta0(part);
+                double const   zeta   = LocalParticle_get_zeta(part);
+                double const   q      = fabs(LocalParticle_get_q0(part))
+                                        * LocalParticle_get_charge_ratio(part);
+                double const   tau    = zeta / beta0;
+                double const   phase  = DEG2RAD  * lag_rf - K_FACTOR * f_rf * tau;
+                double const energy   = q * v_rf * sin(phase);
+                LocalParticle_add_to_energy(part, energy, 1);
+            END_PER_PARTICLE_BLOCK;
+        }
+
+        START_PER_PARTICLE_BLOCK(part0, part);
+            double const gamma0 = LocalParticle_get_gamma0(part);
+            double const eta = alfp - 1.0 / (gamma0 * gamma0);
+            LocalParticle_add_to_zeta(part,
+                -0.5 * eta * slippage_length * LocalParticle_get_delta(part));
+        END_PER_PARTICLE_BLOCK;
+    }
+    else if (mode_flag == 3){ // linear motion fixed RF
+
+        double const alfp =
+            LineSegmentMapData_get_momentum_compaction_factor(el);
+        double const slippage_length =
+            LineSegmentMapData_get_slippage_length(el);
+
+        // Assume there is only one RF term (checked in the Python code)
+        double const v_rf = LineSegmentMapData_get_voltage_rf(el,0);
+        double const f_rf = LineSegmentMapData_get_frequency_rf(el,0);
+
+        double const bucket_length = LocalParticle_get_beta0(part0)*C_LIGHT/f_rf;
+
+        START_PER_PARTICLE_BLOCK(part0, part);
+            double const gamma0 = LocalParticle_get_gamma0(part);
+            double const beta0 = LocalParticle_get_beta0(part);
+            double const q0 = LocalParticle_get_q0(part);
+            double const mass0 = LocalParticle_get_mass0(part);
+            double const eta = alfp - 1.0 / (gamma0 * gamma0);
+            double const E0 = mass0 * gamma0;
+
+            double const qs = sqrt(q0 * fabs(eta) * slippage_length * f_rf * v_rf
+                        / (2 * PI * beta0 * beta0 * beta0 * E0 * C_LIGHT));
+            double const bets = eta * slippage_length / (2 * PI * qs);
+
+            double const sin_s = sin(2 * PI * qs);
+            double const cos_s = cos(2 * PI * qs);
+
+            double shift = bucket_length*floor(LocalParticle_get_zeta(part)/bucket_length+0.5);
+            LocalParticle_add_to_zeta(part,-shift);
+
+            double const new_zeta = cos_s * LocalParticle_get_zeta(part) - bets * sin_s * LocalParticle_get_pzeta(part);
+            double const new_pzeta = sin_s * LocalParticle_get_zeta(part) / bets + cos_s * LocalParticle_get_pzeta(part);
+
+            LocalParticle_set_zeta(part, new_zeta+shift);
+            LocalParticle_update_pzeta(part, new_pzeta);
+        END_PER_PARTICLE_BLOCK;
+    }
+}
+
+GPUFUN
+void uncorrelated_radiation_damping(LocalParticle *part0,
+            LineSegmentMapData el){
+
+    START_PER_PARTICLE_BLOCK(part0, part);
+        LocalParticle_set_x(part,
+            LineSegmentMapData_get_damping_factors(el,0,0)*LocalParticle_get_x(part));
+        LocalParticle_set_px(part,
+            LineSegmentMapData_get_damping_factors(el,1,1)*LocalParticle_get_px(part));
+        LocalParticle_set_y(part,
+            LineSegmentMapData_get_damping_factors(el,2,2)*LocalParticle_get_y(part));
+        LocalParticle_set_py(part,
+            LineSegmentMapData_get_damping_factors(el,3,3)*LocalParticle_get_py(part));
+        LocalParticle_set_zeta(part,
+            LineSegmentMapData_get_damping_factors(el,4,4)*LocalParticle_get_zeta(part));
+        LocalParticle_update_pzeta(part,
+            LineSegmentMapData_get_damping_factors(el,5,5)*LocalParticle_get_pzeta(part));
+    END_PER_PARTICLE_BLOCK;
+}
+
+GPUFUN
+void correlated_radiation_damping(LocalParticle *part0,
+            LineSegmentMapData el){
+
+    uint8_t kill_energy_kick = LocalParticle_check_track_flag(
+                    part0, XS_FLAG_KILL_CAVITY_KICK);
+    START_PER_PARTICLE_BLOCK(part0, part);
+        double in[6];
+        in[0] = LocalParticle_get_x(part);
+        in[1] = LocalParticle_get_px(part);
+        in[2] = LocalParticle_get_y(part);
+        in[3] = LocalParticle_get_py(part);
+        in[4] = LocalParticle_get_zeta(part);
+        in[5] = LocalParticle_get_pzeta(part);
+        double out[6];
+        for(unsigned int i=0;i<6;++i){
+            out[i] = 0;
+            for(unsigned int j=0;j<6;++j){
+                out[i] += LineSegmentMapData_get_damping_factors(el,i,j)*in[j];
+            }
+        }
+        LocalParticle_set_x(part, out[0]);
+        LocalParticle_set_px(part, out[1]);
+        LocalParticle_set_y(part, out[2]);
+        LocalParticle_set_py(part, out[3]);
+        LocalParticle_set_zeta(part, out[4]);
+        if (!kill_energy_kick){
+            LocalParticle_update_pzeta(part,out[5]);}
+    END_PER_PARTICLE_BLOCK;
+}
+
+GPUFUN
+void energy_and_reference_increments(LocalParticle *part0,
+    double const energy_increment, double const energy_ref_increment){
+
+    START_PER_PARTICLE_BLOCK(part0, part);
+        // Change energy without change of reference momentum
+        if (energy_increment !=0){
+        LocalParticle_add_to_energy(part, energy_increment, 1);
+        }
+
+        // Change energy reference
+        // In the transverse plane de change is smoothed, i.e.
+        // both the position and the momentum are scaled,
+        // rather than only the momentum.
+        if (energy_ref_increment != 0){
+            double const old_px = LocalParticle_get_px(part);
+            double const old_py = LocalParticle_get_py(part);
+            double const new_energy0 = LocalParticle_get_mass0(part)
+                *LocalParticle_get_gamma0(part) + energy_ref_increment;
+            double const new_p0c = sqrt(new_energy0*new_energy0
+            -LocalParticle_get_mass0(part)*LocalParticle_get_mass0(part));
+            double const new_beta0 = new_p0c / new_energy0;
+            double const new_gamma0 = new_energy0 / LocalParticle_get_mass0(part);
+            double const geo_emit_factor = sqrt(LocalParticle_get_beta0(part)
+                    *LocalParticle_get_gamma0(part)/new_beta0/new_gamma0);
+            LocalParticle_update_p0c(part, new_p0c); // updates px, py but not in the smoothed way
+            LocalParticle_set_px(part, old_px * geo_emit_factor);
+            LocalParticle_set_py(part, old_py * geo_emit_factor);
+            LocalParticle_scale_x(part, geo_emit_factor);
+            LocalParticle_scale_y(part, geo_emit_factor);
+        }
+    END_PER_PARTICLE_BLOCK;
+}
+
+GPUFUN
+void uncorrelated_gaussian_noise(LocalParticle *part0,
+                    LineSegmentMapData el){
+
+        START_PER_PARTICLE_BLOCK(part0, part);
+            double r = RandomNormal_generate(part);
+            LocalParticle_add_to_x(part,
+                r*LineSegmentMapData_get_gauss_noise_matrix(el,0,0));
+            r = RandomNormal_generate(part);
+            LocalParticle_add_to_px(part,
+                r*LineSegmentMapData_get_gauss_noise_matrix(el,1,1));
+            r = RandomNormal_generate(part);
+            LocalParticle_add_to_y(part,
+                r*LineSegmentMapData_get_gauss_noise_matrix(el,2,2));
+            r = RandomNormal_generate(part);
+            LocalParticle_add_to_py(part,
+                r*LineSegmentMapData_get_gauss_noise_matrix(el,3,3));
+            r = RandomNormal_generate(part);
+            LocalParticle_add_to_zeta(part,
+                r*LineSegmentMapData_get_gauss_noise_matrix(el,4,4));
+            r = RandomNormal_generate(part);
+            double pzeta = LocalParticle_get_pzeta(part);
+            pzeta += r*LineSegmentMapData_get_gauss_noise_matrix(el,5,5);
+            LocalParticle_update_pzeta(part,pzeta);
+        END_PER_PARTICLE_BLOCK;
+}
+
+GPUFUN
+void correlated_gaussian_noise(LocalParticle *part0,
+                    LineSegmentMapData el){
+
+    START_PER_PARTICLE_BLOCK(part0, part);
+        double x[6];
+        double r[6];
+        for(unsigned int i=0;i<6;++i){
+            x[i] = RandomNormal_generate(part);
+            r[i] = 0.0;
+        }
+        for(unsigned int i=0;i<6;++i){
+            for(unsigned int j=0;j<6;++j){
+                r[i] += LineSegmentMapData_get_gauss_noise_matrix(el,i,j)*x[j];
+            }
+        }
+        LocalParticle_add_to_x(part,r[0]);
+        LocalParticle_add_to_px(part,r[1]);
+        LocalParticle_add_to_y(part,r[2]);
+        LocalParticle_add_to_py(part,r[3]);
+        LocalParticle_add_to_zeta(part,r[4]);
+        double pzeta = LocalParticle_get_pzeta(part);
+        LocalParticle_update_pzeta(part,pzeta+r[5]);
+    END_PER_PARTICLE_BLOCK;
+}
+
+GPUFUN
+void LineSegmentMap_track_local_particle(LineSegmentMapData el, LocalParticle* part0){
+
+    remove_closed_orbit(part0,
+        LineSegmentMapData_get_x_ref(el, 0),
+        LineSegmentMapData_get_px_ref(el, 0),
+        LineSegmentMapData_get_y_ref(el, 0),
+        LineSegmentMapData_get_py_ref(el, 0));
+
+    remove_dispersion(part0,
+        LineSegmentMapData_get_dx(el, 0),
+        LineSegmentMapData_get_dpx(el, 0),
+        LineSegmentMapData_get_dy(el, 0),
+        LineSegmentMapData_get_dpy(el, 0));
+
+    transverse_motion(part0, el);
+
+    longitudinal_motion(part0, el);
+
+    energy_and_reference_increments(part0,
+        LineSegmentMapData_get_energy_increment(el),
+        LineSegmentMapData_get_energy_ref_increment(el));
+
+    if (LineSegmentMapData_get_uncorrelated_rad_damping(el) == 1){
+        uncorrelated_radiation_damping(part0,el);
+    }
+    
+    if (LineSegmentMapData_get_correlated_rad_damping(el) == 1){
+        correlated_radiation_damping(part0,el);
+    }
+
+    if (LineSegmentMapData_get_uncorrelated_gauss_noise(el) == 1){
+        uncorrelated_gaussian_noise(part0,el);
+    }
+    
+    if (LineSegmentMapData_get_correlated_gauss_noise(el) == 1){
+        correlated_gaussian_noise(part0,el);
+    }
+
+
+    add_dispersion(part0,
+        LineSegmentMapData_get_dx(el, 1),
+        LineSegmentMapData_get_dpx(el, 1),
+        LineSegmentMapData_get_dy(el, 1),
+        LineSegmentMapData_get_dpy(el, 1));
+
+    add_closed_orbit(part0,
+        LineSegmentMapData_get_x_ref(el, 1),
+        LineSegmentMapData_get_px_ref(el, 1),
+        LineSegmentMapData_get_y_ref(el, 1),
+        LineSegmentMapData_get_py_ref(el, 1));
+
+    double const length = LineSegmentMapData_get_length(el);
+    START_PER_PARTICLE_BLOCK(part0, part);
+        LocalParticle_add_to_s(part, length);
+    END_PER_PARTICLE_BLOCK;
+}
+
+#endif
diff --git a/xtrack/beam_elements/elements_src/linesegmentmap.h b/xtrack/beam_elements/elements_src/linesegmentmap.h
index 9abbb421e..8dc2d8e9c 100644
--- a/xtrack/beam_elements/elements_src/linesegmentmap.h
+++ b/xtrack/beam_elements/elements_src/linesegmentmap.h
@@ -220,11 +220,11 @@ void transverse_motion(LocalParticle *part0, LineSegmentMapData el){
             V[0] = gamma;  V[5] = gamma;
             V[10] = gamma; V[15] = gamma;
             /* C block (top-right) */
-            V[2] = c11;  V[3] = -c12;
-            V[6] = -c21;  V[7] = c22;
+            V[2] = c11;  V[3] = c12;
+            V[6] = c21;  V[7] = c22;
             /* -C^inv block (bottom-left) */
-            V[8]  = -c22; V[9]  =  -c12;
-            V[12] =  -c21; V[13] = -c11;
+            V[8]  = -c22; V[9]  =  c12;
+            V[12] =  c21; V[13] = -c11;
         }
 
         /* Construct 4x4 symplectic matrix R and its inverse from a 2x2 matrix C */
@@ -241,11 +241,11 @@ void transverse_motion(LocalParticle *part0, LineSegmentMapData el){
             V_inv[0] = gamma;  V_inv[5] = gamma;
             V_inv[10] = gamma; V_inv[15] = gamma;
             /* -C block (top-right) */
-            V_inv[2] = -c11;  V_inv[3] =  c12;
-            V_inv[6] =  c21;  V_inv[7] = -c22;
+            V_inv[2] = -c11;  V_inv[3] =  -c12;
+            V_inv[6] =  -c21;  V_inv[7] = -c22;
             /*C^inv block (bottom-left) */
-            V_inv[8]  =  c22; V_inv[9]  =  c12;
-            V_inv[12] =  c21; V_inv[13] =  c11;
+            V_inv[8]  =  c22; V_inv[9]  =  -c12;
+            V_inv[12] =  -c21; V_inv[13] =  c11;
         }
 
 
diff --git a/xtrack/particles/rng_src/.__afsCCC4 b/xtrack/particles/rng_src/.__afsCCC4
new file mode 100644
index 000000000..86b34e7fb
--- /dev/null
+++ b/xtrack/particles/rng_src/.__afsCCC4
@@ -0,0 +1,64 @@
+// copyright ############################### //
+// This file is part of the Xtrack Package.  //
+// Copyright (c) CERN, 2024.                 //
+// ######################################### //
+
+#ifndef XTRACK_BASE_RNG_H
+#define XTRACK_BASE_RNG_H
+
+#include <headers/track.h>
+
+// Only include if compiled standalone, outside of Xsuite
+#ifndef START_PER_PARTICLE_BLOCK
+#include <stdint.h>
+#endif
+
+// Combined LCG-Thausworthe generator from (example 37-4):
+// https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-37-efficient-random-number-generation-and-application
+#define MASK 0xffffffffUL
+#define TAUSWORTHE(s,a,b,c,d) ((((s) &c) <<d) &MASK) ^ (((((s) <<a) &MASK)^(s)) >>b)
+#define LCG(s,A,C) ((((A*(s)) &MASK) + C) &MASK)
+
+GPUFUN
+uint32_t rng_get_int32 (uint32_t *s1, uint32_t *s2, uint32_t *s3, uint32_t *s4 ){
+  *s1 = TAUSWORTHE (*s1, 13, 19, 4294967294UL, 12);  // p1=2^31-1
+  *s2 = TAUSWORTHE (*s2, 2, 25, 4294967288UL, 4);    // p2=2^30-1
+  *s3 = TAUSWORTHE (*s3, 3, 11, 4294967280UL, 17);   // p3=2^28-1
+  *s4 = LCG(*s4, 1664525, 1013904223UL);             // p4=2^32
+
+  // Combined period is lcm(p1,p2,p3,p4) ~ 2^121
+  return ((*s1) ^ (*s2) ^ (*s3) ^ (*s4));
+}
+
+#ifndef TWO_TO_32
+#define TWO_TO_32 4294967296.0
+#endif
+
+GPUFUN
+double rng_get (uint32_t *s1, uint32_t *s2, uint32_t *s3, uint32_t *s4 ){
+
+  return rng_get_int32(s1, s2, s3, s4) / TWO_TO_32; // uniform in [0, 1) 1e10 resolution
+
+}
+
+GPUFUN
+void rng_set (uint32_t *s1, uint32_t *s2, uint32_t *s3, uint32_t *s4, uint32_t s ){
+  *s1 = LCG (s, 69069, 0);
+  if (*s1 < 2) *s1 += 2UL;
+  *s2 = LCG (*s1, 69069, 0);
+  if (*s2 < 8) *s2 += 8UL;
+  *s3 = LCG (*s2, 69069, 0);
+  if (*s3 < 16) *s3 += 16UL;
+  *s4 = LCG (*s3, 69069, 0);
+
+  /* "warm it up" */
+  rng_get (s1, s2, s3, s4);
+  rng_get (s1, s2, s3, s4);
+  rng_get (s1, s2, s3, s4);
+  rng_get (s1, s2, s3, s4);
+  rng_get (s1, s2, s3, s4);
+  rng_get (s1, s2, s3, s4);
+  return;
+}
+
+#endif /* XTRACK_BASE_RNG_H */

From f14a025fd65a925f94a05cd9689b9037fa35cae1 Mon Sep 17 00:00:00 2001
From: NetchAna <vaibhavi.gawas@cern.ch>
Date: Wed, 25 Mar 2026 14:54:54 +0100
Subject: [PATCH 3/4] Fix coupling matrix sign + add test for coupled
 LineSegmentMap

---
 tests/test_coupled_line_segment_map.py       | 203 ++++++
 xtrack/beam_elements/elements_src/.__afs0C09 | 160 -----
 xtrack/beam_elements/elements_src/.__afs1C1F | 323 ----------
 xtrack/beam_elements/elements_src/.__afs2C4D | 157 -----
 xtrack/beam_elements/elements_src/.__afs5693 | 624 ------------------
 xtrack/beam_elements/elements_src/.__afs7D11 | 131 ----
 xtrack/beam_elements/elements_src/.__afsB12B | 447 -------------
 xtrack/beam_elements/elements_src/.__afsEE4E | 422 ------------
 xtrack/beam_elements/elements_src/.__afsF1AA | 635 -------------------
 xtrack/particles/rng_src/.__afsCCC4          |  64 --
 10 files changed, 203 insertions(+), 2963 deletions(-)
 create mode 100644 tests/test_coupled_line_segment_map.py
 delete mode 100644 xtrack/beam_elements/elements_src/.__afs0C09
 delete mode 100644 xtrack/beam_elements/elements_src/.__afs1C1F
 delete mode 100644 xtrack/beam_elements/elements_src/.__afs2C4D
 delete mode 100644 xtrack/beam_elements/elements_src/.__afs5693
 delete mode 100644 xtrack/beam_elements/elements_src/.__afs7D11
 delete mode 100644 xtrack/beam_elements/elements_src/.__afsB12B
 delete mode 100644 xtrack/beam_elements/elements_src/.__afsEE4E
 delete mode 100644 xtrack/beam_elements/elements_src/.__afsF1AA
 delete mode 100644 xtrack/particles/rng_src/.__afsCCC4

diff --git a/tests/test_coupled_line_segment_map.py b/tests/test_coupled_line_segment_map.py
new file mode 100644
index 000000000..e00d2137e
--- /dev/null
+++ b/tests/test_coupled_line_segment_map.py
@@ -0,0 +1,203 @@
+import numpy as np
+import xobjects as xo
+import xtrack as xt
+import xpart as xp
+
+
+# ============================================================
+# RANDOM TEST SETTINGS
+# ============================================================
+
+N_TESTS = 5
+RDT_MAX = 1e-2
+
+
+# ============================================================
+# FCC PARAMETERS (FROM YOU)
+# ============================================================
+
+energy = 120
+p0c = 120e9
+n_ips = 4
+
+beta_x = 240e-3
+beta_y = 1.0e-3
+
+alpha_x = 0.0
+alpha_y = 0.0
+
+Qx = 398.150 / n_ips
+Qy = 398.220 / n_ips
+Qs = 0.0334 / n_ips
+
+phi = 15e-3
+k2_factor = 0.50
+
+sigma_delta_bs = 0.00176
+sigma_z_bs = 5.59e-3
+beta_s_bs = sigma_z_bs / sigma_delta_bs
+
+
+# sextupole optics (can tune later)
+beta_x_sext = 3.0
+beta_y_sext = 500.0
+
+alpha_x_sext = 0.0
+alpha_y_sext = 0.0
+
+
+context = xo.ContextCpu()
+
+
+# ============================================================
+# RDT → Cbar
+# ============================================================
+
+#D. Sagan and D. Rubin. Linear analysis of coupled lattices.
+#$Phys. Rev. ST Accel. Beams, 2:074001, Jul 1999.
+
+#R. Calaga, R. Tomás, and A. Franchi. Betatron coupling:
+#Merging Hamiltonian and matrix approaches. Phys. Rev. ST
+#Accel. Beams, 8:034001, Mar 2005
+
+#Vaibhavi Gawas, Frank Zimmermann, Rogelio Tomas Garcia,
+#and Xavier Buffat. MODELLING LINEAR COUPLING
+#FOR FCC-ee IN XSUITE . In Proc. IPAC’26.
+
+
+def rdts_to_Cbar(f1001, f1010):
+
+    re1010, im1010 = np.real(f1010), np.imag(f1010)
+    re1001, im1001 = np.real(f1001), np.imag(f1001)
+
+    S11 = (im1010 + im1001)
+    S22 = (im1001 - im1010)
+    S12 = (-re1010 + re1001)
+    S21 = (-re1010 - re1001)
+
+    S = np.array([[S11, S12], [S21, S22]])
+
+    detS = np.linalg.det(S)
+    gamma = 1.0 / np.sqrt(1.0 + 4.0 * detS)
+
+    return 2 * gamma * S
+
+
+def Cbar_to_Cphys(Cbar):
+
+    C11 = Cbar[0, 0] * np.sqrt(beta_x / beta_y)
+    C12 = Cbar[0, 1] * np.sqrt(beta_x * beta_y)
+    C21 = Cbar[1, 0] / np.sqrt(beta_x * beta_y)
+    C22 = Cbar[1, 1] * np.sqrt(beta_y / beta_x)
+
+    return C11, C12, C21, C22
+
+
+# ============================================================
+# TEST LOOP
+# ============================================================
+
+for i in range(N_TESTS):
+
+    # -------- random RDTs ----------
+    r1001 = np.random.uniform(0, RDT_MAX)
+    i1001 = np.random.uniform(0, RDT_MAX)
+    r1010 = np.random.uniform(0, RDT_MAX)
+    i1010 = np.random.uniform(0, RDT_MAX)
+
+    f1001 = r1001 + 1j * i1001
+    f1010 = r1010 + 1j * i1010
+
+    Cbar = rdts_to_Cbar(f1001, f1010)
+    c11, c12, c21, c22 = Cbar_to_Cphys(Cbar)
+
+    # -------- sext strengths ----------
+    k2_left = k2_factor / (2 * phi * beta_y * beta_y_sext) * np.sqrt(beta_x / beta_x_sext)
+    k2_right = k2_left
+
+    # -------- arcs ----------
+    arc_left = xt.LineSegmentMap(
+        _context=context,
+        qx=0.0,
+        qy=0.25,
+        qs=0.0,
+        betx=[beta_x, beta_x_sext],
+        bety=[beta_y, beta_y_sext],
+        alfx=[alpha_x, alpha_x_sext],
+        alfy=[alpha_y, alpha_y_sext],
+        c11=[c11, 0.0],
+        c12=[c12, 0.0],
+        c21=[c21, 0.0],
+        c22=[c22, 0.0],
+        bets=beta_s_bs,
+    )
+
+    arc_mid = xt.LineSegmentMap(
+        _context=context,
+        qx=Qx,
+        qy=Qy - 0.5,
+        qs=Qs,
+        betx=[beta_x_sext, beta_x_sext],
+        bety=[beta_y_sext, beta_y_sext],
+        alfx=[alpha_x_sext, alpha_x_sext],
+        alfy=[alpha_y_sext, alpha_y_sext],
+        bets=beta_s_bs,
+    )
+
+    arc_right = xt.LineSegmentMap(
+        _context=context,
+        qx=0.0,
+        qy=0.25,
+        qs=0.0,
+        betx=[beta_x_sext, beta_x],
+        bety=[beta_y_sext, beta_y],
+        alfx=[alpha_x_sext, alpha_x],
+        alfy=[alpha_y_sext, alpha_y],
+        c11=[0.0, c11],
+        c12=[0.0, c12],
+        c21=[0.0, c21],
+        c22=[0.0, c22],
+        bets=beta_s_bs,
+    )
+
+    sext_left = xt.Multipole(order=2, knl=[0, 0, k2_left], _context=context)
+    sext_right = xt.Multipole(order=2, knl=[0, 0, -k2_right], _context=context)
+
+    line = xt.Line(elements=[
+        arc_mid,
+        sext_right,
+        arc_right,
+        arc_left,
+        sext_left
+    ])
+
+    line.particle_ref = xp.Particles(
+        _context=context,
+        p0c=p0c,
+        mass0=xp.ELECTRON_MASS_EV
+    )
+
+    line.build_tracker(_context = context, use_prebuilt_kernels = False)
+
+    tw = line.twiss(method="4d", coupling_edw_teng=True)
+
+    f1001_tw = tw.f1001[3]
+    f1010_tw = tw.f1010[3]
+
+    print("\n==============================")
+    print("TEST", i + 1)
+
+    print("INPUT")
+    print("f1001 =", f1001)
+    print("f1010 =", f1010)
+
+    print("\nTWISS")
+    print("f1001 =", f1001_tw)
+    print("f1010 =", f1010_tw)
+
+    ## Xsuite twiss.py outputs -np.conj(rdt) of the rdt which should match the input
+    print("\nnegative of Conjugate")
+    print("f1001 =",-np.conj(f1001_tw))
+    print("f1010 =",-np.conj(f1010_tw))
+
+    
\ No newline at end of file
diff --git a/xtrack/beam_elements/elements_src/.__afs0C09 b/xtrack/beam_elements/elements_src/.__afs0C09
deleted file mode 100644
index 26edc6fa1..000000000
--- a/xtrack/beam_elements/elements_src/.__afs0C09
+++ /dev/null
@@ -1,160 +0,0 @@
-// copyright ############################### //
-// This file is part of the Xtrack Package.  //
-// Copyright (c) CERN, 2025.                 //
-// ######################################### //
-#ifndef XTRACK_TRACK_MAGNET_EDGE_H
-#define XTRACK_TRACK_MAGNET_EDGE_H
-
-#include <headers/track.h>
-#include <beam_elements/elements_src/track_dipole_edge_linear.h>
-#include <beam_elements/elements_src/track_yrotation.h>
-#include <beam_elements/elements_src/track_wedge.h>
-#include <beam_elements/elements_src/track_mult_fringe.h>
-#include <beam_elements/elements_src/track_dipole_fringe.h>
-
-
-GPUFUN
-void track_magnet_edge_particles(
-    LocalParticle* part0,
-    const int8_t model,  // 0: linear, 1: full, 2: dipole-only, 3: ax ay cancellation
-    const uint8_t is_exit,
-    const double half_gap,
-    const double* knorm,
-    const double* kskew,
-    const int64_t k_order,
-    const double* knl,
-    const double* ksl,
-    const double factor_knl_ksl,
-    const int64_t kl_order,
-    const double ksol,
-    const double x0_solenoid,
-    const double y0_solenoid,
-    const double length,
-    const double face_angle,
-    const double face_angle_feed_down,
-    const double fringe_integral,
-    const double factor_for_backtrack // -1 for backtracking, 1 for forward tracking
-) {
-    double k0 = 0;
-    if (k_order > -1) k0 += knorm[0];
-    if (fabs(length) > 1e-10 && kl_order > -1) k0 += factor_knl_ksl * knl[0] / length;
-
-    // Assume we are coming from or going to a drift
-    if (is_exit) {
-        START_PER_PARTICLE_BLOCK(part0, part);
-            LocalParticle_set_ax(part, 0.);
-            LocalParticle_set_ay(part, 0.);
-        END_PER_PARTICLE_BLOCK;
-    }
-    else {
-        START_PER_PARTICLE_BLOCK(part0, part);
-            LocalParticle_set_ax(part, -0.5 * ksol * (LocalParticle_get_y(part) - y0_solenoid));
-            LocalParticle_set_ay(part, 0.5 * ksol * (LocalParticle_get_x(part) - x0_solenoid));
-        END_PER_PARTICLE_BLOCK;
-    }
-
-    if (model == 0) {  // Linear model
-        // Calculate coefficients for x and y to compute the px and py kicks
-        double r21, r43;
-        compute_dipole_edge_linear_coefficients(
-            k0, face_angle, face_angle_feed_down, half_gap, fringe_integral,
-            &r21, &r43
-        );
-
-        r21 = r21 * factor_for_backtrack;
-        r43 = r43 * factor_for_backtrack;
-
-        START_PER_PARTICLE_BLOCK(part0, part);
-            DipoleEdgeLinear_single_particle(part, r21, r43);
-        END_PER_PARTICLE_BLOCK;
-        return;
-    }
-    else if (model == 1 || model == 2) { // Full model
-
-        if (factor_for_backtrack < 0) {
-            START_PER_PARTICLE_BLOCK(part0, part);
-                LocalParticle_kill_particle(part, -32);
-            END_PER_PARTICLE_BLOCK;
-        }
-
-        uint8_t should_rotate = 0;
-        double sin_ = 0, cos_ = 1, tan_ = 0;
-        if (fabs(face_angle) > 10e-10) {
-            should_rotate = 1;
-            sin_ = sin(face_angle);
-            cos_ = cos(face_angle);
-            tan_ = tan(face_angle);
-        }
-
-        if (is_exit) k0 = -k0;
-
-        #define MAGNET_Y_ROTATE(PART) \
-            if (should_rotate) YRotation_single_particle((PART), -sin_, cos_, -tan_)
-
-        #define MAGNET_DIPOLE_FRINGE(PART) \
-            DipoleFringe_single_particle((PART), fringe_integral, half_gap, k0)
-
-        #define MAGNET_MULTIPOLE_FRINGE(PART) \
-            MultFringe_track_single_particle( \
-                (PART), \
-                knorm, \
-                kskew, \
-                k_order, \
-                knl, \
-                ksl, \
-                kl_order, \
-                length / factor_knl_ksl, \
-                is_exit, \
-                /* min_order */ 1 \
-            );
-        // Above, I use the length to rescale knl and ksl. Here I am relying on
-        // the fact that the length is only used to obtain kn and ks in
-        // MultFringe_track_single_particle. To be remembered if the fringe
-        // model changes!
-
-        #define MAGNET_WEDGE(PART) \
-            if (should_rotate & (k_order >= 0)) Wedge_single_particle((PART), -face_angle, knorm[0])
-
-        #define MAGNET_QUAD_WEDGE(PART) \
-            if (should_rotate & (k_order >= 1)) Quad_wedge_single_particle((PART), -face_angle, knorm[1])
-
-        if (is_exit == 0){ // entry
-            START_PER_PARTICLE_BLOCK(part0, part);
-                MAGNET_Y_ROTATE(part);
-                MAGNET_DIPOLE_FRINGE(part);
-                if (model == 1){
-                    MAGNET_MULTIPOLE_FRINGE(part);
-                }
-                if (model == 1){
-                    MAGNET_QUAD_WEDGE(part);
-                }
-                MAGNET_WEDGE(part);
-            END_PER_PARTICLE_BLOCK;
-        }
-        else { // exit
-            START_PER_PARTICLE_BLOCK(part0, part);
-                MAGNET_WEDGE(part);
-                if (model == 1){
-                    MAGNET_QUAD_WEDGE(part);
-                }
-                if (model == 1){
-                    MAGNET_MULTIPOLE_FRINGE(part);
-                }
-                MAGNET_DIPOLE_FRINGE(part);
-                MAGNET_Y_ROTATE(part);
-            END_PER_PARTICLE_BLOCK;
-        }
-
-        #undef MAGNET_Y_ROTATE
-        #undef MAGNET_DIPOLE_FRINGE
-        #undef MAGNET_MULTIPOLE_FRINGE
-        #undef MAGNET_WEDGE
-        #undef MAGNET_QUAD_WEDGE
-    }
-    else if (model == 3) { // only ax ay cancellation (already done above)
-        // do nothing
-    }
-    // If model is not 0 or 1, do nothing
-}
-
-#endif // XTRACK_TRACK_MAGNET_EDGE_H
diff --git a/xtrack/beam_elements/elements_src/.__afs1C1F b/xtrack/beam_elements/elements_src/.__afs1C1F
deleted file mode 100644
index c400d0797..000000000
--- a/xtrack/beam_elements/elements_src/.__afs1C1F
+++ /dev/null
@@ -1,323 +0,0 @@
-// copyright ############################### //
-// This file is part of the Xtrack Package.  //
-// Copyright (c) CERN, 2023.                 //
-// ######################################### //
-#ifndef XTRACK_TRACK_MAGNET_KICK_H
-#define XTRACK_TRACK_MAGNET_KICK_H
-
-#include <headers/track.h>
-
-
-GPUFUN
-void kick_simple_single_particle(
-    LocalParticle* part,
-    int64_t order,
-    double inv_factorial,
-    const double* knl,
-    const double* ksl,
-    double factor,
-    double kick_weight
-);
-
-
-GPUFUN
-void track_magnet_kick_single_particle(
-    LocalParticle* part,
-    double length,
-    int64_t order,
-    double inv_factorial_order,
-    GPUGLMEM const double* knl,
-    GPUGLMEM const double* ksl,
-    double const factor_knl_ksl,
-    double kick_weight,
-    double k0,
-    double k1,
-    double k2,
-    double k3,
-    double k0s,
-    double k1s,
-    double k2s,
-    double k3s,
-    double h,
-    double hxl,
-    double k0_h_correction,
-    double k1_h_correction,
-    uint8_t rot_frame
-){
-
-    double const chi = LocalParticle_get_chi(part);
-    double const x = LocalParticle_get_x(part);
-    double const y = LocalParticle_get_y(part);
-
-    double knl_main[4] = {k0, k1, k2, k3};
-    double ksl_main[4] = {k0s, k1s, k2s, k3s};
-
-    for (int index = 0; index < 4; index++) {
-        knl_main[index] = knl_main[index] * length;
-        ksl_main[index] = ksl_main[index] * length;
-    }
-
-    // multipolar kick
-    kick_simple_single_particle(
-        part,
-        order,
-        inv_factorial_order,
-        knl,
-        ksl,
-        factor_knl_ksl,
-        kick_weight
-    );
-
-    kick_simple_single_particle(
-        part,
-        /* order */ 3,
-        /* inv_factorial_order */ 1. / (3 * 2),
-        knl_main,
-        ksl_main,
-        /* factor_knl_ksl */ 1,
-        kick_weight
-    );
-
-    // Correct for the curvature
-    double dpx = 0;
-    double dpy = 0;
-    double dzeta = 0;
-
-    if (rot_frame) {
-        double const hl = h * length * kick_weight + hxl * kick_weight;
-        dpx += hl * (1. + LocalParticle_get_delta(part));
-        double const rv0v = 1./LocalParticle_get_rvv(part);
-        dzeta += -rv0v * hl * x;
-    }
-
-    double htot = h;
-    if (length != 0) {
-        htot += hxl / length;
-    }
-
-    // Correct for the curvature
-    // k0h correction can be computed from this term in the hamiltonian
-    // H = 1/2 h k0 x^2
-    // (see MAD 8 physics manual, eq. 5.15, and apply Hamilton's eq. dp/ds = -dH/dx)
-    double k0l_mult = 0;
-    if (order >= 0) {
-        k0l_mult = knl[0] * factor_knl_ksl;
-    }
-    dpx += -chi * (k0_h_correction  *length + k0l_mult) * kick_weight * htot * x;
-
-    // k1h correction can be computed from this term in the hamiltonian
-    // H = 1/3 hk1 x^3 - 1/2 hk1 xy^2
-    // (see MAD 8 physics manual, eq. 5.15, and apply Hamilton's eq. dp/ds = -dH/dx)
-    double k1l_mult = 0;
-    if (order >= 1) {
-        k1l_mult = knl[1] * factor_knl_ksl;
-    }
-    dpx += htot * chi * (k1_h_correction * length + k1l_mult) * kick_weight * (-x * x + 0.5 * y * y);
-    dpy += htot * chi * (k1_h_correction * length  + k1l_mult) * kick_weight * x * y;
-
-    LocalParticle_add_to_px(part, dpx);
-    LocalParticle_add_to_py(part, dpy);
-    LocalParticle_add_to_zeta(part, dzeta);
-
-}
-
-
-
-GPUFUN
-uint8_t kick_is_inactive(
-    int64_t order,
-    GPUGLMEM const double* knl,
-    GPUGLMEM const double* ksl,
-    double k0,
-    double k1,
-    double k2,
-    double k3,
-    double k0s,
-    double k1s,
-    double k2s,
-    double k3s,
-    double h
-){
-    if (h != 0) return 0;
-    if (k0 != 0) return 0;
-    if (k1 != 0) return 0;
-    if (k2 != 0) return 0;
-    if (k3 != 0) return 0;
-    if (k0s != 0) return 0;
-    if (k1s != 0) return 0;
-    if (k2s != 0) return 0;
-    if (k3s != 0) return 0;
-
-    for (int index = order; index >= 0; index--) {
-        if (knl[index] != 0) return 0;
-        if (ksl[index] != 0) return 0;
-    }
-
-    return 1;
-
-}
-
-GPUFUN
-void kick_simple_single_coordinates(
-    double const x,
-    double const y,
-    double const chi,
-    int64_t order,
-    double inv_factorial,
-    const double* knl,
-    const double* ksl,
-    double factor,
-    double kick_weight,
-    double *dpx,
-    double *dpy
-) {
-
-    int64_t index = order;
-
-    double dpx_mul = chi * knl[index] * factor * inv_factorial;
-    double dpy_mul = chi * ksl[index] * factor * inv_factorial;
-
-    while( index > 0 )
-    {
-        double const zre = dpx_mul * x - dpy_mul * y;
-        double const zim = dpx_mul * y + dpy_mul * x;
-
-        inv_factorial *= index;
-        index -= 1;
-
-        double this_knl = chi * knl[index] * factor;
-        double this_ksl = chi * ksl[index] * factor;
-
-        dpx_mul = this_knl * inv_factorial + zre;
-        dpy_mul = this_ksl * inv_factorial + zim;
-    }
-
-    dpx_mul = -dpx_mul; // rad
-
-    *dpx = kick_weight * dpx_mul;
-    *dpy = kick_weight * dpy_mul;
-}
-
-
-GPUFUN
-void kick_simple_single_particle(
-    LocalParticle* part,
-    int64_t order,
-    double inv_factorial,
-    const double* knl,
-    const double* ksl,
-    double factor,
-    double kick_weight
-) {
-    double const chi = LocalParticle_get_chi(part);
-    double const x = LocalParticle_get_x(part);
-    double const y = LocalParticle_get_y(part);
-
-    double dpx, dpy;
-
-    kick_simple_single_coordinates(
-        x,
-        y,
-        chi,
-        order,
-        inv_factorial,
-        knl,
-        ksl,
-        factor,
-        kick_weight,
-        &dpx,
-        &dpy);
-
-    LocalParticle_add_to_px(part, dpx);
-    LocalParticle_add_to_py(part, dpy);
-}
-
-GPUFUN
-void evaluate_field_from_strengths(
-    double const p0c,
-    double const q0,
-    double const x,
-    double const y,
-    double length,
-    int64_t order,
-    double inv_factorial_order,
-    GPUGLMEM const double* knl,
-    GPUGLMEM const double* ksl,
-    double const factor_knl_ksl,
-    double k0,
-    double k1,
-    double k2,
-    double k3,
-    double k0s,
-    double k1s,
-    double k2s,
-    double k3s,
-    double ks,
-    double dks_ds,
-    double x0_solenoid,
-    double y0_solenoid,
-    double *Bx_T,
-    double *By_T,
-    double *Bz_T
-){
-    if (length == 0.0) {
-        *Bx_T = 0.0;
-        *By_T = 0.0;
-        *Bz_T = 0.0;
-        return;
-    }
-
-    double knl_main[4] = {k0, k1, k2, k3};
-    double ksl_main[4] = {k0s, k1s, k2s, k3s};
-
-    for (int index = 0; index < 4; index++) {
-        knl_main[index] = knl_main[index] * length;
-        ksl_main[index] = ksl_main[index] * length;
-    }
-
-    // multipolar kick
-    double dpx_mul = 0.;
-    double dpy_mul = 0.;
-    kick_simple_single_coordinates(
-        x,
-        y,
-        1., // chi
-        order,
-        inv_factorial_order,
-        knl,
-        ksl,
-        factor_knl_ksl,
-        1., // kick_weight
-        &dpx_mul,
-        &dpy_mul);
-
-
-    // main kick
-    double dpx_main=0.;
-    double dpy_main=0.;
-    kick_simple_single_coordinates(
-        x,
-        y,
-        1., // chi
-        3, // order
-        1. / (3 * 2), //inv_factorial_order
-        knl_main,
-        ksl_main,
-        1, // factor_knl_ksl,
-        1., // kick_weight
-        &dpx_main,
-        &dpy_main);
-
-    double const dpx = dpx_mul + dpx_main;
-    double const dpy = dpy_mul + dpy_main;
-
-    double const brho_0 = p0c / C_LIGHT / q0; // [T m]
-
-
-    *Bx_T = dpy * brho_0 / length - 0.5 * dks_ds * brho_0 * (x - x0_solenoid); // [T]
-    *By_T = -dpx * brho_0 / length - 0.5 * dks_ds * brho_0 * (y - y0_solenoid); // [T]
-    *Bz_T = ks * brho_0; // [T]
-
-}
-
-#endif
\ No newline at end of file
diff --git a/xtrack/beam_elements/elements_src/.__afs2C4D b/xtrack/beam_elements/elements_src/.__afs2C4D
deleted file mode 100644
index 7a6c2624f..000000000
--- a/xtrack/beam_elements/elements_src/.__afs2C4D
+++ /dev/null
@@ -1,157 +0,0 @@
-#ifndef XTRACK_TRACK_MAGNET_CONFIGURE_H
-#define XTRACK_TRACK_MAGNET_CONFIGURE_H
-
-#define H_TOLERANCE (1e-8)
-
-GPUFUN
-void configure_tracking_model(
-    int64_t model,
-    double k0,
-    double k1,
-    double h,
-    double ks,
-    double* k0_drift,
-    double* k1_drift,
-    double* h_drift,
-    double* ks_drift,
-    double* k0_kick,
-    double* k1_kick,
-    double* h_kick,
-    double* k0_h_correction,
-    double* k1_h_correction,
-    int8_t* kick_rot_frame,
-    int8_t* out_drift_model
-){
-
-    // model = 0 or 1 : adaptive
-    // model = 2: bend-kick-bend
-    // model = 3: rot-kick-rot
-    // model = 4: mat-kick-mat (previously called `expanded`)
-    // model = 5: drift-kick-drift-exact
-    // model = 6: drift-kick-drift-expanded
-    // model = -1: kick only (not exposed in python)
-    // model = -2: sol-kick-sol (not exposed in python)
-
-    if (model==1){
-        model = 3; // backward compatibility
-    }
-
-    int8_t h_is_zero = (fabs(h) < H_TOLERANCE);
-    int64_t drift_model = 0;
-
-    if (model == 2){ // bend-kick-bend
-        if (h_is_zero){
-            drift_model = 5; // bend without h
-        }
-        else{
-            drift_model = 4; // bend with h
-        }
-    }
-    else if(model == 3){ // rot-kick-rot
-        if (h_is_zero){
-            drift_model = 1; // drift exact
-        }
-        else{
-            drift_model = 2; // polar drift
-        }
-    }
-    else if(model == 4){ // mat-kick-mat
-        drift_model = 3; // expanded
-    }
-    else if(model == 5){ // drift-kick-drift-exact
-        drift_model = 1; // drift exact
-    }
-    else if(model == 6){ // drift-kick-drift-expanded
-        drift_model = 0; // drift expanded
-    }
-    else if(model == -1){ // kick only
-        drift_model = -1;
-    }
-    else if(model == -2){ // sol-kick-sol
-        drift_model = 6; // solenoid
-    }
-    else{
-        // This should never happen, but just in case
-        drift_model = 99999999;
-    }
-
-    if (drift_model == -1 || drift_model == 0 || drift_model == 1){ // drift expanded, drift exact, kick only
-        *k0_drift = 0.0;
-        *k1_drift = 0.0;
-        *h_drift = 0.0;
-        *ks_drift = 0.0;
-        *k0_kick = k0;
-        *k1_kick = k1;
-        *h_kick = h;
-        *k0_h_correction = k0;
-        *k1_h_correction = k1;
-        *kick_rot_frame = 1;
-    }
-    else if (drift_model == 2){ // polar drift
-        *k0_drift = 0.0;
-        *k1_drift = 0.0;
-        *h_drift = h;
-        *ks_drift = 0.0;
-        *k0_kick = k0;
-        *k1_kick = k1;
-        *h_kick = h;
-        *k0_h_correction = k0;
-        *k1_h_correction = k1;
-        *kick_rot_frame = 0;
-    }
-    else if (drift_model == 3){ // expanded dipole-quadrupole
-        *k0_drift = k0;
-        *k1_drift = k1;
-        *h_drift = h;
-        *ks_drift = 0.0;
-        *k0_kick = 0.0;
-        *k1_kick = 0.0;
-        *h_kick = h;
-        *k0_h_correction = 0.;
-        *k1_h_correction = k1;
-        *kick_rot_frame = 0;
-    }
-    else if (drift_model == 4){ // bend with h
-        *k0_drift = k0;
-        *k1_drift = 0.0;
-        *h_drift = h;
-        *ks_drift = 0.0;
-        *k0_kick = 0.0;
-        *k1_kick = k1;
-        *h_kick = h;
-        *k0_h_correction = 0.;
-        *k1_h_correction = k1;
-        *kick_rot_frame = 0;
-    }
-    else if (drift_model == 5){ // bend without h
-        *k0_drift = k0;
-        *k1_drift = 0.0;
-        *h_drift = 0.0;
-        *ks_drift = 0.0;
-        *k0_kick = 0.0;
-        *k1_kick = k1;
-        *h_kick = 0.0;
-        *k0_h_correction = 0.;
-        *k1_h_correction = 0.;
-        *kick_rot_frame = 0;
-    }
-    else if (drift_model == 6){ // solenoid
-        *k0_drift = 0.0;
-        *k1_drift = 0.0;
-        *h_drift = 0.0;
-        *ks_drift = ks;
-        *k0_kick = k0;
-        *k1_kick = k1;
-        *h_kick = h;
-        *k0_h_correction = k0;
-        *k1_h_correction = k1;
-        *kick_rot_frame = 1;
-    }
-
-
-    *out_drift_model = drift_model;
-}
-
-#undef H_TOLERANCE
-
-#endif // XTRACK_TRACK_MAGNET_CONFIGURE_H
\ No newline at end of file
diff --git a/xtrack/beam_elements/elements_src/.__afs5693 b/xtrack/beam_elements/elements_src/.__afs5693
deleted file mode 100644
index 299a4799a..000000000
--- a/xtrack/beam_elements/elements_src/.__afs5693
+++ /dev/null
@@ -1,624 +0,0 @@
-// copyright ############################### //
-// This file is part of the Xtrack Package.  //
-// Copyright (c) CERN, 2023.                 //
-// ######################################### //
-#ifndef XTRACK_TRACK_MAGNET_H
-#define XTRACK_TRACK_MAGNET_H
-
-#include <headers/track.h>
-#include <beam_elements/elements_src/track_magnet_kick.h>
-#include <beam_elements/elements_src/track_magnet_drift.h>
-#include <beam_elements/elements_src/track_magnet_edge.h>
-#include <beam_elements/elements_src/track_magnet_configure.h>
-
-#ifndef XTRACK_MULTIPOLE_NO_SYNRAD
-#include <beam_elements/elements_src/track_magnet_radiation.h>
-#endif
-
-#define H_TOLERANCE (1e-8)
-
-#ifndef VSWAP
-    #define VSWAP(a, b) { double tmp = a; a = b; b = tmp; }
-#endif
-
-
-GPUFUN
-void track_magnet_body_single_particle(
-    LocalParticle* part,
-    const double length,
-    const int64_t order,
-    const double inv_factorial_order,
-    GPUGLMEM const double* knl,
-    GPUGLMEM const double* ksl,
-    const double factor_knl_ksl,
-    const int64_t num_multipole_kicks,
-    const int8_t kick_rot_frame,
-    const int8_t drift_model,
-    const int8_t integrator,
-    const double k0_drift,
-    const double k1_drift,
-    const double ks_drift,
-    const double h_drift,
-    const double k0_kick,
-    const double k1_kick,
-    const double h_kick,
-    const double hxl,
-    const double k0_h_correction,
-    const double k1_h_correction,
-    const double k2,
-    const double k3,
-    const double k0s,
-    const double k1s,
-    const double k2s,
-    const double k3s,
-    const double dks_ds,
-    const double x0_solenoid,
-    const double y0_solenoid,
-    const int64_t radiation_flag,
-    const int64_t spin_flag,
-    SynchrotronRadiationRecordData radiation_record,
-    double* dp_record_exit,
-    double* dpx_record_exit,
-    double* dpy_record_exit
-) {
-
-    #define MAGNET_KICK(part, weight) \
-        track_magnet_kick_single_particle(\
-            part, length, order, inv_factorial_order, \
-            knl, ksl, factor_knl_ksl, (weight), \
-            k0_kick, k1_kick, k2, k3, k0s, k1s, k2s, k3s, h_kick,\
-            hxl, k0_h_correction, k1_h_correction, kick_rot_frame\
-        )
-
-    #define MAGNET_DRIFT(part, dlength) \
-        track_magnet_drift_single_particle(\
-            part, (dlength), k0_drift, k1_drift, ks_drift, h_drift,\
-            x0_solenoid, y0_solenoid, drift_model\
-        )
-
-    #ifdef XTRACK_MULTIPOLE_NO_SYNRAD
-        #define WITH_RADIATION(ll, code)\
-        {\
-            code;\
-        }
-    #else
-        #define WITH_RADIATION(ll, code) \
-        { \
-            double const old_x = LocalParticle_get_x(part); \
-            double const old_y = LocalParticle_get_y(part); \
-            double const old_zeta = LocalParticle_get_zeta(part); \
-            double const old_kin_px = LocalParticle_get_px(part) - LocalParticle_get_ax(part);\
-            double const old_kin_py = LocalParticle_get_py(part) - LocalParticle_get_ay(part);\
-            code; \
-            if ((radiation_flag || spin_flag) && length > 0){ \
-                double h_for_rad = h_kick + hxl / length; \
-                if (fabs(h_drift) > 0){ h_for_rad = h_drift; } \
-                double Bx_T, By_T, Bz_T; \
-                double const p0c = LocalParticle_get_p0c(part); \
-                double const q0 = LocalParticle_get_q0(part); \
-                double const new_x = LocalParticle_get_x(part); \
-                double const new_y = LocalParticle_get_y(part); \
-                double const new_kin_px = LocalParticle_get_px(part) - LocalParticle_get_ax(part);\
-                double const new_kin_py = LocalParticle_get_py(part) - LocalParticle_get_ay(part);\
-                double const mean_x = 0.5 * (old_x + new_x); \
-                double const mean_y = 0.5 * (old_y + new_y); \
-                double const mean_kin_px = 0.5 * (old_kin_px + new_kin_px); \
-                double const mean_kin_py = 0.5 * (old_kin_py + new_kin_py); \
-                evaluate_field_from_strengths( \
-                    p0c, \
-                    q0, \
-                    mean_x, \
-                    mean_y, \
-                    length, \
-                    order, \
-                    inv_factorial_order, \
-                    knl, \
-                    ksl, \
-                    factor_knl_ksl, \
-                    k0_drift + k0_kick, \
-                    k1_drift + k1_kick, \
-                    k2, \
-                    k3, \
-                    k0s, \
-                    k1s, \
-                    k2s, \
-                    k3s, \
-                    ks_drift, \
-                    dks_ds, \
-                    x0_solenoid, \
-                    y0_solenoid, \
-                    &Bx_T, \
-                    &By_T, \
-                    &Bz_T \
-                ); \
-                double const dzeta = LocalParticle_get_zeta(part) - old_zeta; \
-                double const rvv = LocalParticle_get_rvv(part); \
-                double l_path = rvv * (ll - dzeta); \
-                if (spin_flag){ \
-                    magnet_spin( \
-                        part, \
-                        Bx_T, \
-                        By_T, \
-                        Bz_T, \
-                        h_for_rad, \
-                        ll, \
-                        l_path); \
-                } \
-                if (radiation_flag){ \
-                    double const B_perp_T = compute_b_perp_mod( \
-                        mean_kin_px, \
-                        mean_kin_py, \
-                        LocalParticle_get_delta(part), \
-                        Bx_T, \
-                        By_T, \
-                        Bz_T); \
-                    magnet_radiation( \
-                        part, \
-                        B_perp_T, \
-                        ll, \
-                        l_path, \
-                        radiation_flag, \
-                        radiation_record, \
-                        dp_record_exit, dpx_record_exit, dpy_record_exit \
-                    ); \
-                } \
-            }\
-        }
-    #endif
-
-    if (num_multipole_kicks == 0 && k0_kick == 0 && k1_kick == 0 && h_kick == 0) { //only drift
-        WITH_RADIATION(length,
-            MAGNET_DRIFT(part, length);
-        )
-    }
-    else{
-
-        
-    // START GENERATED INTEGRATION CODE
-
-    if (integrator == 1){ // TEAPOT
-
-        WITH_RADIATION(length,
-            const double kick_weight = 1. / num_multipole_kicks;
-            double edge_drift_weight = 0.5;
-            double inside_drift_weight = 0;
-            if (num_multipole_kicks > 1) {
-                edge_drift_weight = 1. / (2 * (1 + num_multipole_kicks));
-                inside_drift_weight = (
-                    ((double) num_multipole_kicks)
-                        / ((double)(num_multipole_kicks*num_multipole_kicks) - 1));
-            }
-
-            MAGNET_DRIFT(part, edge_drift_weight*length);
-            for (int i_kick=0; i_kick<num_multipole_kicks - 1; i_kick++) {
-                MAGNET_KICK(part, kick_weight);
-                MAGNET_DRIFT(part, inside_drift_weight*length);
-            }
-            MAGNET_KICK(part, kick_weight);
-            MAGNET_DRIFT(part, edge_drift_weight*length);
-        )
-
-    }
-    else if (integrator==3){ // uniform
-
-        const double kick_weight = 1. / num_multipole_kicks;
-        const double drift_weight = kick_weight;
-
-        for (int i_kick=0; i_kick<num_multipole_kicks; i_kick++) {
-            WITH_RADIATION(drift_weight*length,
-                MAGNET_DRIFT(part, 0.5*drift_weight*length);
-                MAGNET_KICK(part, kick_weight);
-                MAGNET_DRIFT(part, 0.5*drift_weight*length);
-            )
-        }
-
-    }
-    else if (integrator==2){ // YOSHIDA 4
-
-        const int64_t n_kicks_yoshida = 7;
-        const int64_t num_slices = (num_multipole_kicks / n_kicks_yoshida
-                                + (num_multipole_kicks % n_kicks_yoshida != 0));
-
-        const double slice_length = length / (num_slices);
-        const double kick_weight = 1. / num_slices;
-        const double d_yoshida[] =
-                     // From MAD-NG
-                     {3.922568052387799819591407413100e-01,
-                      5.100434119184584780271052295575e-01,
-                      -4.710533854097565531482416645304e-01,
-                      6.875316825251809316199569366290e-02};
-                    //  {0x1.91abc4988937bp-2, 0x1.052468fb75c74p-1, // same in hex
-                    //  -0x1.e25bd194051b9p-2, 0x1.199cec1241558p-4 };
-                    //  {1/8.0, 1/8.0, 1/8.0, 1/8.0}; // Uniform, for debugging
-        const double k_yoshida[] =
-                     // From MAD-NG
-                     {7.845136104775599639182814826199e-01,
-                      2.355732133593569921359289764951e-01,
-                      -1.177679984178870098432412305556e+00,
-                      1.315186320683906284756403692882e+00};
-                    //  {0x1.91abc4988937bp-1, 0x1.e2743579895b4p-3, // same in hex
-                    //  -0x1.2d7c6f7933b93p+0, 0x1.50b00cfb7be3ep+0 };
-                    //  {1/7.0, 1/7.0, 1/7.0, 1/7.0}; // Uniform, for debugging
-
-            for (int ii = 0; ii < num_slices; ii++) {
-                WITH_RADIATION(slice_length,
-                    MAGNET_DRIFT(part, slice_length * d_yoshida[0]);
-                    MAGNET_KICK(part, kick_weight * k_yoshida[0]);
-                    MAGNET_DRIFT(part, slice_length * d_yoshida[1]);
-                    MAGNET_KICK(part, kick_weight * k_yoshida[1]);
-                    MAGNET_DRIFT(part, slice_length * d_yoshida[2]);
-                    MAGNET_KICK(part, kick_weight * k_yoshida[2]);
-                    MAGNET_DRIFT(part, slice_length * d_yoshida[3]);
-                    MAGNET_KICK(part, kick_weight * k_yoshida[3]);
-                    MAGNET_DRIFT(part, slice_length * d_yoshida[3]);
-                    MAGNET_KICK(part, kick_weight * k_yoshida[2]);
-                    MAGNET_DRIFT(part, slice_length * d_yoshida[2]);
-                    MAGNET_KICK(part, kick_weight * k_yoshida[1]);
-                    MAGNET_DRIFT(part, slice_length * d_yoshida[1]);
-                    MAGNET_KICK(part, kick_weight * k_yoshida[0]);
-                    MAGNET_DRIFT(part, slice_length * d_yoshida[0]);
-                ) // WITH_RADIATION
-            }
-    } // integrator if
-
-    // END GENERATED INTEGRATION CODE
-
-
-    }
-
-    #undef MAGNET_KICK
-    #undef MAGNET_DRIFT
-    #undef WITH_RADIATION
-
-}
-
-
-GPUFUN
-void track_magnet_particles(
-    double const weight,
-    LocalParticle* part0,
-    double length,
-    int64_t order,
-    double inv_factorial_order,
-    GPUGLMEM const double* knl,
-    GPUGLMEM const double* ksl,
-    int64_t num_multipole_kicks,
-    int8_t model,
-    int8_t default_model,
-    int8_t integrator,
-    int8_t default_integrator,
-    int64_t radiation_flag,
-    int64_t radiation_flag_parent,
-    SynchrotronRadiationRecordData radiation_record,
-    double delta_taper,
-    double h,
-    double hxl,
-    double k0,
-    double k1,
-    double k2,
-    double k3,
-    double k0s,
-    double k1s,
-    double k2s,
-    double k3s,
-    double ks,
-    double dks_ds,
-    double x0_solenoid,
-    double y0_solenoid,
-    int64_t rbend_model, // -1: not used, 0: auto, 1: curved body, 2: straight body
-    int64_t rbend_compensate_sagitta,
-    double rbend_shift,
-    double rbend_angle_diff,
-    double length_straight,
-    int64_t body_active,
-    int64_t edge_entry_active,
-    int64_t edge_exit_active,
-    int64_t edge_entry_model,
-    int64_t edge_exit_model,
-    double edge_entry_angle,
-    double edge_exit_angle,
-    double edge_entry_angle_fdown,
-    double edge_exit_angle_fdown,
-    double edge_entry_fint,
-    double edge_exit_fint,
-    double edge_entry_hgap,
-    double edge_exit_hgap
-) {
-
-    double factor_knl_ksl = 1.0;
-
-    // Used only for rbend-straight-body
-    double theta_in = 0.;
-    double theta_out = 0.;
-    double cos_theta_in = 1.;
-    double sin_theta_in = 0.;
-    double cos_theta_out = 1.;
-    double sin_theta_out = 0.;
-    double length_curved = 0.;
-    double x0_mid = 0.;
-    double x0_in = 0.;
-    double x0_out = 0.;
-
-    if (rbend_model == 0){
-        // auto mode, curved body
-        rbend_model = 1;
-    }
-
-    double const angle = h * length;
-    if (rbend_model == 1){
-        // curved body
-        edge_entry_angle += (angle - rbend_angle_diff) / 2.0;
-        edge_exit_angle += (angle + rbend_angle_diff) / 2.0;
-    }
-    else if (rbend_model == 2){
-        // straight body
-        theta_in = (angle - rbend_angle_diff) / 2.0;
-        if (fabs(theta_in) > 1e-10){
-            sin_theta_in = sin(theta_in);
-            cos_theta_in = cos(theta_in);
-        }
-        theta_out = (angle + rbend_angle_diff) / 2.0;
-        if (fabs(theta_out) > 1e-10){
-            sin_theta_out = sin(theta_out);
-            cos_theta_out = cos(theta_out);
-        }
-
-        length_curved = length;
-        length = length_straight;
-
-        x0_mid -= rbend_shift;
-
-        if (rbend_compensate_sagitta && fabs(angle) > 1e-10){
-            // shift by half the sagitta
-            double cos_rbha = cos(angle / 2.);
-            x0_mid += 0.5 / h * (1 - cos_rbha);
-        }
-
-        x0_in = x0_mid;
-        x0_out = x0_mid;
-        if (fabs(angle) > 1e-10){
-            double const px0_in = sin(theta_in);
-            double const px0_mid = px0_in - h * length_straight / 2;
-            double const sqrt_mid = sqrt(1 - px0_mid * px0_mid);
-            x0_in -= 1/h *(sqrt_mid - cos_theta_in);
-            x0_out += 1/h * (cos_theta_out - sqrt_mid);
-        }
-        ;
-        h = 0; // treat magnet as straight
-        // We are entering the fringe with an angle, the linear fringe
-        // needs to come from the  expansion around the right angle
-        edge_entry_angle_fdown += theta_in;
-        edge_exit_angle_fdown += theta_out;
-    }
-
-    double core_length, core_length_curved,factor_knl_ksl_body,
-           factor_knl_ksl_edge,factor_backtrack_edge;
-    // Backtracking
-    if (LocalParticle_check_track_flag(part0, XS_FLAG_BACKTRACK)) {
-        core_length = -length * weight;
-        core_length_curved = -length_curved * weight;
-        factor_knl_ksl_body = -factor_knl_ksl * weight;
-        factor_knl_ksl_edge = factor_knl_ksl; // Edge has a specific factor for backtracking
-        factor_backtrack_edge = -1.;
-        hxl = -hxl;
-        VSWAP(edge_entry_active, edge_exit_active);
-        VSWAP(edge_entry_model, edge_exit_model);
-        VSWAP(edge_entry_angle, edge_exit_angle);
-        VSWAP(edge_entry_angle_fdown, edge_exit_angle_fdown);
-        VSWAP(edge_entry_fint, edge_exit_fint);
-        VSWAP(edge_entry_hgap, edge_exit_hgap);
-        VSWAP(theta_in, theta_out);
-        VSWAP(cos_theta_in, cos_theta_out);
-        VSWAP(sin_theta_in, sin_theta_out);
-        theta_in = -theta_in;
-        theta_out = -theta_out;
-        sin_theta_in = -sin_theta_in;
-        sin_theta_out = -sin_theta_out;
-        VSWAP(x0_in, x0_out);
-    }
-    else {
-        core_length = length * weight;
-        core_length_curved = length_curved * weight;
-        factor_knl_ksl_body = factor_knl_ksl * weight;
-        factor_knl_ksl_edge = factor_knl_ksl;
-        factor_backtrack_edge = 1.;
-    }
-
-    #ifndef XTRACK_MULTIPOLE_NO_SYNRAD
-        if (radiation_flag == 10){ // from parent
-            radiation_flag = radiation_flag_parent;
-        }
-    #endif
-
-    // Tapering
-    if (LocalParticle_check_track_flag(part0, XS_FLAG_SR_TAPER)){
-        part0->ipart = 0;
-        delta_taper = LocalParticle_get_delta(part0); // I can use part0 because
-                                                      // there is only one particle
-                                                      // when doing the tapering
-    }
-
-    #ifndef XTRACK_MULTIPOLE_NO_SYNRAD
-        if (radiation_flag){
-            // knl and ksl are scaled by the called functions below using factor_knl_ksl
-            factor_knl_ksl_body *= (1. + delta_taper);
-            factor_knl_ksl_edge *= (1. + delta_taper);
-            // k0, k1, k2, k3, k0s, k1s, k2s, k3s are scaled directly here
-            k0 *= (1 + delta_taper);
-            k1 *= (1 + delta_taper);
-            k2 *= (1 + delta_taper);
-            k3 *= (1 + delta_taper);
-            k0s *= (1 + delta_taper);
-            k1s *= (1 + delta_taper);
-            k2s *= (1 + delta_taper);
-            k3s *= (1 + delta_taper);
-        }
-    #endif
-
-    if (edge_entry_active){
-
-        if (rbend_model == 2){
-            // straight body --> curvature in the edges
-            START_PER_PARTICLE_BLOCK(part0, part);
-                YRotation_single_particle(part, -sin_theta_in, cos_theta_in,
-                                          -sin_theta_in/cos_theta_in);
-                LocalParticle_add_to_x(part, x0_in);
-            END_PER_PARTICLE_BLOCK;
-        }
-
-        double knorm[] = {k0, k1, k2, k3};
-        double kskew[] = {k0s, k1s, k2s, k3s};
-
-        track_magnet_edge_particles(
-            part0,
-            edge_entry_model,
-            0, // is_exit
-            edge_entry_hgap,
-            knorm,
-            kskew,
-            3, // k_order,
-            knl,
-            ksl,
-            factor_knl_ksl_edge,
-            order,
-            ks,
-            x0_solenoid,
-            y0_solenoid,
-            length,
-            edge_entry_angle,
-            edge_entry_angle_fdown,
-            edge_entry_fint,
-            factor_backtrack_edge
-        );
-    }
-
-    if (body_active){
-
-        if (integrator == 0){
-            integrator = default_integrator;
-        }
-        if (model == 0){
-            model = default_model;
-        }
-        if (model==-1){ // kick only
-            integrator = 3; // uniform
-            num_multipole_kicks = 1;
-        }
-
-        // Adjust the number of multipole kicks based on the weight
-        if (weight != 1.0 && num_multipole_kicks > 0){
-            num_multipole_kicks = (int64_t) ceil(num_multipole_kicks * weight);
-        }
-
-        // Compute the number of kicks for auto mode
-        if (num_multipole_kicks == 0) { // num_multipole_kicks = 0 means auto mode
-            // If there are active kicks the number of kicks is guessed. Otherwise,
-            // only the drift is performed.
-            if (!kick_is_inactive(order, knl, ksl, k0, k1, k2, k3, k0s, k1s, k2s, k3s, h)){
-                if (fabs(h) < 1e-8){
-                    num_multipole_kicks = 1; // straight magnet, one multipole kick in the middle
-                }
-                else{
-                    double b_circum = 2 * 3.14159 / fabs(h);
-                    num_multipole_kicks = fabs(core_length) / b_circum / 0.5e-3; // 0.5 mrad per kick (on average)
-                    if (num_multipole_kicks < 1){
-                        num_multipole_kicks = 1;
-                    }
-                }
-            }
-        }
-
-        double k0_drift=0, k1_drift=0, h_drift=0, ks_drift=0;
-        double k0_kick=0, k1_kick=0, h_kick=0;
-        double k0_h_correction=0, k1_h_correction=0;
-        int8_t kick_rot_frame=0;
-        int8_t drift_model=0;
-        configure_tracking_model(
-            model,
-            k0,
-            k1,
-            h,
-            ks,
-            &k0_drift,
-            &k1_drift,
-            &h_drift,
-            &ks_drift,
-            &k0_kick,
-            &k1_kick,
-            &h_kick,
-            &k0_h_correction,
-            &k1_h_correction,
-            &kick_rot_frame,
-            &drift_model
-        );
-
-        double dp_record_exit, dpx_record_exit, dpy_record_exit;
-
-        START_PER_PARTICLE_BLOCK(part0, part);
-            track_magnet_body_single_particle(
-                part, core_length, order, inv_factorial_order,
-                knl, ksl,
-                factor_knl_ksl_body,
-                num_multipole_kicks, kick_rot_frame, drift_model, integrator,
-                k0_drift, k1_drift, ks_drift, h_drift,
-                k0_kick, k1_kick, h_kick, hxl,
-                k0_h_correction, k1_h_correction,
-                k2, k3, k0s, k1s, k2s, k3s,
-                dks_ds,
-                x0_solenoid, y0_solenoid,
-                radiation_flag,
-                1, // spin_flag
-                radiation_record,
-                &dp_record_exit, &dpx_record_exit, &dpy_record_exit
-            );
-        END_PER_PARTICLE_BLOCK;
-
-        if (rbend_model == 2){
-            // straight body --> correct s to match the curved frame
-            double const ds = (core_length_curved - core_length);
-            START_PER_PARTICLE_BLOCK(part0, part);
-                LocalParticle_add_to_s(part, ds);
-                LocalParticle_add_to_zeta(part, ds);
-            END_PER_PARTICLE_BLOCK;
-        }
-    }
-
-    if (edge_exit_active){
-        double knorm[] = {k0, k1, k2, k3};
-        double kskew[] = {k0s, k1s, k2s, k3s};
-
-        track_magnet_edge_particles(
-            part0,
-            edge_exit_model,
-            1, // is_exit
-            edge_exit_hgap,
-            knorm,
-            kskew,
-            3, // k_order,
-            knl,
-            ksl,
-            factor_knl_ksl_edge,
-            order,
-            ks,
-            x0_solenoid,
-            y0_solenoid,
-            length,
-            edge_exit_angle,
-            edge_exit_angle_fdown,
-            edge_exit_fint,
-            factor_backtrack_edge
-        );
-
-        if (rbend_model == 2){
-            // straight body --> curvature in the edges
-            START_PER_PARTICLE_BLOCK(part0, part);
-                LocalParticle_add_to_x(part, -x0_out); // shift by half sagitta
-                YRotation_single_particle(part, -sin_theta_out, cos_theta_out,
-                    -sin_theta_out/cos_theta_out);
-            END_PER_PARTICLE_BLOCK;
-        }
-    }
-
-}
-
-#endif
\ No newline at end of file
diff --git a/xtrack/beam_elements/elements_src/.__afs7D11 b/xtrack/beam_elements/elements_src/.__afs7D11
deleted file mode 100644
index 78ec3403a..000000000
--- a/xtrack/beam_elements/elements_src/.__afs7D11
+++ /dev/null
@@ -1,131 +0,0 @@
-// copyright ############################### //
-// This file is part of the Xtrack Package.  //
-// Copyright (c) CERN, 2024.                 //
-// ######################################### //
-#ifndef XTRACK_TRACK_MULT_FRINGE_H
-#define XTRACK_TRACK_MULT_FRINGE_H
-
-#include <headers/track.h>
-
-// This functionality is ported from MAD-NG
-
-GPUFUN
-void MultFringe_track_single_particle(
-    LocalParticle* part,  // Particle to be tracked
-    const double* kn,  // Normal components; array of length `order`
-    const double* ks,  // Skew components; array of length `order`
-    int64_t k_order,  // Order components
-    const double* knl,  // Second set of normal components; array of length kl_order
-    const double* ksl,  // Second set of skey components; array of length kl_order
-    int64_t kl_order,  // Order of the fringe
-    const double length, // Effective length of the magnet corresponding to knl, ksl
-    const uint8_t is_exit,  // If truthy it's the exit fringe, otherwise the entry
-    uint64_t min_order  // Minimum order of the fringe, ignore the lower components
-) {
-    if (k_order == -1 && kl_order == -1) return;
-
-    if (LocalParticle_check_track_flag(part, XS_FLAG_BACKTRACK)) {
-        LocalParticle_kill_particle(part, -32);
-        return;
-    }
-
-    const double beta0 = LocalParticle_get_beta0(part);
-    const double q = LocalParticle_get_q0(part) * LocalParticle_get_charge_ratio(part);
-    const double direction = is_exit ? -1 : 1;
-
-    // Particle coordinates
-    const double x = LocalParticle_get_x(part);
-    const double px = LocalParticle_get_px(part);
-    const double y = LocalParticle_get_y(part);
-    const double py = LocalParticle_get_py(part);
-    const double t = LocalParticle_get_zeta(part) / beta0;
-    const double pt = LocalParticle_get_ptau(part);
-
-    const double rpp = LocalParticle_get_rpp(part);
-
-    double rx = 1;
-    double ix = 0;
-    double fx = 0;
-    double fxx = 0;
-    double fxy = 0;
-    double fy = 0;
-    double fyx = 0;
-    double fyy = 0;
-
-    uint32_t order = (k_order > kl_order) ? k_order : kl_order;
-    double inv_factorial = 1;
-
-    for (uint32_t ii = 0; ii <= order; ii++)
-    {
-        if (ii > 1) inv_factorial /= ii;
-        double component = ii + 1;
-        double drx = rx;
-        double dix = ix;
-        rx = drx * x - dix * y;
-        ix = drx * y + dix * x;
-
-        double kn_total = 0;
-        double ks_total = 0;
-
-        if (ii >= min_order) {
-            if (ii <= k_order) {
-                kn_total += kn[ii] * inv_factorial;
-                ks_total += ks[ii] * inv_factorial;
-            }
-            if (ii <= kl_order && length != 0.) {
-                kn_total += knl[ii] / length * inv_factorial;
-                ks_total += ksl[ii] / length * inv_factorial;
-            }
-        }
-
-        double nj = -q * direction / (4 * (component + 1));
-        double nf = (component + 2) / component;
-        double kj = kn_total;
-        double ksj = ks_total;
-        double u, v, du, dv;
-
-
-        if (ii == 0) {
-            u = nj * (-ksj * ix);
-            v = nj * (ksj * rx);
-            du = nj * (-ksj * dix);
-            dv = nj * (ksj * drx);
-        } else {
-            u = nj * (kj * rx - ksj * ix);
-            v = nj * (kj * ix + ksj * rx);
-            du = nj * (kj * drx - ksj * dix);
-            dv = nj * (kj * dix + ksj * drx);
-        }
-
-        double dux = component * du;
-        double dvx = component * dv;
-        double duy = -component * dv;
-        double dvy = component * du;
-
-        fx = fx + u * x + nf * v * y;
-        fy = fy + u * y - nf * v * x;
-        fxx = fxx + dux * x + nf * dvx * y + u;
-        fyy = fyy + duy * y - nf * dvy * x + u;
-        fxy = fxy + duy * x + nf * (dvy * y + v);
-        fyx = fyx + dux * y - nf * (dvx * x + v);
-
-    }
-
-    double a = 1 - fxx * rpp;
-    double b = -fyx * rpp;
-    double c = -fxy * rpp;
-    double d = 1 - fyy * rpp;
-    double det = (a * d - b * c);
-
-    double new_px = (d * px - b * py) / det;
-    double new_py = (a * py - c * px) / det;
-    double delta_t = (1 / beta0 + pt) * (new_px * fx + new_py * fy) * POW3(rpp);
-
-    LocalParticle_add_to_x(part, -fx * rpp);
-    LocalParticle_add_to_y(part, -fy * rpp);
-    LocalParticle_set_px(part, new_px);
-    LocalParticle_set_py(part, new_py);
-    LocalParticle_set_zeta(part, (t + delta_t) * beta0);
-}
-
-#endif // XTRACK_TRACK_MULT_FRINGE_H
diff --git a/xtrack/beam_elements/elements_src/.__afsB12B b/xtrack/beam_elements/elements_src/.__afsB12B
deleted file mode 100644
index fd745eabc..000000000
--- a/xtrack/beam_elements/elements_src/.__afsB12B
+++ /dev/null
@@ -1,447 +0,0 @@
-// copyright ############################### //
-// This file is part of the Xtrack Package.  //
-// Copyright (c) CERN, 2023.                 //
-// ######################################### //
-#ifndef XTRACK_TRACK_MAGNET_DRIFT_H
-#define XTRACK_TRACK_MAGNET_DRIFT_H
-
-#include <headers/track.h>
-
-#define IS_ZERO(X) (fabs(X) < 1e-9)
-
-GPUFUN
-void track_expanded_drift_single_particle(LocalParticle* part, double length){
-    double const rpp    = LocalParticle_get_rpp(part);
-    double const rv0v    = 1./LocalParticle_get_rvv(part);
-    double const xp     = LocalParticle_get_px(part) * rpp;
-    double const yp     = LocalParticle_get_py(part) * rpp;
-    double const dzeta  = 1 - rv0v * ( 1. + ( xp*xp + yp*yp ) / 2. );
-
-    LocalParticle_add_to_x(part, xp * length );
-    LocalParticle_add_to_y(part, yp * length );
-    LocalParticle_add_to_s(part, length);
-    LocalParticle_add_to_zeta(part, length * dzeta );
-}
-
-
-GPUFUN
-void track_exact_drift_single_particle(LocalParticle* part, double length){
-    double const px = LocalParticle_get_px(part);
-    double const py = LocalParticle_get_py(part);
-    double const rv0v    = 1./LocalParticle_get_rvv(part);
-    double const one_plus_delta = 1. + LocalParticle_get_delta(part);
-
-    double const one_over_pz = 1./sqrt(one_plus_delta*one_plus_delta
-                                       - px * px - py * py);
-    double const dzeta = 1 - rv0v * one_plus_delta * one_over_pz;
-
-    LocalParticle_add_to_x(part, px * one_over_pz * length);
-    LocalParticle_add_to_y(part, py * one_over_pz * length);
-    LocalParticle_add_to_zeta(part, dzeta * length);
-    LocalParticle_add_to_s(part, length);
-}
-
-GPUFUN
-void track_polar_drift_single_particle(
-    LocalParticle* part,  // LocalParticle to track
-    const double length,  // length of the element
-    const double h        // curvature
-) {
-
-    // Based on SUBROUTINE Sprotr in PTC and curex_drift in MAD-NG
-
-    const double rvv = LocalParticle_get_rvv(part);
-    // Particle coordinates
-    const double x = LocalParticle_get_x(part);
-    const double y = LocalParticle_get_y(part);
-    const double px = LocalParticle_get_px(part);
-    const double py = LocalParticle_get_py(part);
-    const double s = length;
-
-    const double one_plus_delta = LocalParticle_get_delta(part) + 1.0;
-    const double pz = sqrt(POW2(one_plus_delta) - POW2(px) - POW2(py));
-
-    double new_x, new_px, new_y, delta_ell;
-
-    // Polar drift
-    double const rho = 1 / h;
-    const double ca = cos(h * s);
-    const double sa = sin(h * s);
-    const double sa2 = sin(0.5 * h * s);
-    const double _pz = 1 / pz;
-    const double pxt = px * _pz;
-    const double _ptt = 1 / (ca - sa * pxt);
-    const double pst = (x + rho) * sa * _pz * _ptt;
-
-    new_x = (x + rho * (2 * sa2 * sa2 + sa * pxt)) * _ptt;
-    new_px = ca * px + sa * pz;
-    new_y = y + pst * py;
-    delta_ell = one_plus_delta * (x + rho) * sa / ca / pz / (1 - px * sa / ca / pz);
-
-    // Update Particles object
-    LocalParticle_set_x(part, new_x);
-    LocalParticle_set_px(part, new_px);
-    LocalParticle_set_y(part, new_y);
-    LocalParticle_add_to_zeta(part, length - delta_ell / rvv);
-    LocalParticle_add_to_s(part, s);
-}
-
-
-GPUFUN
-void track_expanded_combined_dipole_quad_single_particle(
-    LocalParticle* part,  // LocalParticle to track
-    const double length,  // length of the element
-    const double k0_,     // normal dipole strength
-    const double k1_,     // normal quadrupole strength
-    const double h        // curvature
-) {
-    // From madx: https://github.com/MethodicalAcceleratorDesign/MAD-X/blob/8695bd422dc403a01aa185e9fea16603bbd5b3e1/src/trrun.f90#L4320
-    // Particle coordinates
-    const double x = LocalParticle_get_x(part);
-    const double y = LocalParticle_get_y(part);
-    const double px = LocalParticle_get_px(part);
-    const double py = LocalParticle_get_py(part);
-    const double rvv = LocalParticle_get_rvv(part);
-
-    // In MAD-X (delta + 1) is computed:
-    // const double delta_plus_1 = sqrt(pt*pt + 2.0*pt*beti + 1.0);
-    const double delta_plus_1 = LocalParticle_get_delta(part) + 1;
-    const double chi = LocalParticle_get_chi(part);
-
-    const double k0 = chi * k0_ / delta_plus_1;
-    const double k1 = chi * k1_ / delta_plus_1;
-
-    const double Kx = k0 * h + k1;
-    const double Ky = -k1;
-
-    double x_, px_, y_, py_, Sx, Sy, Cx, Cy;
-
-    if (Kx > 0.0) {
-        double sqrt_Kx = sqrt(Kx);
-        Sx = sin(sqrt_Kx * length) / sqrt_Kx;
-        Cx = cos(sqrt_Kx * length);
-    }
-    else if (Kx < 0.0) {
-        double sqrt_Kx = sqrt(-Kx); // the imaginary part
-        Sx = sinh(sqrt_Kx * length) / sqrt_Kx; // sin(ix) = i sinh(x)
-        Cx = cosh(sqrt_Kx * length); // cos(ix) = cosh(x)
-    }
-    else { // Kx == 0.0
-        Sx = length;
-        Cx = 1.0;
-    }
-
-    if (Ky > 0.0) {
-        double sqrt_Ky = sqrt(Ky);
-        Sy = sin(sqrt_Ky * length) / sqrt_Ky;
-        Cy = cos(sqrt_Ky * length);
-    }
-    else if (Ky < 0.0) {
-        double sqrt_Ky = sqrt(-Ky); // the imaginary part
-        Sy = sinh(sqrt_Ky * length) / sqrt_Ky; // sin(ix) = i sinh(x)
-        Cy = cosh(sqrt_Ky * length);  // cos(ix) = cosh(x)
-    }
-    else { // Ky == 0.0
-        Sy = length;
-        Cy = 1.0;
-    }
-
-    // useful quantities
-    const double xp = px / delta_plus_1;
-    const double yp = py / delta_plus_1;
-    const double A = -Kx * x - k0 + h;
-    const double B = xp;
-    const double C = -Ky * y;
-    const double D = yp;
-
-    // transverse map
-    x_ = x * Cx + xp * Sx;
-    y_ = y * Cy + yp * Sy;
-    px_ = (A * Sx + B * Cx) * delta_plus_1;
-    py_ = (C * Sy + D * Cy) * delta_plus_1;
-
-    if (NONZERO(Kx))
-        x_ = x_ + (k0 - h) * (Cx - 1.0) / Kx;
-    else
-        x_ = x_ - (k0 - h) * 0.5 * POW2(length);
-
-    // longitudinal map
-    double length_ = length; // will be the total path length traveled by the particle
-    if (NONZERO(Kx)) {
-        length_ -= (h * ((Cx - 1.0) * xp + Sx * A + length * (k0 - h))) / Kx;
-        length_ += 0.5 * (
-            - (POW2(A) * Cx * Sx) / (2.0 * Kx) \
-            + (POW2(B) * Cx * Sx) / 2.0 \
-            + (POW2(A) * length) / (2.0 * Kx) \
-            + (POW2(B) * length) / 2.0 \
-            - (A * B * POW2(Cx)) / Kx \
-            + (A * B) / Kx
-        );
-    }
-    else {
-        length_ += h * length * (
-            3.0 * length * xp \
-            + 6.0 * x \
-            - (k0 - h) * POW2(length)
-        ) / 6.0;
-        length_ += 0.5 * (POW2(B)) * length;
-    }
-
-    if (NONZERO(Ky)) {
-        length_ += 0.5 * (
-            - (POW2(C) * Cy * Sy) / (2.0 * Ky) \
-            + (POW2(D) * Cy * Sy) / 2.0 \
-            + (POW2(C) * length) / (2.0 * Ky) \
-            + (POW2(D) * length) / 2.0 \
-            - (C * D * POW2(Cy)) / Ky \
-            + (C * D) / Ky
-        );
-    }
-    else {
-        length_ += 0.5 * POW2(D) * length;
-    }
-
-    const double dzeta = length - length_ / rvv;
-
-    LocalParticle_set_x(part, x_);
-    LocalParticle_set_px(part, px_);
-    LocalParticle_set_y(part, y_);
-    LocalParticle_set_py(part, py_);
-    LocalParticle_add_to_zeta(part, dzeta);
-    LocalParticle_add_to_s(part, length);
-
-}
-
-GPUFUN
-void track_curved_exact_bend_single_particle(
-    LocalParticle* part,  // LocalParticle to track
-    const double length,  // length of the element
-    const double k0,      // normal dipole strength
-    const double h        // curvature
-) {
-
-    // Here we assume that the caller has ensured h != 0
-
-    double const k0_chi = k0 * LocalParticle_get_chi(part);
-
-    if (fabs(k0_chi) < 1e-8) {
-        track_polar_drift_single_particle(part, length, h);
-        return;
-    }
-
-    const double rvv = LocalParticle_get_rvv(part);
-    // Particle coordinates
-    const double x = LocalParticle_get_x(part);
-    const double y = LocalParticle_get_y(part);
-    const double px = LocalParticle_get_px(part);
-    const double py = LocalParticle_get_py(part);
-    const double s = length;
-
-    const double one_plus_delta = LocalParticle_get_delta(part) + 1.0;
-    const double A = 1.0 / sqrt(POW2(one_plus_delta) - POW2(py));
-    const double pz = sqrt(POW2(one_plus_delta) - POW2(px) - POW2(py));
-
-    double new_x, new_px, new_y, delta_ell;
-
-    // The case for non-zero curvature, s is arc length
-    // Useful constants
-    const double C = pz - k0_chi * ((1 / h) + x);
-    new_px = px * cos(s * h) + C * sin(s * h);
-    double const new_pz = sqrt(POW2(one_plus_delta) - POW2(new_px) - POW2(py));
-    // double const d_new_px_ds = new_px / new_pz;
-
-    const double d_new_px_ds = C * h * cos(h * s) - h * px * sin(h * s);
-
-    // Update particle coordinates
-    new_x = (new_pz * h - d_new_px_ds - k0_chi) / (h * k0_chi);
-    const double D = asin(A * px) - asin(A * new_px);
-    new_y = y + ((py * s) / (k0_chi / h)) + (py / k0_chi) * D;
-
-    delta_ell = ((one_plus_delta * s * h) / k0_chi) + (one_plus_delta / k0_chi) * D;
-
-    // Update Particles object
-    LocalParticle_set_x(part, new_x);
-    LocalParticle_set_px(part, new_px);
-    LocalParticle_set_y(part, new_y);
-    LocalParticle_add_to_zeta(part, length - delta_ell / rvv);
-    LocalParticle_add_to_s(part, s);
-}
-
-GPUFUN
-void track_straight_exact_bend_single_particle(
-    LocalParticle* part,  // LocalParticle to track
-    const double length,  // length of the element
-    const double k0       // normal dipole strength
-) {
-
-    // Here we assume that the caller has ensured h != 0
-
-    double const k0_chi = k0 * LocalParticle_get_chi(part);
-
-    if (fabs(k0_chi) < 1e-8) {
-        track_exact_drift_single_particle(part, length);
-        return;
-    }
-
-    const double rvv = LocalParticle_get_rvv(part);
-    // Particle coordinates
-    const double x = LocalParticle_get_x(part);
-    const double y = LocalParticle_get_y(part);
-    const double px = LocalParticle_get_px(part);
-    const double py = LocalParticle_get_py(part);
-    const double s = length;
-
-    const double one_plus_delta = LocalParticle_get_delta(part) + 1.0;
-    const double A = 1.0 / sqrt(POW2(one_plus_delta) - POW2(py));
-    const double pz = sqrt(POW2(one_plus_delta) - POW2(px) - POW2(py));
-
-    double new_x, new_px, new_y, delta_ell;
-
-    // STRAIGHT EXACT BEND
-    // The case for zero curvature -- straight bend, s is Cartesian length
-    new_px = px - k0_chi * s;
-    new_x = x + (sqrt(POW2(one_plus_delta) - POW2(new_px) - POW2(py)) - pz) / k0_chi;
-
-    const double D = asin(A * px) - asin(A * new_px);
-    new_y = y + (py / k0_chi) * D;
-
-    delta_ell = (one_plus_delta / k0_chi) * D;
-
-    // Update Particles object
-    LocalParticle_set_x(part, new_x);
-    LocalParticle_set_px(part, new_px);
-    LocalParticle_set_y(part, new_y);
-    LocalParticle_add_to_zeta(part, length - delta_ell / rvv);
-    LocalParticle_add_to_s(part, s);
-}
-
-GPUFUN
-void track_solenoid_single_particle(
-    LocalParticle* part,
-    double length,
-    double ks,
-    double x0_solenoid,
-    double y0_solenoid
-) {
-    const double sk = ks / 2;
-
-    if (IS_ZERO(sk)) {
-        track_exact_drift_single_particle(part, length);
-        LocalParticle_set_ax(part, 0);
-        LocalParticle_set_ay(part, 0);
-        return;
-    }
-
-    if (IS_ZERO(length)){
-        return;
-    }
-
-    const double skl = sk * length;
-
-    // Particle coordinates
-    const double x = LocalParticle_get_x(part) - x0_solenoid;
-    const double px = LocalParticle_get_px(part);
-    const double y = LocalParticle_get_y(part) - y0_solenoid;
-    const double py = LocalParticle_get_py(part);
-    const double delta = LocalParticle_get_delta(part);
-    const double rvv = LocalParticle_get_rvv(part);
-
-    // set up constants
-    const double pk1 = px + sk * y;
-    const double pk2 = py - sk * x;
-    const double ptr2 = pk1 * pk1 + pk2 * pk2;
-    const double one_plus_delta = 1 + delta;
-    const double one_plus_delta_sq = one_plus_delta * one_plus_delta;
-    const double pz = sqrt(one_plus_delta_sq - ptr2);
-
-    // set up constants
-    const double cosTh = cos(skl / pz);
-    const double sinTh = sin(skl / pz);
-
-    const double si = sin(skl / pz) / sk;
-    const double rps[4] = {
-        cosTh * x + sinTh * y,
-        cosTh * px + sinTh * py,
-        cosTh * y - sinTh * x,
-        cosTh * py - sinTh * px
-    };
-    double const new_x = cosTh * rps[0] + si * rps[1];
-    double const new_px = cosTh * rps[1] - sk * sinTh * rps[0];
-    double const new_y = cosTh * rps[2] + si * rps[3];
-    double const new_py = cosTh * rps[3] - sk * sinTh * rps[2];
-    double const add_to_zeta = length * (1 - one_plus_delta / (pz * rvv));
-    double const new_ax = -0.5 * ks * new_y;
-    double const new_ay = 0.5 * ks * new_x;
-
-    LocalParticle_set_x(part, new_x + x0_solenoid);
-    LocalParticle_set_px(part, new_px);
-    LocalParticle_set_y(part, new_y + y0_solenoid);
-    LocalParticle_set_py(part, new_py);
-    LocalParticle_add_to_zeta(part, add_to_zeta);
-    LocalParticle_add_to_s(part, length);
-    LocalParticle_set_ax(part, new_ax);
-    LocalParticle_set_ay(part, new_ay);
-
-}
-
-
-
-GPUFUN
-void track_magnet_drift_single_particle(
-    LocalParticle* part,  // LocalParticle to track
-    const double length,  // length of the element
-    const double k0,      // normal dipole strength
-    const double k1,      // normal quadrupole strength
-    const double ks,      // solenoid strength
-    const double h,       // curvature
-    const double x0_solenoid,
-    const double y0_solenoid,
-    const int64_t drift_model      // drift model
-) {
-
-    // drift_model = 0 : drift expanded (caller has ensured k0=0, k1=0, h=0)
-    // drift_model = 1 : drift exact (caller has ensured k0=0, k1=0, h=0)
-    // drift_model = 2 : polar drift (caller has ensured k0=0, k1=0, h!=0)
-    // drift_model = 3 : k0, k1, h expanded map (this is general for all possible values)
-    // drift_model = 4 : bend with h (caller has ensured k1=0, h!=0)
-    // drift_model = 5 : bend without h (caller has ensured k1=0, h=0)
-    // drift_model = -1 : no drift
-
-    if (drift_model == -1) {
-        return;
-    }
-
-    if (length == 0.0) {
-        return;
-    }
-    switch (drift_model) {
-        case 0:
-            track_expanded_drift_single_particle(part, length);
-            break;
-        case 1:
-            track_exact_drift_single_particle(part, length);
-            break;
-        case 2:
-            track_polar_drift_single_particle(part, length, h);
-            break;
-        case 3:
-            track_expanded_combined_dipole_quad_single_particle(part, length, k0, k1, h);
-            break;
-        case 4:
-            track_curved_exact_bend_single_particle(part, length, k0, h);
-            break;
-        case 5:
-            track_straight_exact_bend_single_particle(part, length, k0);
-            break;
-        case 6:
-            track_solenoid_single_particle(part, length, ks, x0_solenoid, y0_solenoid);
-            break;
-        default:
-            break;
-    }
-
-}
-
-#undef IS_ZERO
-
-#endif
diff --git a/xtrack/beam_elements/elements_src/.__afsEE4E b/xtrack/beam_elements/elements_src/.__afsEE4E
deleted file mode 100644
index 639e04581..000000000
--- a/xtrack/beam_elements/elements_src/.__afsEE4E
+++ /dev/null
@@ -1,422 +0,0 @@
-// copyright ############################### //
-// This file is part of the Xtrack Package.  //
-// Copyright (c) CERN, 2025.                 //
-// ######################################### //
-#ifndef XTRACK_TRACK_MISALIGNMENT_H
-#define XTRACK_TRACK_MISALIGNMENT_H
-
-#include <headers/track.h>
-#include <beam_elements/elements_src/track_xrotation.h>
-#include <beam_elements/elements_src/track_yrotation.h>
-#include <beam_elements/elements_src/track_srotation.h>
-#include <beam_elements/elements_src/track_drift.h>
-#include <beam_elements/elements_src/track_xyshift.h>
-
-
-#define IF_NONZERO(VALUE, EXPRESSION) { if (VALUE != 0.0) { EXPRESSION; } }
-#define LOOP_PARTICLES(PART0, CODE) {\
-    START_PER_PARTICLE_BLOCK(PART0, part) { \
-        CODE ; \
-    } END_PER_PARTICLE_BLOCK; }
-#define Y_ROTATE(PART0, THETA) IF_NONZERO(THETA, LOOP_PARTICLES(PART0, YRotation_single_particle(part, sin(THETA), cos(THETA), tan(THETA))))
-// Flip some signs so that the input matches the MAD-X survey convention
-#define X_ROTATE(PART0, PHI) IF_NONZERO(PHI, LOOP_PARTICLES(PART0, XRotation_single_particle(part, -sin(PHI), cos(PHI), -tan(PHI))))
-#define S_ROTATE(PART0, PSI) IF_NONZERO(PSI, LOOP_PARTICLES(PART0, SRotation_single_particle(part, sin(PSI), cos(PSI))))
-#define XY_SHIFT(PART0, DX, DY) LOOP_PARTICLES(PART0, XYShift_single_particle(part, DX, DY))
-#define S_SHIFT(PART0, DS) IF_NONZERO(DS, LOOP_PARTICLES(PART0, { \
-        Drift_single_particle_exact(part, DS); \
-        LocalParticle_add_to_zeta(part, -DS); \
-        LocalParticle_add_to_s(part, -DS); \
-    }))
-
-
-GPUFUN void matrix_multiply_4x4(const double[4][4], const double[4][4], double[4][4]);
-GPUFUN void matrix_rigid_affine_inverse(const double[4][4], double[4][4]);
-
-
-GPUFUN
-void track_misalignment_entry_straight(
-    LocalParticle* part0,  // LocalParticle to track
-    double dx,  // misalignment in x
-    double dy,  // misalignment in y
-    double ds,  // misalignment in s
-    double theta, // rotation around y, yaw, positive s to x
-    double phi,  // rotation around x, pitch, positive s to y
-    double psi_no_frame,  // rotation around s, roll, positive y to x
-    double anchor, // anchor of the misalignment as offset in m from entry
-    double length,  // length of the misaligned element
-    double psi_with_frame,  // psi_with_frame of the element, positive s to x
-    int8_t backtrack
-) {
-    // Silence the warning about unused variable length
-    (void)length; // kept for API consistency with track_misalignment_exit_straight
-
-    const double mis_x = dx - anchor * cos(phi) * sin(theta);
-    const double mis_y = dy - anchor * sin(phi);
-    const double mis_s = ds - anchor * (cos(phi) * cos(theta) - 1);
-
-    // Apply transformations
-    if (!backtrack){
-        XY_SHIFT(part0, mis_x, mis_y);
-        S_SHIFT(part0, mis_s);
-        Y_ROTATE(part0, theta);
-        X_ROTATE(part0, phi);
-        S_ROTATE(part0, psi_no_frame);
-        S_ROTATE(part0, psi_with_frame);
-    } else {
-        S_ROTATE(part0, -psi_with_frame);
-        S_ROTATE(part0, -psi_no_frame);
-        X_ROTATE(part0, -phi);
-        Y_ROTATE(part0, -theta);
-        S_SHIFT(part0, -mis_s);
-        XY_SHIFT(part0, -mis_x, -mis_y);
-    }
-}
-
-
-GPUFUN
-void track_misalignment_exit_straight(
-    LocalParticle* part0,  // LocalParticle to track
-    double dx,  // misalignment in x
-    double dy,  // misalignment in y
-    double ds,  // misalignment in s
-    double theta, // rotation around y, yaw, positive s to x
-    double phi,  // rotation around x, pitch, positive s to y
-    double psi_no_frame,  // rotation around s, roll, positive y to x
-    double anchor, // anchor of the misalignment as offset in m from entry
-    double length,  // length of the misaligned element
-    double psi_with_frame,  // psi_with_frame of the element, positive s to x
-    int8_t backtrack
-) {
-    const double neg_part_length = anchor - length;
-    const double mis_x = neg_part_length * cos(phi) * sin(theta) - dx;
-    const double mis_y = neg_part_length * sin(phi) - dy;
-    const double mis_s = neg_part_length * (cos(phi) * cos(theta) - 1) - ds;
-
-    // Apply transformations
-    if (!backtrack){
-        S_ROTATE(part0, -psi_with_frame);
-        S_ROTATE(part0, -psi_no_frame);
-        X_ROTATE(part0, -phi);
-        Y_ROTATE(part0, -theta);
-        S_SHIFT(part0, mis_s);
-        XY_SHIFT(part0, mis_x, mis_y);
-    } else {
-        XY_SHIFT(part0, -mis_x, -mis_y);
-        S_SHIFT(part0, -mis_s);
-        Y_ROTATE(part0, theta);
-        X_ROTATE(part0, phi);
-        S_ROTATE(part0, psi_no_frame);
-        S_ROTATE(part0, psi_with_frame);
-    }
-}
-
-
-GPUFUN
-void track_misalignment_entry_curved(
-    LocalParticle* part0,  // LocalParticle to track
-    double dx,  // misalignment in x
-    double dy,  // misalignment in y
-    double ds,  // misalignment in s
-    double theta, // rotation around y, yaw, positive s to x
-    double phi,  // rotation around x, pitch, positive s to y
-    double psi_no_frame,  // rotation around s, roll, positive y to x
-    double anchor, // anchor of the misalignment as offset in m from entry
-    double length,  // length of the misaligned element
-    double angle,  // angle by which the element bends the reference frame
-    double h,  // curvature, only used when length == 0
-    double psi_with_frame,  // psi_with_frame of the element, positive s to x
-    int8_t backtrack
-) {
-    // Handle as a straight case if a thick element has no bending.
-    // For a thin element, we still perform the curved calculation if h is given.
-    if (angle == 0.0 && (length != 0.0 || h == 0.0)) {
-        track_misalignment_entry_straight(part0, dx, dy, ds, theta, phi,
-            psi_no_frame, anchor, length, psi_with_frame, backtrack);
-        return;
-    }
-    // Precompute trigonometric functions
-    const double s_phi = sin(phi), c_phi = cos(phi);
-    const double s_theta = sin(theta), c_theta = cos(theta);
-    const double s_psi = sin(psi_no_frame), c_psi = cos(psi_no_frame);
-
-    /* We need to compute the transformation that takes us from the aligned
-       frame to the entry of the element in the misaligned frame:
-
-       misaligned_entry = matrix_first_part @ misalignment_matrix @ inv(matrix_first_part)
-       
-       where:
-       - misalignment_matrix is the matrix that applies the misalignment
-       - matrix_first_part is the matrix that takes us from the entry of the
-         element to the anchor of the misalignment
-    */
-
-    // Misalignment matrix
-    const double misalignment_matrix[4][4] = {
-        {
-            -s_phi * s_psi * s_theta + c_psi * c_theta,
-            -c_psi * s_phi * s_theta - c_theta * s_psi,
-            c_phi * s_theta,
-            dx
-        },
-        {
-            c_phi * s_psi,
-            c_phi * c_psi,
-            s_phi, dy
-        },
-        {
-            -c_theta * s_phi * s_psi - c_psi * s_theta,
-            -c_psi * c_theta * s_phi + s_psi * s_theta,
-            c_phi * c_theta,
-            ds
-        },
-        {0, 0, 0, 1}
-    };
-
-    // Compute matrix that takes us from the reference point of the misalignment
-    // to the entry of the element
-    if (length != 0.0) h = angle / length;
-    const double part_angle = anchor * h;
-    const double delta_x_first_part = (cos(part_angle) - 1) * cos(psi_with_frame) / h;
-    const double delta_y_first_part = (cos(part_angle) - 1) * sin(psi_with_frame) / h;
-    const double delta_s_first_part = sin(part_angle) / h;
-
-    const double matrix_first_part[4][4] = {
-            {
-                (cos(part_angle) - 1) * POW2(cos(psi_with_frame)) + 1,
-                (cos(part_angle) - 1) * cos(psi_with_frame) * sin(psi_with_frame),
-                -cos(psi_with_frame) * sin(part_angle),
-                delta_x_first_part
-            },
-            {
-                (cos(part_angle) - 1) * cos(psi_with_frame) * sin(psi_with_frame),
-                (cos(part_angle) - 1) * POW2(sin(psi_with_frame)) + 1,
-                -sin(part_angle) * sin(psi_with_frame),
-                delta_y_first_part
-            },
-            {
-                cos(psi_with_frame) * sin(part_angle),
-                sin(part_angle) * sin(psi_with_frame),
-                cos(part_angle),
-                delta_s_first_part
-            },
-            {0, 0, 0, 1}
-        };
-
-    double inv_matrix_first_part[4][4];
-    matrix_rigid_affine_inverse(matrix_first_part, inv_matrix_first_part);
-
-    // Compute the transformation that takes us from the aligned frame to the
-    // entry of the element in the misaligned frame
-    double misaligned_entry[4][4], temp[4][4];
-    matrix_multiply_4x4(matrix_first_part, misalignment_matrix, temp);
-    matrix_multiply_4x4(temp, inv_matrix_first_part, misaligned_entry);
-
-    // Extract the basic transformations from the misalignment matrix
-    const double mis_x = misaligned_entry[0][3];
-    const double mis_y = misaligned_entry[1][3];
-    const double mis_s = misaligned_entry[2][3];
-
-    const double rot_theta = atan2(misaligned_entry[0][2], misaligned_entry[2][2]);
-    const double rot_phi = atan2(misaligned_entry[1][2], sqrt(misaligned_entry[1][0] * misaligned_entry[1][0] + misaligned_entry[1][1] * misaligned_entry[1][1]));
-    const double rot_psi = atan2(misaligned_entry[1][0], misaligned_entry[1][1]);
-
-    // Apply transformations
-    if (!backtrack) {
-        XY_SHIFT(part0, mis_x, mis_y);
-        S_SHIFT(part0, mis_s);
-        Y_ROTATE(part0, rot_theta);
-        X_ROTATE(part0, rot_phi);
-        S_ROTATE(part0, rot_psi);
-        S_ROTATE(part0, psi_with_frame);
-    } else {
-        S_ROTATE(part0, -psi_with_frame);
-        S_ROTATE(part0, -rot_psi);
-        X_ROTATE(part0, -rot_phi);
-        Y_ROTATE(part0, -rot_theta);
-        S_SHIFT(part0, -mis_s);
-        XY_SHIFT(part0, -mis_x, -mis_y);
-    }
-}
-
-
-GPUFUN
-void track_misalignment_exit_curved(
-    LocalParticle* part0,  // LocalParticle to track
-    double dx,  // misalignment in x
-    double dy,  // misalignment in y
-    double ds,  // misalignment in s
-    double theta, // rotation around y, yaw, positive s to x
-    double phi,  // rotation around x, pitch, positive s to y
-    double psi_no_frame,  // rotation around s, roll, positive y to x
-    double anchor, // anchor of the misalignment as a fraction of the length
-    double length,  // length of the misaligned element
-    double angle,  // angle by which the element bends the reference frame
-    double h,  // curvature, only used when length == 0
-    double psi_with_frame,  // psi_with_frame of the element, positive s to x
-    int8_t backtrack  // whether to backtrack the particle
-) {
-    // Handle as a straight case if a thick element has no bending.
-    // For a thin element, we still perform the curved calculation if h is given.
-    if (angle == 0.0 && (length != 0.0 || h == 0.0)) {
-        track_misalignment_exit_straight(
-            part0, dx, dy, ds, theta, phi, psi_no_frame, anchor, length,
-            psi_with_frame, backtrack);
-        return;
-    }
-    // Precompute trigonometric functions
-    double s_phi = sin(phi), c_phi = cos(phi);
-    double s_theta = sin(theta), c_theta = cos(theta);
-    double s_psi = sin(psi_no_frame), c_psi = cos(psi_no_frame);
-
-    /* We need to compute the transformation that takes us from the misaligned
-       frame to the aligned frame:
-
-       realign = inv(matrix_second_part) * inv(misalignment_matrix) * matrix_second_part
-
-       where:
-       - misalignment_matrix is the matrix that applies the misalignment
-       - matrix_second_part is the matrix that takes us from the frame in the
-         middle of the element (anchor) to the end of the element
-    */
-
-    // Misalignment matrix
-    const double misalignment_matrix[4][4] = {
-        {
-            -s_phi * s_psi * s_theta + c_psi * c_theta,
-            -c_psi * s_phi * s_theta - c_theta * s_psi,
-            c_phi * s_theta,
-            dx
-        },
-        {
-            c_phi * s_psi,
-            c_phi * c_psi,
-            s_phi,
-            dy
-        },
-        {
-            -c_theta * s_phi * s_psi - c_psi * s_theta,
-            -c_psi * c_theta * s_phi + s_psi * s_theta,
-            c_phi * c_theta,
-            ds
-        },
-        {0, 0, 0, 1}
-    };
-
-    // Compute the inverse of the misalignment matrix
-    double inv_misalignment_matrix[4][4];
-    matrix_rigid_affine_inverse(misalignment_matrix, inv_misalignment_matrix);
-
-    // Compute the inverse of the matrix that takes us from the point of the
-    // misalignment to the exit of the element.
-    if (length != 0.0) h = angle / length;
-    const double part_angle = angle - h * anchor;
-
-    const double s_part_angle = sin(part_angle), c_part_angle = cos(part_angle);
-    const double s_tilt = sin(psi_with_frame), c_tilt = cos(psi_with_frame);
-
-    const double delta_x_second_part = (c_part_angle - 1) * c_tilt / h;
-    const double delta_y_second_part = (c_part_angle - 1) * s_tilt / h;
-    const double delta_s_second_part = s_part_angle / h;
-
-    const double matrix_second_part[4][4] = {
-        {
-            (c_part_angle - 1) * POW2(c_tilt) + 1,
-            (c_part_angle - 1) * c_tilt * s_tilt,
-            c_tilt * -s_part_angle,
-            delta_x_second_part
-        },
-        {
-            (c_part_angle - 1) * c_tilt * s_tilt,
-            (c_part_angle - 1) * POW2(s_tilt) + 1,
-            -s_part_angle * s_tilt,
-            delta_y_second_part
-        },
-        {
-            c_tilt * s_part_angle,
-            s_part_angle * s_tilt,
-            c_part_angle,
-            delta_s_second_part
-        },
-        {0, 0, 0, 1}
-    };
-
-    double inv_matrix_second_part[4][4];
-    matrix_rigid_affine_inverse(matrix_second_part, inv_matrix_second_part);
-
-    // Compute the realignment matrix
-    double realign[4][4], temp[4][4];
-    matrix_multiply_4x4(inv_matrix_second_part, inv_misalignment_matrix, temp);
-    matrix_multiply_4x4(temp, matrix_second_part, realign);
-
-    // Extract the basic transformations from the realignment matrix
-    double mis_x = realign[0][3];
-    double mis_y = realign[1][3];
-    double mis_s = realign[2][3];
-
-    double rot_theta = atan2(realign[0][2], realign[2][2]);
-    double rot_phi = atan2(realign[1][2], sqrt(realign[1][0] * realign[1][0] + realign[1][1] * realign[1][1]));
-    double rot_psi = atan2(realign[1][0], realign[1][1]);
-
-    // Apply transformations
-    if (!backtrack){
-        S_ROTATE(part0, -psi_with_frame);
-        XY_SHIFT(part0, mis_x, mis_y);
-        S_SHIFT(part0, mis_s);
-        Y_ROTATE(part0, rot_theta);
-        X_ROTATE(part0, rot_phi);
-        S_ROTATE(part0, rot_psi);
-    } else {
-        S_ROTATE(part0, -rot_psi);
-        X_ROTATE(part0, -rot_phi);
-        Y_ROTATE(part0, -rot_theta);
-        S_SHIFT(part0, -mis_s);
-        XY_SHIFT(part0, -mis_x, -mis_y);
-        S_ROTATE(part0, psi_with_frame);
-    }
-}
-
-
-GPUFUN
-void matrix_multiply_4x4(const double a[4][4], const double b[4][4], double result[4][4]) {
-    // Multiply two 4x4 matrices `a` and `b`, and store the result in `result`.
-    for (int i = 0; i < 4; i++) {
-        result[i][0] = a[i][0] * b[0][0] + a[i][1] * b[1][0] + a[i][2] * b[2][0] + a[i][3] * b[3][0];
-        result[i][1] = a[i][0] * b[0][1] + a[i][1] * b[1][1] + a[i][2] * b[2][1] + a[i][3] * b[3][1];
-        result[i][2] = a[i][0] * b[0][2] + a[i][1] * b[1][2] + a[i][2] * b[2][2] + a[i][3] * b[3][2];
-        result[i][3] = a[i][0] * b[0][3] + a[i][1] * b[1][3] + a[i][2] * b[2][3] + a[i][3] * b[3][3];
-    }
-}
-
-
-GPUFUN
-void matrix_rigid_affine_inverse(const double m[4][4], double inv_m[4][4]) {
-    // Compute the inverse `inv_m` of a rigid affine transformation matrix `m`.
-    const double m00 = m[0][0], m01 = m[0][1], m02 = m[0][2], m03 = m[0][3];
-    const double m10 = m[1][0], m11 = m[1][1], m12 = m[1][2], m13 = m[1][3];
-    const double m20 = m[2][0], m21 = m[2][1], m22 = m[2][2], m23 = m[2][3];
-
-    // Invert the rotation part of `m`: as it's rigid, it's simply a transpose
-    inv_m[0][0] = m00;
-    inv_m[0][1] = m10;
-    inv_m[0][2] = m20;
-    inv_m[1][0] = m01;
-    inv_m[1][1] = m11;
-    inv_m[1][2] = m21;
-    inv_m[2][0] = m02;
-    inv_m[2][1] = m12;
-    inv_m[2][2] = m22;
-
-    // Compute the translation part, -m[:3, :3]^{-1} @ m[:3, 3]
-    inv_m[0][3] = -m00 * m03 - m10 * m13 - m20 * m23;
-    inv_m[1][3] = -m01 * m03 - m11 * m13 - m21 * m23;
-    inv_m[2][3] = -m02 * m03 - m12 * m13 - m22 * m23;
-
-    // Fill the last row
-    inv_m[3][0] = 0;
-    inv_m[3][1] = 0;
-    inv_m[3][2] = 0;
-    inv_m[3][3] = 1;
-}
-
-#endif  // XTRACK_TRACK_MISALIGNMENT_H
diff --git a/xtrack/beam_elements/elements_src/.__afsF1AA b/xtrack/beam_elements/elements_src/.__afsF1AA
deleted file mode 100644
index 12e929645..000000000
--- a/xtrack/beam_elements/elements_src/.__afsF1AA
+++ /dev/null
@@ -1,635 +0,0 @@
-// copyright ############################### //
-// This file is part of the Xtrack Package.  //
-// Copyright (c) CERN, 2021.                 //
-// ######################################### //
-
-#ifndef XTRACK_LINESEGMENTMAP_H
-#define XTRACK_LINESEGMENTMAP_H
-
-#include <headers/track.h>
-
-
-GPUFUN
-void remove_closed_orbit(
-    LocalParticle* part0, double const x_ref, double const px_ref,
-    double const y_ref, double const py_ref){
-
-    START_PER_PARTICLE_BLOCK(part0, part);
-        // Remove closed orbit
-        LocalParticle_add_to_x(part, -x_ref);
-        LocalParticle_add_to_px(part, -px_ref);
-        LocalParticle_add_to_y(part, -y_ref);
-        LocalParticle_add_to_py(part, -py_ref);
-    END_PER_PARTICLE_BLOCK;
-}
-
-GPUFUN
-void add_closed_orbit(
-    LocalParticle* part0, double const x_ref, double const px_ref,
-    double const y_ref, double const py_ref){
-
-    START_PER_PARTICLE_BLOCK(part0, part);
-        // Add closed orbit
-        LocalParticle_add_to_x(part, x_ref);
-        LocalParticle_add_to_px(part, px_ref);
-        LocalParticle_add_to_y(part, y_ref);
-        LocalParticle_add_to_py(part, py_ref);
-    END_PER_PARTICLE_BLOCK;
-}
-
-GPUFUN
-void remove_dispersion(
-    LocalParticle* part0, double const dx_0, double const dpx_0,
-    double const dy_0, double const dpy_0){
-
-    START_PER_PARTICLE_BLOCK(part0, part);
-        // Remove dispersion
-        // Symplecticity correction (not working, to be investigated)
-        // LocalParticle_add_to_zeta(part, (
-        //     dpx_0 * LocalParticle_get_x(part)
-        //     - dx_0 * LocalParticle_get_px(part)
-        //     + dpy_0 * LocalParticle_get_y(part)
-        //     - dy_0 * LocalParticle_get_py(part)
-        //     )/LocalParticle_get_rvv(part));
-        double const delta = LocalParticle_get_delta(part);
-        LocalParticle_add_to_x(part, -dx_0 * delta);
-        LocalParticle_add_to_px(part, -dpx_0 * delta);
-        LocalParticle_add_to_y(part, -dy_0 * delta);
-        LocalParticle_add_to_py(part, -dpy_0 * delta);
-    END_PER_PARTICLE_BLOCK;
-}
-
-GPUFUN
-void add_dispersion(
-    LocalParticle* part0, double const dx_1, double const dpx_1,
-    double const dy_1, double const dpy_1){
-
-    START_PER_PARTICLE_BLOCK(part0, part);
-        // Add dispersion
-        // Symplecticity correction (not working, to be investigated)
-        // LocalParticle_add_to_zeta(part, (
-        //     dpx_1 * LocalParticle_get_x(part)
-        //     - dx_1 * LocalParticle_get_px(part)
-        //     + dpy_1 * LocalParticle_get_y(part)
-        //     - dy_1 * LocalParticle_get_py(part)
-        //     )/LocalParticle_get_rvv(part));
-        double const delta = LocalParticle_get_delta(part);
-        LocalParticle_add_to_x(part, dx_1 * delta);
-        LocalParticle_add_to_px(part, dpx_1 * delta);
-        LocalParticle_add_to_y(part, dy_1 * delta);
-        LocalParticle_add_to_py(part, dpy_1 * delta);
-    END_PER_PARTICLE_BLOCK;
-}
-
-GPUFUN
-void transverse_motion(LocalParticle *part0, LineSegmentMapData el){
-
-    double const qx = LineSegmentMapData_get_qx(el);
-    double const qy = LineSegmentMapData_get_qy(el);
-    double const det_xx = LineSegmentMapData_get_det_xx(el);
-    double const det_xy = LineSegmentMapData_get_det_xy(el);
-    double const det_yx = LineSegmentMapData_get_det_yx(el);
-    double const det_yy = LineSegmentMapData_get_det_yy(el);
-    double const alfx_0 = LineSegmentMapData_get_alfx(el, 0);
-    double const betx_0 = LineSegmentMapData_get_betx(el, 0);
-    double const alfy_0 = LineSegmentMapData_get_alfy(el, 0);
-    double const bety_0 = LineSegmentMapData_get_bety(el, 0);
-    double const alfx_1 = LineSegmentMapData_get_alfx(el, 1);
-    double const betx_1 = LineSegmentMapData_get_betx(el, 1);
-    double const alfy_1 = LineSegmentMapData_get_alfy(el, 1);
-    double const bety_1 = LineSegmentMapData_get_bety(el, 1);
-    int64_t const ndqx = LineSegmentMapData_len_coeffs_dqx(el);
-    int64_t const ndqy = LineSegmentMapData_len_coeffs_dqy(el);
-
-    double const c11_0= LineSegmentMapData_get_c11(el,0);
-    double const c12_0= LineSegmentMapData_get_c12(el,0);
-    double const c21_0= LineSegmentMapData_get_c21(el,0);
-    double const c22_0= LineSegmentMapData_get_c22(el,0);
-
-    double const c11_1= LineSegmentMapData_get_c11(el,1);;
-    double const c12_1= LineSegmentMapData_get_c12(el,1);;
-    double const c21_1= LineSegmentMapData_get_c21(el,1);;
-    double const c22_1= LineSegmentMapData_get_c22(el,1);
-
-    double const coupling_0 = c11_0 + c22_0 + c12_0 + c21_0;
-    double const coupling_1 = c11_1 + c22_1 + c12_1 + c21_1;
-
-    int64_t detuning;
-    double sin_x = 0;
-    double cos_x = 0;
-    double sin_y = 0;
-    double cos_y = 0;
-
-    int64_t any_chroma = 0;
-
-    for (int i_dqx=0; i_dqx<ndqx; i_dqx++){
-        any_chroma += LineSegmentMapData_get_coeffs_dqx(el, 0) != 0;
-    }
-
-    for (int i_dqy=0; i_dqy<ndqy; i_dqy++){
-        any_chroma += LineSegmentMapData_get_coeffs_dqy(el, 0) != 0;
-    }
-
-    if (any_chroma ||
-        det_xx != 0.0 || det_xy != 0.0 || det_yx != 0.0 || det_yy != 0.0){
-        detuning = 1;
-    }
-    else{
-        detuning = 0;
-        sin_x = sin(2 * PI * qx);
-        cos_x = cos(2 * PI * qx);
-        sin_y = sin(2 * PI * qy);
-        cos_y = cos(2 * PI * qy);
-    }
-
-    double const sqrt_betprod_x = sqrt(betx_1 * betx_0);
-    double const sqrt_betprod_y = sqrt(bety_1 * bety_0);
-
-    double const sqrt_betratio_x = sqrt(betx_1 / betx_0);
-    double const sqrt_betratio_y = sqrt(bety_1 / bety_0);
-
-    START_PER_PARTICLE_BLOCK(part0, part);
-        if (detuning){
-            double const J_x = 0.5 * (
-                (1.0 + alfx_0 * alfx_0) / betx_0
-                    * LocalParticle_get_x(part)*LocalParticle_get_x(part)
-                + 2 * alfx_0
-                    * LocalParticle_get_x(part)*LocalParticle_get_px(part)
-                + betx_0
-                    * LocalParticle_get_px(part)*LocalParticle_get_px(part));
-            double const J_y = 0.5 * (
-                (1.0 + alfy_0 * alfy_0) /bety_0
-                    * LocalParticle_get_y(part)*LocalParticle_get_y(part)
-                + 2*alfy_0
-                    * LocalParticle_get_y(part)*LocalParticle_get_py(part)
-                + bety_0
-                    * LocalParticle_get_py(part)*LocalParticle_get_py(part));
-            double phase = 2*PI*(qx + det_xx * J_x + det_xy * J_y);
-            for (int i_dqx=1; i_dqx<ndqx; i_dqx++){
-                phase += 2*PI*(LineSegmentMapData_get_coeffs_dqx(el, i_dqx) *
-                               pow(LocalParticle_get_delta(part), (double)i_dqx));
-            }
-            cos_x = cos(phase);
-            sin_x = sin(phase);
-            phase = 2*PI*(qy + det_yx * J_x + det_yy * J_y);
-            for (int i_dqy=1; i_dqy<ndqy; i_dqy++){
-                phase += 2*PI*(LineSegmentMapData_get_coeffs_dqy(el, i_dqy) *
-                               pow(LocalParticle_get_delta(part), (double)i_dqy));
-            }
-            cos_y = cos(phase);
-            sin_y = sin(phase);
-        }
-
-        double const M00_x = sqrt_betratio_x*(cos_x+alfx_0*sin_x);
-        double const M01_x = sqrt_betprod_x*sin_x;
-        double const M10_x = ((alfx_0-alfx_1)*cos_x
-                    -(1+alfx_0*alfx_1)*sin_x
-                    )/sqrt_betprod_x;
-        double const M11_x = (cos_x-alfx_1*sin_x)/sqrt_betratio_x;
-        double const M00_y = sqrt_betratio_y*(cos_y+alfy_0*sin_y);
-        double const M01_y = sqrt_betprod_y*sin_y;
-        double const M10_y = ((alfy_0-alfy_1)*cos_y
-                    -(1+alfy_0*alfy_1)*sin_y
-                    )/sqrt_betprod_y;
-        double const M11_y = (cos_y-alfy_1*sin_y)/sqrt_betratio_y;
-
-        if (coupling_0 == 0.0 && coupling_1 == 0.0) {
-            // No coupling matrices, use simple transfer matrix
-            double const x_out = M00_x*LocalParticle_get_x(part) + M01_x * LocalParticle_get_px(part);
-            double const px_out = M10_x*LocalParticle_get_x(part) + M11_x * LocalParticle_get_px(part);
-            double const y_out = M00_y*LocalParticle_get_y(part) + M01_y * LocalParticle_get_py(part);
-            double const py_out = M10_y*LocalParticle_get_y(part) + M11_y * LocalParticle_get_py(part);
-
-            LocalParticle_set_x(part, x_out);
-            LocalParticle_set_px(part, px_out);
-            LocalParticle_set_y(part, y_out);
-            LocalParticle_set_py(part, py_out);
-        }
-
-        /* Construct 4x4 symplectic matrix R and its inverse from a 2x2 matrix C */
-        void compute_matrices_V(double const c11, double const c12,
-                          double const c21, double const c22,
-                          double const gamma,
-                          double V[16]) 
-        {
-        /* Initialize R and Rinv to zero */
-            for(int i=0;i<16;i++) {
-            V[i] = 0.0;
-            }
-            /* Identity 2x2 blocks (top-left and bottom-right) */
-            V[0] = gamma;  V[5] = gamma;
-            V[10] = gamma; V[15] = gamma;
-            /* C block (top-right) */
-            V[2] = c11;  V[3] = -c12;
-            V[6] = -c21;  V[7] = c22;
-            /* -C^inv block (bottom-left) */
-            V[8]  = -c22; V[9]  =  -c12;
-            V[12] =  -c21; V[13] = -c11;
-        }
-
-        /* Construct 4x4 symplectic matrix R and its inverse from a 2x2 matrix C */
-        void compute_matrices_Vinv(double const c11, double const c12,
-                          double const c21, double const c22,
-                          double const gamma,
-                          double V_inv[16]) 
-        {
-            /* Initialize R and Rinv to zero */
-            for(int i=0;i<16;i++) {
-                V_inv[i] = 0.0;
-            }
-            /* Identity 2x2 blocks (top-left and bottom-right) */
-            V_inv[0] = gamma;  V_inv[5] = gamma;
-            V_inv[10] = gamma; V_inv[15] = gamma;
-            /* -C block (top-right) */
-            V_inv[2] = -c11;  V_inv[3] =  c12;
-            V_inv[6] =  c21;  V_inv[7] = -c22;
-            /*C^inv block (bottom-left) */
-            V_inv[8]  =  c22; V_inv[9]  =  c12;
-            V_inv[12] =  c21; V_inv[13] =  c11;
-        }
-
-
-        if (coupling_0 != 0.0 || coupling_1 != 0.0) {
-        
-            double const det_0 = c11_0 * c22_0 - c12_0 * c21_0;
-            double const det_1 = c11_1 * c22_1 - c12_1 * c21_1;
-
-            double const gamma_0 = sqrt(1 - det_0);
-            double const gamma_1 = sqrt(1 - det_1);
-
-            double V_1[16], V_0inv[16];
-            double M[16];
-            double T[16];
-            double temp[16];
-
-            compute_matrices_V(c11_1, c12_1, c21_1, c22_1, gamma_1, V_1);  
-            compute_matrices_Vinv(c11_0, c12_0, c21_0, c22_0, gamma_0, V_0inv);
-            
-            //the final transformation is V_1 * M * V_0inv
-
-            for(int i=0;i<16;i++) {
-                M[i] = 0.0;
-            }
-            M[0] = M00_x; M[1] = M01_x; M[2]=0.0; M[3]=0.0;
-            M[4] = M10_x; M[5] = M11_x; M[6]=0.0; M[7]=0.0;
-            M[8] = 0.0; M[9] = 0.0; M[10]=M00_y; M[11]=M01_y;
-            M[12]= 0.0; M[13]= 0.0; M[14]=M10_y; M[15]=M11_y;       
-
-            for(int i=0;i<16;i++) {
-                T[i] = 0.0;
-            }
-
-            for(int i=0;i<16;i++) {
-                temp[i] = 0.0;
-            }
-
-            // temp = M * V_0inv
-            for (int i=0;i<4;i++) {
-                for (int j=0;j<4;j++) {
-                    temp[i*4+j] = 0.0;
-                    for (int k=0;k<4;k++)
-                        temp[i*4+j] += M[i*4+k]*V_0inv[k*4+j];
-                }
-            }
-            // T = V_1 * temp
-            for (int i=0;i<4;i++) {
-                for (int j=0;j<4;j++) {
-                    T[i*4+j] = 0.0;
-                    for (int k=0;k<4;k++)
-                        T[i*4+j] += V_1[i*4+k]*temp[k*4+j];
-                }
-            }
-
-                double const old_x = LocalParticle_get_x(part);
-                double const old_px = LocalParticle_get_px(part);
-                double const old_y = LocalParticle_get_y(part);
-                double const old_py = LocalParticle_get_py(part);
-
-                double const x_out = (
-                    T[0]*old_x + T[1]*old_px + T[2]*old_y + T[3]*old_py);
-                double const px_out = (
-                    T[4]*old_x + T[5]*old_px + T[6]*old_y + T[7]*old_py);
-                double const y_out = (
-                    T[8]*old_x + T[9]*old_px + T[10]*old_y + T[11]*old_py);
-                double const py_out = (
-                    T[12]*old_x + T[13]*old_px + T[14]*old_y + T[15]*old_py);
-
-                LocalParticle_set_x(part, x_out);
-                LocalParticle_set_px(part, px_out);
-                LocalParticle_set_y(part, y_out);
-                LocalParticle_set_py(part, py_out);
-            END_PER_PARTICLE_BLOCK;
-    }
-}
-
-
-GPUFUN
-void longitudinal_motion(LocalParticle *part0,
-                         LineSegmentMapData el){
-
-    int64_t const mode_flag = LineSegmentMapData_get_longitudinal_mode_flag(el);
-
-    if (mode_flag==1){ // linear motion fixed qs
-        double const qs = LineSegmentMapData_get_qs(el);
-        double const bets = LineSegmentMapData_get_bets(el);
-        double const bucket_length = LineSegmentMapData_get_bucket_length(el)*LocalParticle_get_beta0(part0)*C_LIGHT;
-        double const sin_s = sin(2 * PI * qs);
-        double const cos_s = cos(2 * PI * qs);
-        uint8_t kill_energy_kick = LocalParticle_check_track_flag(
-                    part0, XS_FLAG_KILL_CAVITY_KICK);
-        START_PER_PARTICLE_BLOCK(part0, part);
-            // We set cos_s = 999 if long map is to be skipped
-            double shift = 0.0;
-            if (bucket_length > 0.0) {
-                shift = bucket_length*floor(LocalParticle_get_zeta(part)/bucket_length+0.5);
-                LocalParticle_add_to_zeta(part,-shift);
-            }
-            double const new_zeta = cos_s * LocalParticle_get_zeta(part) - bets * sin_s * LocalParticle_get_pzeta(part);
-            double const new_pzeta = sin_s * LocalParticle_get_zeta(part) / bets + cos_s * LocalParticle_get_pzeta(part);
-
-            LocalParticle_set_zeta(part, new_zeta+shift);
-            if (!kill_energy_kick){
-                LocalParticle_update_pzeta(part, new_pzeta);
-            }
-        END_PER_PARTICLE_BLOCK;
-    }
-    else if (mode_flag==2){ // non-linear motion
-
-        double const alfp =
-            LineSegmentMapData_get_momentum_compaction_factor(el);
-        double const slippage_length =
-            LineSegmentMapData_get_slippage_length(el);
-
-        START_PER_PARTICLE_BLOCK(part0, part);
-            double const gamma0 = LocalParticle_get_gamma0(part);
-            double const eta = alfp - 1.0 / (gamma0 * gamma0);
-            LocalParticle_add_to_zeta(part,
-                -0.5 * eta * slippage_length * LocalParticle_get_delta(part));
-        END_PER_PARTICLE_BLOCK;
-
-        int64_t const n_rf = LineSegmentMapData_len_voltage_rf(el);
-        for (int i_rf=0; i_rf<n_rf; i_rf++){
-
-            double const v_rf = LineSegmentMapData_get_voltage_rf(el,i_rf);
-            double const f_rf = LineSegmentMapData_get_frequency_rf(el,i_rf);
-            double const lag_rf = LineSegmentMapData_get_lag_rf(el,i_rf);
-
-            if (f_rf == 0) continue;
-
-            START_PER_PARTICLE_BLOCK(part0, part);
-                double const K_FACTOR = ( ( double )2.0 *PI ) / C_LIGHT;
-                double const   beta0  = LocalParticle_get_beta0(part);
-                double const   zeta   = LocalParticle_get_zeta(part);
-                double const   q      = fabs(LocalParticle_get_q0(part))
-                                        * LocalParticle_get_charge_ratio(part);
-                double const   tau    = zeta / beta0;
-                double const   phase  = DEG2RAD  * lag_rf - K_FACTOR * f_rf * tau;
-                double const energy   = q * v_rf * sin(phase);
-                LocalParticle_add_to_energy(part, energy, 1);
-            END_PER_PARTICLE_BLOCK;
-        }
-
-        START_PER_PARTICLE_BLOCK(part0, part);
-            double const gamma0 = LocalParticle_get_gamma0(part);
-            double const eta = alfp - 1.0 / (gamma0 * gamma0);
-            LocalParticle_add_to_zeta(part,
-                -0.5 * eta * slippage_length * LocalParticle_get_delta(part));
-        END_PER_PARTICLE_BLOCK;
-    }
-    else if (mode_flag == 3){ // linear motion fixed RF
-
-        double const alfp =
-            LineSegmentMapData_get_momentum_compaction_factor(el);
-        double const slippage_length =
-            LineSegmentMapData_get_slippage_length(el);
-
-        // Assume there is only one RF term (checked in the Python code)
-        double const v_rf = LineSegmentMapData_get_voltage_rf(el,0);
-        double const f_rf = LineSegmentMapData_get_frequency_rf(el,0);
-
-        double const bucket_length = LocalParticle_get_beta0(part0)*C_LIGHT/f_rf;
-
-        START_PER_PARTICLE_BLOCK(part0, part);
-            double const gamma0 = LocalParticle_get_gamma0(part);
-            double const beta0 = LocalParticle_get_beta0(part);
-            double const q0 = LocalParticle_get_q0(part);
-            double const mass0 = LocalParticle_get_mass0(part);
-            double const eta = alfp - 1.0 / (gamma0 * gamma0);
-            double const E0 = mass0 * gamma0;
-
-            double const qs = sqrt(q0 * fabs(eta) * slippage_length * f_rf * v_rf
-                        / (2 * PI * beta0 * beta0 * beta0 * E0 * C_LIGHT));
-            double const bets = eta * slippage_length / (2 * PI * qs);
-
-            double const sin_s = sin(2 * PI * qs);
-            double const cos_s = cos(2 * PI * qs);
-
-            double shift = bucket_length*floor(LocalParticle_get_zeta(part)/bucket_length+0.5);
-            LocalParticle_add_to_zeta(part,-shift);
-
-            double const new_zeta = cos_s * LocalParticle_get_zeta(part) - bets * sin_s * LocalParticle_get_pzeta(part);
-            double const new_pzeta = sin_s * LocalParticle_get_zeta(part) / bets + cos_s * LocalParticle_get_pzeta(part);
-
-            LocalParticle_set_zeta(part, new_zeta+shift);
-            LocalParticle_update_pzeta(part, new_pzeta);
-        END_PER_PARTICLE_BLOCK;
-    }
-}
-
-GPUFUN
-void uncorrelated_radiation_damping(LocalParticle *part0,
-            LineSegmentMapData el){
-
-    START_PER_PARTICLE_BLOCK(part0, part);
-        LocalParticle_set_x(part,
-            LineSegmentMapData_get_damping_factors(el,0,0)*LocalParticle_get_x(part));
-        LocalParticle_set_px(part,
-            LineSegmentMapData_get_damping_factors(el,1,1)*LocalParticle_get_px(part));
-        LocalParticle_set_y(part,
-            LineSegmentMapData_get_damping_factors(el,2,2)*LocalParticle_get_y(part));
-        LocalParticle_set_py(part,
-            LineSegmentMapData_get_damping_factors(el,3,3)*LocalParticle_get_py(part));
-        LocalParticle_set_zeta(part,
-            LineSegmentMapData_get_damping_factors(el,4,4)*LocalParticle_get_zeta(part));
-        LocalParticle_update_pzeta(part,
-            LineSegmentMapData_get_damping_factors(el,5,5)*LocalParticle_get_pzeta(part));
-    END_PER_PARTICLE_BLOCK;
-}
-
-GPUFUN
-void correlated_radiation_damping(LocalParticle *part0,
-            LineSegmentMapData el){
-
-    uint8_t kill_energy_kick = LocalParticle_check_track_flag(
-                    part0, XS_FLAG_KILL_CAVITY_KICK);
-    START_PER_PARTICLE_BLOCK(part0, part);
-        double in[6];
-        in[0] = LocalParticle_get_x(part);
-        in[1] = LocalParticle_get_px(part);
-        in[2] = LocalParticle_get_y(part);
-        in[3] = LocalParticle_get_py(part);
-        in[4] = LocalParticle_get_zeta(part);
-        in[5] = LocalParticle_get_pzeta(part);
-        double out[6];
-        for(unsigned int i=0;i<6;++i){
-            out[i] = 0;
-            for(unsigned int j=0;j<6;++j){
-                out[i] += LineSegmentMapData_get_damping_factors(el,i,j)*in[j];
-            }
-        }
-        LocalParticle_set_x(part, out[0]);
-        LocalParticle_set_px(part, out[1]);
-        LocalParticle_set_y(part, out[2]);
-        LocalParticle_set_py(part, out[3]);
-        LocalParticle_set_zeta(part, out[4]);
-        if (!kill_energy_kick){
-            LocalParticle_update_pzeta(part,out[5]);}
-    END_PER_PARTICLE_BLOCK;
-}
-
-GPUFUN
-void energy_and_reference_increments(LocalParticle *part0,
-    double const energy_increment, double const energy_ref_increment){
-
-    START_PER_PARTICLE_BLOCK(part0, part);
-        // Change energy without change of reference momentum
-        if (energy_increment !=0){
-        LocalParticle_add_to_energy(part, energy_increment, 1);
-        }
-
-        // Change energy reference
-        // In the transverse plane de change is smoothed, i.e.
-        // both the position and the momentum are scaled,
-        // rather than only the momentum.
-        if (energy_ref_increment != 0){
-            double const old_px = LocalParticle_get_px(part);
-            double const old_py = LocalParticle_get_py(part);
-            double const new_energy0 = LocalParticle_get_mass0(part)
-                *LocalParticle_get_gamma0(part) + energy_ref_increment;
-            double const new_p0c = sqrt(new_energy0*new_energy0
-            -LocalParticle_get_mass0(part)*LocalParticle_get_mass0(part));
-            double const new_beta0 = new_p0c / new_energy0;
-            double const new_gamma0 = new_energy0 / LocalParticle_get_mass0(part);
-            double const geo_emit_factor = sqrt(LocalParticle_get_beta0(part)
-                    *LocalParticle_get_gamma0(part)/new_beta0/new_gamma0);
-            LocalParticle_update_p0c(part, new_p0c); // updates px, py but not in the smoothed way
-            LocalParticle_set_px(part, old_px * geo_emit_factor);
-            LocalParticle_set_py(part, old_py * geo_emit_factor);
-            LocalParticle_scale_x(part, geo_emit_factor);
-            LocalParticle_scale_y(part, geo_emit_factor);
-        }
-    END_PER_PARTICLE_BLOCK;
-}
-
-GPUFUN
-void uncorrelated_gaussian_noise(LocalParticle *part0,
-                    LineSegmentMapData el){
-
-        START_PER_PARTICLE_BLOCK(part0, part);
-            double r = RandomNormal_generate(part);
-            LocalParticle_add_to_x(part,
-                r*LineSegmentMapData_get_gauss_noise_matrix(el,0,0));
-            r = RandomNormal_generate(part);
-            LocalParticle_add_to_px(part,
-                r*LineSegmentMapData_get_gauss_noise_matrix(el,1,1));
-            r = RandomNormal_generate(part);
-            LocalParticle_add_to_y(part,
-                r*LineSegmentMapData_get_gauss_noise_matrix(el,2,2));
-            r = RandomNormal_generate(part);
-            LocalParticle_add_to_py(part,
-                r*LineSegmentMapData_get_gauss_noise_matrix(el,3,3));
-            r = RandomNormal_generate(part);
-            LocalParticle_add_to_zeta(part,
-                r*LineSegmentMapData_get_gauss_noise_matrix(el,4,4));
-            r = RandomNormal_generate(part);
-            double pzeta = LocalParticle_get_pzeta(part);
-            pzeta += r*LineSegmentMapData_get_gauss_noise_matrix(el,5,5);
-            LocalParticle_update_pzeta(part,pzeta);
-        END_PER_PARTICLE_BLOCK;
-}
-
-GPUFUN
-void correlated_gaussian_noise(LocalParticle *part0,
-                    LineSegmentMapData el){
-
-    START_PER_PARTICLE_BLOCK(part0, part);
-        double x[6];
-        double r[6];
-        for(unsigned int i=0;i<6;++i){
-            x[i] = RandomNormal_generate(part);
-            r[i] = 0.0;
-        }
-        for(unsigned int i=0;i<6;++i){
-            for(unsigned int j=0;j<6;++j){
-                r[i] += LineSegmentMapData_get_gauss_noise_matrix(el,i,j)*x[j];
-            }
-        }
-        LocalParticle_add_to_x(part,r[0]);
-        LocalParticle_add_to_px(part,r[1]);
-        LocalParticle_add_to_y(part,r[2]);
-        LocalParticle_add_to_py(part,r[3]);
-        LocalParticle_add_to_zeta(part,r[4]);
-        double pzeta = LocalParticle_get_pzeta(part);
-        LocalParticle_update_pzeta(part,pzeta+r[5]);
-    END_PER_PARTICLE_BLOCK;
-}
-
-GPUFUN
-void LineSegmentMap_track_local_particle(LineSegmentMapData el, LocalParticle* part0){
-
-    remove_closed_orbit(part0,
-        LineSegmentMapData_get_x_ref(el, 0),
-        LineSegmentMapData_get_px_ref(el, 0),
-        LineSegmentMapData_get_y_ref(el, 0),
-        LineSegmentMapData_get_py_ref(el, 0));
-
-    remove_dispersion(part0,
-        LineSegmentMapData_get_dx(el, 0),
-        LineSegmentMapData_get_dpx(el, 0),
-        LineSegmentMapData_get_dy(el, 0),
-        LineSegmentMapData_get_dpy(el, 0));
-
-    transverse_motion(part0, el);
-
-    longitudinal_motion(part0, el);
-
-    energy_and_reference_increments(part0,
-        LineSegmentMapData_get_energy_increment(el),
-        LineSegmentMapData_get_energy_ref_increment(el));
-
-    if (LineSegmentMapData_get_uncorrelated_rad_damping(el) == 1){
-        uncorrelated_radiation_damping(part0,el);
-    }
-    
-    if (LineSegmentMapData_get_correlated_rad_damping(el) == 1){
-        correlated_radiation_damping(part0,el);
-    }
-
-    if (LineSegmentMapData_get_uncorrelated_gauss_noise(el) == 1){
-        uncorrelated_gaussian_noise(part0,el);
-    }
-    
-    if (LineSegmentMapData_get_correlated_gauss_noise(el) == 1){
-        correlated_gaussian_noise(part0,el);
-    }
-
-
-    add_dispersion(part0,
-        LineSegmentMapData_get_dx(el, 1),
-        LineSegmentMapData_get_dpx(el, 1),
-        LineSegmentMapData_get_dy(el, 1),
-        LineSegmentMapData_get_dpy(el, 1));
-
-    add_closed_orbit(part0,
-        LineSegmentMapData_get_x_ref(el, 1),
-        LineSegmentMapData_get_px_ref(el, 1),
-        LineSegmentMapData_get_y_ref(el, 1),
-        LineSegmentMapData_get_py_ref(el, 1));
-
-    double const length = LineSegmentMapData_get_length(el);
-    START_PER_PARTICLE_BLOCK(part0, part);
-        LocalParticle_add_to_s(part, length);
-    END_PER_PARTICLE_BLOCK;
-}
-
-#endif
diff --git a/xtrack/particles/rng_src/.__afsCCC4 b/xtrack/particles/rng_src/.__afsCCC4
deleted file mode 100644
index 86b34e7fb..000000000
--- a/xtrack/particles/rng_src/.__afsCCC4
+++ /dev/null
@@ -1,64 +0,0 @@
-// copyright ############################### //
-// This file is part of the Xtrack Package.  //
-// Copyright (c) CERN, 2024.                 //
-// ######################################### //
-
-#ifndef XTRACK_BASE_RNG_H
-#define XTRACK_BASE_RNG_H
-
-#include <headers/track.h>
-
-// Only include if compiled standalone, outside of Xsuite
-#ifndef START_PER_PARTICLE_BLOCK
-#include <stdint.h>
-#endif
-
-// Combined LCG-Thausworthe generator from (example 37-4):
-// https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-37-efficient-random-number-generation-and-application
-#define MASK 0xffffffffUL
-#define TAUSWORTHE(s,a,b,c,d) ((((s) &c) <<d) &MASK) ^ (((((s) <<a) &MASK)^(s)) >>b)
-#define LCG(s,A,C) ((((A*(s)) &MASK) + C) &MASK)
-
-GPUFUN
-uint32_t rng_get_int32 (uint32_t *s1, uint32_t *s2, uint32_t *s3, uint32_t *s4 ){
-  *s1 = TAUSWORTHE (*s1, 13, 19, 4294967294UL, 12);  // p1=2^31-1
-  *s2 = TAUSWORTHE (*s2, 2, 25, 4294967288UL, 4);    // p2=2^30-1
-  *s3 = TAUSWORTHE (*s3, 3, 11, 4294967280UL, 17);   // p3=2^28-1
-  *s4 = LCG(*s4, 1664525, 1013904223UL);             // p4=2^32
-
-  // Combined period is lcm(p1,p2,p3,p4) ~ 2^121
-  return ((*s1) ^ (*s2) ^ (*s3) ^ (*s4));
-}
-
-#ifndef TWO_TO_32
-#define TWO_TO_32 4294967296.0
-#endif
-
-GPUFUN
-double rng_get (uint32_t *s1, uint32_t *s2, uint32_t *s3, uint32_t *s4 ){
-
-  return rng_get_int32(s1, s2, s3, s4) / TWO_TO_32; // uniform in [0, 1) 1e10 resolution
-
-}
-
-GPUFUN
-void rng_set (uint32_t *s1, uint32_t *s2, uint32_t *s3, uint32_t *s4, uint32_t s ){
-  *s1 = LCG (s, 69069, 0);
-  if (*s1 < 2) *s1 += 2UL;
-  *s2 = LCG (*s1, 69069, 0);
-  if (*s2 < 8) *s2 += 8UL;
-  *s3 = LCG (*s2, 69069, 0);
-  if (*s3 < 16) *s3 += 16UL;
-  *s4 = LCG (*s3, 69069, 0);
-
-  /* "warm it up" */
-  rng_get (s1, s2, s3, s4);
-  rng_get (s1, s2, s3, s4);
-  rng_get (s1, s2, s3, s4);
-  rng_get (s1, s2, s3, s4);
-  rng_get (s1, s2, s3, s4);
-  rng_get (s1, s2, s3, s4);
-  return;
-}
-
-#endif /* XTRACK_BASE_RNG_H */

From 1150f36c199164150f97295f7185376d9a847d57 Mon Sep 17 00:00:00 2001
From: NetchAna <vaibhavi.gawas@cern.ch>
Date: Wed, 25 Mar 2026 14:55:13 +0100
Subject: [PATCH 4/4] Ignore AFS temp and generated C files

---
 .gitignore | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.gitignore b/.gitignore
index 292503b26..4495f13d0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,3 +16,5 @@ xtrack/prebuilt_kernels
 !xtrack/prebuilt_kernels/*.py
 checkpoint_restart.dat
 
+.__afs*
+*.c
