From 3497dec1bdee035abc75af4e4ad6ca012ad689a9 Mon Sep 17 00:00:00 2001 From: OpenAI Date: Thu, 25 Jun 2026 07:35:41 +0000 Subject: [PATCH] Enable inverse tracking for dipole fringe edges --- tests/test_elements_thick.py | 58 ++++- .../beam_elements/elements_src/dipoleedge.h | 24 +- .../beam_elements/elements_src/magnet_edge.h | 1 + .../track_dipole_edge_nonlinear.h | 35 +++ .../elements_src/track_dipole_fringe.h | 233 ++++++++++++++++++ .../beam_elements/elements_src/track_magnet.h | 13 +- .../elements_src/track_magnet.template.h | 13 +- .../elements_src/track_magnet_edge.h | 90 ++++++- 8 files changed, 433 insertions(+), 34 deletions(-) diff --git a/tests/test_elements_thick.py b/tests/test_elements_thick.py index 2a61036fd..03ad1d53a 100644 --- a/tests/test_elements_thick.py +++ b/tests/test_elements_thick.py @@ -1138,6 +1138,37 @@ def test_fringe_implementations(test_context): xo.assert_allclose(np.linalg.det(R_ptc), 1, rtol=0, atol=1e-8) # Symplecticity check +@for_all_test_contexts +def test_full_dipole_edge_backtrack(test_context): + + skip_if_forbid_compile() + + fringe = xt.DipoleEdge( + k=0.12, fint=100, hgap=0.035, model='full') + + line = xt.Line(elements=[fringe]) + line.particle_ref = xp.Particles( + mass0=xp.PROTON_MASS_EV, beta0=0.5) + line.reset_s_at_end_turn = False + line.build_tracker(_context=test_context) + + p0 = line.build_particles(px=0.5, py=0.001, y=0.01, delta=0.1) + p_test = p0.copy(_context=test_context) + + line.track(p_test) + assert np.all(p_test.state == 1) + line.track(p_test, backtrack=True) + + p0.move(_context=xo.context_default) + p_test.move(_context=xo.context_default) + + assert np.all(p_test.state == 1) + for coordinate in ['x', 'px', 'y', 'py', 'zeta', 'delta', 's']: + xo.assert_allclose( + getattr(p_test, coordinate), getattr(p0, coordinate), + rtol=0, atol=1e-12) + + @for_all_test_contexts def test_backtrack_with_bend_quadrupole_and_cfm(test_context): @@ -1188,8 +1219,8 @@ def test_backtrack_with_bend_quadrupole_and_cfm(test_context): xo.assert_allclose(p2.zeta, p0.zeta, atol=1e-15, rtol=0) xo.assert_allclose(p2.delta, p0.delta, atol=1e-15, rtol=0) - # Same for dipole edge - de = xt.DipoleEdge(e1=0.1, k=3, fint=0.3) + # Same for a full dipole edge + de = xt.DipoleEdge(e1=0.1, k=3, fint=0.3, model='full') line = xt.Line(elements=[de]) line.particle_ref = xp.Particles(mass0=xp.PROTON_MASS_EV, beta0=0.5) line.reset_s_at_end_turn = False @@ -1198,16 +1229,27 @@ def test_backtrack_with_bend_quadrupole_and_cfm(test_context): zeta=0.05, delta=0.01) p1 = p0.copy(_context=test_context) line.track(p1) - p1.move(_context=xo.context_default) assert np.all(p1.state == 1) - line.configure_bend_model(edge='full') p2 = p1.copy(_context=test_context) line.track(p2, backtrack=True) - p2.move(_context=xo.context_default) - assert np.all(p2.state == -32) - # Same for combined function magnet - cfm = xt.Bend(length=1.0, k1=0.2, angle=0.1) + p0.move(_context=xo.context_default) + p2.move(_context=xo.context_default) + assert np.all(p2.state == 1) + xo.assert_allclose(p2.s, p0.s, atol=1e-14, rtol=0) + xo.assert_allclose(p2.x, p0.x, atol=1e-14, rtol=0) + xo.assert_allclose(p2.px, p0.px, atol=1e-14, rtol=0) + xo.assert_allclose(p2.y, p0.y, atol=1e-14, rtol=0) + xo.assert_allclose(p2.py, p0.py, atol=1e-14, rtol=0) + xo.assert_allclose(p2.zeta, p0.zeta, atol=1e-14, rtol=0) + xo.assert_allclose(p2.delta, p0.delta, atol=1e-14, rtol=0) + + # Same for a combined-function magnet with dipole-only full edges + cfm = xt.Bend( + length=1.0, k1=0.2, angle=0.1, + edge_entry_model='dipole-only', edge_exit_model='dipole-only', + edge_entry_fint=0.5, edge_exit_fint=0.5, + edge_entry_hgap=0.03, edge_exit_hgap=0.03) line = xt.Line(elements=[cfm]) line.particle_ref = xp.Particles(mass0=xp.PROTON_MASS_EV, beta0=0.5) line.reset_s_at_end_turn = False diff --git a/xtrack/beam_elements/elements_src/dipoleedge.h b/xtrack/beam_elements/elements_src/dipoleedge.h index 4b8362bbd..cdd9f1f34 100644 --- a/xtrack/beam_elements/elements_src/dipoleedge.h +++ b/xtrack/beam_elements/elements_src/dipoleedge.h @@ -52,22 +52,26 @@ void DipoleEdge_track_local_particle(DipoleEdgeData el, LocalParticle* part0){ } else if (model == 1){ - if (LocalParticle_check_track_flag(part0, XS_FLAG_BACKTRACK)){ - START_PER_PARTICLE_BLOCK(part0, part); - LocalParticle_kill_particle(part, -32); - END_PER_PARTICLE_BLOCK; - return; - } - double const e1 = DipoleEdgeData_get_e1(el); double const fint = DipoleEdgeData_get_fint(el); double const hgap = DipoleEdgeData_get_hgap(el); double const k = DipoleEdgeData_get_k(el); int64_t const side = DipoleEdgeData_get_side(el); + int64_t const backtrack = + LocalParticle_check_track_flag(part0, XS_FLAG_BACKTRACK); - START_PER_PARTICLE_BLOCK(part0, part); - DipoleEdgeNonLinear_single_particle(part, k, e1, fint, hgap, side); - END_PER_PARTICLE_BLOCK; + if (backtrack){ + START_PER_PARTICLE_BLOCK(part0, part); + DipoleEdgeNonLinear_single_particle_backtrack( + part, k, e1, fint, hgap, side); + END_PER_PARTICLE_BLOCK; + } + else{ + START_PER_PARTICLE_BLOCK(part0, part); + DipoleEdgeNonLinear_single_particle( + part, k, e1, fint, hgap, side); + END_PER_PARTICLE_BLOCK; + } } diff --git a/xtrack/beam_elements/elements_src/magnet_edge.h b/xtrack/beam_elements/elements_src/magnet_edge.h index b7e3ac2f1..6de0cc411 100644 --- a/xtrack/beam_elements/elements_src/magnet_edge.h +++ b/xtrack/beam_elements/elements_src/magnet_edge.h @@ -38,6 +38,7 @@ void MagnetEdge_track_local_particle(MagnetEdgeData el, LocalParticle* part0) part0, model, is_exit, + is_exit, // physical_is_exit half_gap, knorm, kskew, diff --git a/xtrack/beam_elements/elements_src/track_dipole_edge_nonlinear.h b/xtrack/beam_elements/elements_src/track_dipole_edge_nonlinear.h index ea09741c9..5b535007f 100644 --- a/xtrack/beam_elements/elements_src/track_dipole_edge_nonlinear.h +++ b/xtrack/beam_elements/elements_src/track_dipole_edge_nonlinear.h @@ -43,4 +43,39 @@ void DipoleEdgeNonLinear_single_particle(LocalParticle* part, } } + +GPUFUN +void DipoleEdgeNonLinear_single_particle_backtrack(LocalParticle* part, + double const k, double const e1, double const fint, double const hgap, + int64_t const side +){ + + double sin_, cos_, tan_; + if (fabs(e1) < 10e-10) { + sin_ = -999.0; cos_ = -999.0; tan_ = -999.0; + } + else{ + sin_ = sin(e1); cos_ = cos(e1); tan_ = tan(e1); + } + + if (side == 0){ // entry: invert wedge, fringe, rotation + if (sin_ > -99.){ + Wedge_single_particle(part, e1, k); + } + DipoleFringe_single_particle_backtrack(part, fint, hgap, k); + if (sin_ > -99.){ + YRotation_single_particle(part, sin_, cos_, tan_); + } + } + else if (side == 1){ // exit: invert rotation, fringe, wedge + if (sin_ > -99.){ + YRotation_single_particle(part, sin_, cos_, tan_); + } + DipoleFringe_single_particle_backtrack(part, fint, hgap, -k); + if (sin_ > -99.){ + Wedge_single_particle(part, e1, k); + } + } +} + #endif \ No newline at end of file diff --git a/xtrack/beam_elements/elements_src/track_dipole_fringe.h b/xtrack/beam_elements/elements_src/track_dipole_fringe.h index 388dbc174..dd3370b74 100644 --- a/xtrack/beam_elements/elements_src/track_dipole_fringe.h +++ b/xtrack/beam_elements/elements_src/track_dipole_fringe.h @@ -83,6 +83,227 @@ void DipoleFringe_single_particle( LocalParticle_set_zeta(part, new_zeta); } +/* + * Fixed-cost approximate inverse of the MAD-NG dipole-fringe map. + * + * The exact inverse reduces to a scalar equation for the entrance py. The + * quadratic predictor expands fi0 through py^2 and is followed by one + * analytic Newton correction evaluated with the exact fi0. There is no + * convergence loop. + */ + +typedef struct DipoleFringeInverseTerms_s { + double b0; + double c3; + double relp; + double tfac; + double fi0; + double kx; + double ky; + double kz; +} DipoleFringeInverseTerms; + +GPUFUN +int64_t DipoleFringe_inverse_terms( + const double px, + const double pt, + const double delta, + const double beta0, + const double chi, + const double fint, + const double hgap, + const double k0, + const double py, + DipoleFringeInverseTerms* terms +) { + const double fh = hgap * fint; + const double fsad = (fh > 10e-10) ? 1. / (72. * fh) : 0.; + const double b0 = k0 * chi; + const double dpp = POW2(1. + delta); + const double pz_sq = dpp - POW2(px) - POW2(py); + if (!(pz_sq > 0.) || !isfinite(pz_sq)) return 1; + + const double pz = sqrt(pz_sq); + const double _pz = 1. / pz; + const double relp = 1. / sqrt(dpp); + const double tfac = -(1. / beta0 + pt); + const double c2 = 2. * b0 * fh; + const double c3 = POW2(b0) * fsad * relp; + + const double xp = px * _pz; + const double yp = py * _pz; + const double xyp = xp * yp; + const double yp2 = 1. + POW2(yp); + const double xp2 = POW2(xp); + const double _yp2 = 1. / yp2; + + const double fi0 = atan(xp * _yp2) + - c2 * (1. + xp2 * (1. + yp2)) * pz; + const double cos2_fi0 = POW2(cos(fi0)); + if (!(cos2_fi0 > 0.) || !isfinite(cos2_fi0)) return 2; + + const double co2 = b0 / cos2_fi0; + const double co1 = co2 / (1. + POW2(xp * _yp2)) * _yp2; + const double co3 = co2 * c2; + + const double fi1 = co1 - 2. * co3 * xp * (1. + yp2) * pz; + const double fi2 = -2. * co1 * xyp * _yp2 + - 2. * co3 * xp * xyp * pz; + const double fi3 = -co3 * (1. + xp2 * (1. + yp2)); + + const double kx = fi1 * (1. + xp2) * _pz + + fi2 * xyp * _pz - fi3 * xp; + const double ky = fi1 * xyp * _pz + + fi2 * yp2 * _pz - fi3 * yp; + const double kz = fi1 * tfac * xp * POW2(_pz) + + fi2 * tfac * yp * POW2(_pz) - fi3 * tfac * _pz; + + if (!isfinite(b0) || !isfinite(c3) || !isfinite(relp) + || !isfinite(tfac) || !isfinite(fi0) || !isfinite(kx) + || !isfinite(ky) || !isfinite(kz)) return 3; + + terms->b0 = b0; + terms->c3 = c3; + terms->relp = relp; + terms->tfac = tfac; + terms->fi0 = fi0; + terms->kx = kx; + terms->ky = ky; + terms->kz = kz; + return 0; +} + +GPUFUN +int64_t DipoleFringe_single_particle_inverse_quadratic_corrected( + LocalParticle* part, + const double fint, + const double hgap, + const double k0 +) { + if (fabs(k0) < 10e-10) return 0; + + const double beta0 = LocalParticle_get_beta0(part); + if (beta0 == 0. || !isfinite(beta0)) return 1; + + const double exit_x = LocalParticle_get_x(part); + const double px = LocalParticle_get_px(part); + const double exit_y = LocalParticle_get_y(part); + const double exit_py = LocalParticle_get_py(part); + const double exit_zeta = LocalParticle_get_zeta(part); + const double pt = LocalParticle_get_ptau(part); + const double delta = LocalParticle_get_delta(part); + const double chi = LocalParticle_get_chi(part); + + const double fh = hgap * fint; + const double fsad = (fh > 10e-10) ? 1. / (72. * fh) : 0.; + const double b0 = k0 * chi; + const double c2 = 2. * b0 * fh; + const double dpp = POW2(1. + delta); + const double px2 = POW2(px); + const double a = dpp - px2; + if (!(a > 0.) || !isfinite(a)) return 2; + + const double s = sqrt(a); + const double relp = 1. / sqrt(dpp); + const double c3 = POW2(b0) * fsad * relp; + const double phi_c = atan(px / s) - c2 * (dpp + px2) / s; + const double phi_2 = -px / (2. * dpp * s) + + c2 * (dpp - 5. * px2) / (2. * POW3(s)); + const double tan_phi_c = tan(phi_c); + const double sec2_phi_c = 1. + POW2(tan_phi_c); + const double eta = exit_py + 4. * c3 * POW3(exit_y) + + b0 * tan_phi_c * exit_y; + const double alpha = b0 * exit_y * sec2_phi_c * phi_2; + double discriminant = 1. - 4. * alpha * eta; + + double discriminant_scale = fabs(4. * alpha * eta); + if (discriminant_scale < 1.) discriminant_scale = 1.; + const double discriminant_tolerance = 64. * DBL_EPSILON + * discriminant_scale; + if (!isfinite(tan_phi_c) || !isfinite(eta) || !isfinite(alpha) + || !isfinite(discriminant) + || discriminant < -discriminant_tolerance) return 3; + if (discriminant < 0.) discriminant = 0.; + + double entry_py = eta; + if (alpha != 0.) { + entry_py = 2. * eta / (1. + sqrt(discriminant)); + } + if (!isfinite(entry_py) || !(a - POW2(entry_py) > 0.)) return 3; + + DipoleFringeInverseTerms terms; + if (DipoleFringe_inverse_terms( + px, pt, delta, beta0, chi, fint, hgap, k0, entry_py, + &terms) != 0) return 4; + + /* One fixed analytic correction of the quadratic predictor. */ + const double pz = sqrt(a - POW2(entry_py)); + const double u = px * pz / a; + const double du_dpy = -px * entry_py / (a * pz); + const double dshape_dpy = -entry_py / pz + + 4. * px2 * entry_py / POW3(pz) + + 3. * px2 * POW3(entry_py) / (POW4(pz) * pz); + const double dfi0_dpy = du_dpy / (1. + POW2(u)) + - c2 * dshape_dpy; + const double tan_fi0 = tan(terms.fi0); + const double residual = entry_py + - 4. * terms.c3 * POW3(exit_y) + - terms.b0 * tan_fi0 * exit_y - exit_py; + const double residual_derivative = 1. + - terms.b0 * exit_y * (1. + POW2(tan_fi0)) * dfi0_dpy; + + if (isfinite(residual) && isfinite(residual_derivative) + && fabs(residual_derivative) > 1e-14) { + const double corrected_py = entry_py - residual / residual_derivative; + if (isfinite(corrected_py) && a - POW2(corrected_py) > 0.) { + DipoleFringeInverseTerms corrected_terms; + if (DipoleFringe_inverse_terms( + px, pt, delta, beta0, chi, fint, hgap, k0, + corrected_py, &corrected_terms) == 0) { + entry_py = corrected_py; + terms = corrected_terms; + } + } + } + + const double exit_y2 = POW2(exit_y); + const double entry_y = exit_y - 0.5 * terms.ky * exit_y2; + const double entry_x = exit_x - 0.5 * terms.kx * exit_y2; + const double entry_zeta = exit_zeta - beta0 * ( + 0.5 * terms.kz * exit_y2 + + terms.c3 * POW4(exit_y) * POW2(terms.relp) * terms.tfac); + + /* Select the same square-root branch as the direct map. */ + const double branch_value = 1. - terms.ky * exit_y; + double branch_scale = fabs(terms.ky * exit_y); + if (branch_scale < 1.) branch_scale = 1.; + if (branch_value < -64. * DBL_EPSILON * branch_scale) return 5; + + if (!isfinite(entry_x) || !isfinite(entry_y) + || !isfinite(entry_py) || !isfinite(entry_zeta)) return 6; + + LocalParticle_set_x(part, entry_x); + LocalParticle_set_y(part, entry_y); + LocalParticle_set_py(part, entry_py); + LocalParticle_set_zeta(part, entry_zeta); + return 0; +} + +GPUFUN +void DipoleFringe_single_particle_backtrack( + LocalParticle* part, + const double fint, + const double hgap, + const double k0 +) { + const int64_t status = + DipoleFringe_single_particle_inverse_quadratic_corrected( + part, fint, hgap, k0); + if (status != 0) { + LocalParticle_kill_particle(part, -32); + } +} + #endif // no XTRACK_FRINGE_FROM_PTC @@ -173,6 +394,18 @@ void DipoleFringe_single_particle( LocalParticle_set_py(part, new_py); LocalParticle_add_to_zeta(part, -d_tau * beta0); // PTC uses tau = ct } + +GPUFUN +void DipoleFringe_single_particle_backtrack( + LocalParticle* part, + const double fint, + const double hgap, + const double k0 +) { + if (fabs(k0) < 10e-10) return; + LocalParticle_kill_particle(part, -32); +} + #endif // XTRACK_FRINGE_FROM_PTC diff --git a/xtrack/beam_elements/elements_src/track_magnet.h b/xtrack/beam_elements/elements_src/track_magnet.h index 597809b1c..77cc2599c 100644 --- a/xtrack/beam_elements/elements_src/track_magnet.h +++ b/xtrack/beam_elements/elements_src/track_magnet.h @@ -410,6 +410,8 @@ void track_magnet_particles( double core_length, core_length_curved,factor_knl_ksl_body, factor_knl_ksl_edge,factor_backtrack_edge; + uint8_t edge_entry_is_exit = 0; + uint8_t edge_exit_is_exit = 1; // Backtracking if (LocalParticle_check_track_flag(part0, XS_FLAG_BACKTRACK)) { core_length = -length * weight; @@ -424,6 +426,11 @@ void track_magnet_particles( VSWAP(edge_entry_angle_fdown, edge_exit_angle_fdown); VSWAP(edge_entry_fint, edge_exit_fint); VSWAP(edge_entry_hgap, edge_exit_hgap); + { + const uint8_t tmp_is_exit = edge_entry_is_exit; + edge_entry_is_exit = edge_exit_is_exit; + edge_exit_is_exit = tmp_is_exit; + } VSWAP(theta_in, theta_out); VSWAP(cos_theta_in, cos_theta_out); VSWAP(sin_theta_in, sin_theta_out); @@ -489,7 +496,8 @@ void track_magnet_particles( track_magnet_edge_particles( part0, edge_entry_model, - 0, // is_exit + 0, // is_exit in the current tracking direction + edge_entry_is_exit, // physical_is_exit edge_entry_hgap, knorm, kskew, @@ -612,7 +620,8 @@ void track_magnet_particles( track_magnet_edge_particles( part0, edge_exit_model, - 1, // is_exit + 1, // is_exit in the current tracking direction + edge_exit_is_exit, // physical_is_exit edge_exit_hgap, knorm, kskew, diff --git a/xtrack/beam_elements/elements_src/track_magnet.template.h b/xtrack/beam_elements/elements_src/track_magnet.template.h index dc3766bb5..9ab0cc791 100644 --- a/xtrack/beam_elements/elements_src/track_magnet.template.h +++ b/xtrack/beam_elements/elements_src/track_magnet.template.h @@ -329,6 +329,8 @@ void track_magnet_particles( double core_length, core_length_curved,factor_knl_ksl_body, factor_knl_ksl_edge,factor_backtrack_edge; + uint8_t edge_entry_is_exit = 0; + uint8_t edge_exit_is_exit = 1; // Backtracking if (LocalParticle_check_track_flag(part0, XS_FLAG_BACKTRACK)) { core_length = -length * weight; @@ -343,6 +345,11 @@ void track_magnet_particles( VSWAP(edge_entry_angle_fdown, edge_exit_angle_fdown); VSWAP(edge_entry_fint, edge_exit_fint); VSWAP(edge_entry_hgap, edge_exit_hgap); + { + const uint8_t tmp_is_exit = edge_entry_is_exit; + edge_entry_is_exit = edge_exit_is_exit; + edge_exit_is_exit = tmp_is_exit; + } VSWAP(theta_in, theta_out); VSWAP(cos_theta_in, cos_theta_out); VSWAP(sin_theta_in, sin_theta_out); @@ -408,7 +415,8 @@ void track_magnet_particles( track_magnet_edge_particles( part0, edge_entry_model, - 0, // is_exit + 0, // is_exit in the current tracking direction + edge_entry_is_exit, // physical_is_exit edge_entry_hgap, knorm, kskew, @@ -531,7 +539,8 @@ void track_magnet_particles( track_magnet_edge_particles( part0, edge_exit_model, - 1, // is_exit + 1, // is_exit in the current tracking direction + edge_exit_is_exit, // physical_is_exit edge_exit_hgap, knorm, kskew, diff --git a/xtrack/beam_elements/elements_src/track_magnet_edge.h b/xtrack/beam_elements/elements_src/track_magnet_edge.h index be0a921b4..f31828518 100644 --- a/xtrack/beam_elements/elements_src/track_magnet_edge.h +++ b/xtrack/beam_elements/elements_src/track_magnet_edge.h @@ -13,11 +13,37 @@ #include "xtrack/beam_elements/elements_src/track_dipole_fringe.h" +GPUFUN +int64_t MagnetEdge_has_multipole_fringe( + const double* knorm, + const double* kskew, + const int64_t k_order, + GPUGLMEM const double* knl, + GPUGLMEM const double* ksl, + const double factor_knl_ksl, + const int64_t kl_order, + const double length +) { + for (int64_t ii = 1; ii <= k_order; ++ii) { + if (knorm[ii] != 0. || kskew[ii] != 0.) return 1; + } + + if (length != 0. && factor_knl_ksl != 0.) { + for (int64_t ii = 1; ii <= kl_order; ++ii) { + if (knl[ii] != 0. || ksl[ii] != 0.) return 1; + } + } + + return 0; +} + + 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 uint8_t is_exit, // Exit in the current tracking direction + const uint8_t physical_is_exit, // Exit in the forward tracking direction const double half_gap, const double* knorm, const double* kskew, @@ -78,12 +104,19 @@ void track_magnet_edge_particles( END_PER_PARTICLE_BLOCK; return; } - else if (model == 1 || model == 2) { // Full model - - if (factor_for_backtrack < 0) { + else if (model == 1 || model == 2) { // Full or dipole-only model + + const int64_t backtrack = factor_for_backtrack < 0.; + if (backtrack && model == 1 + && MagnetEdge_has_multipole_fringe( + knorm, kskew, k_order, knl, ksl, + factor_knl_ksl, kl_order, length)) { + /* The inverse implemented here covers the dipole fringe. The + * higher-order MAD-NG multipole fringe still has no inverse. */ START_PER_PARTICLE_BLOCK(part0, part); LocalParticle_kill_particle(part, -32); END_PER_PARTICLE_BLOCK; + return; } uint8_t should_rotate = 0; @@ -95,14 +128,21 @@ void track_magnet_edge_particles( tan_ = tan(face_angle); } - if (is_exit) k0 = -k0; + if (physical_is_exit) k0 = -k0; #define MAGNET_Y_ROTATE(PART) \ if (should_rotate) YRotation_single_particle((PART), -sin_, cos_, -tan_) + #define MAGNET_Y_ROTATE_BACKTRACK(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_DIPOLE_FRINGE_BACKTRACK(PART) \ + DipoleFringe_single_particle_backtrack( \ + (PART), fringe_integral, half_gap, k0) + #define MAGNET_MULTIPOLE_FRINGE(PART) \ MultFringe_track_single_particle( \ (PART), \ @@ -113,7 +153,7 @@ void track_magnet_edge_particles( ksl, \ kl_order, \ length / factor_knl_ksl, \ - is_exit, \ + physical_is_exit, \ /* min_order */ 1 \ ); // Above, I use the length to rescale knl and ksl. Here I am relying on @@ -124,41 +164,67 @@ void track_magnet_edge_particles( #define MAGNET_WEDGE(PART) \ if (should_rotate & (k_order >= 0)) Wedge_single_particle((PART), -face_angle, knorm[0]) + #define MAGNET_WEDGE_BACKTRACK(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 + #define MAGNET_QUAD_WEDGE_BACKTRACK(PART) \ + if (should_rotate & (k_order >= 1)) Quad_wedge_single_particle((PART), face_angle, knorm[1]) + + if (!backtrack && physical_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 + else if (!backtrack) { // 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; } + else if (physical_is_exit == 0) { // inverse of an entry edge + START_PER_PARTICLE_BLOCK(part0, part); + MAGNET_WEDGE_BACKTRACK(part); + if (model == 1){ + MAGNET_QUAD_WEDGE_BACKTRACK(part); + } + MAGNET_DIPOLE_FRINGE_BACKTRACK(part); + MAGNET_Y_ROTATE_BACKTRACK(part); + END_PER_PARTICLE_BLOCK; + } + else { // inverse of an exit edge + START_PER_PARTICLE_BLOCK(part0, part); + MAGNET_Y_ROTATE_BACKTRACK(part); + MAGNET_DIPOLE_FRINGE_BACKTRACK(part); + if (model == 1){ + MAGNET_QUAD_WEDGE_BACKTRACK(part); + } + MAGNET_WEDGE_BACKTRACK(part); + END_PER_PARTICLE_BLOCK; + } #undef MAGNET_Y_ROTATE + #undef MAGNET_Y_ROTATE_BACKTRACK #undef MAGNET_DIPOLE_FRINGE + #undef MAGNET_DIPOLE_FRINGE_BACKTRACK #undef MAGNET_MULTIPOLE_FRINGE #undef MAGNET_WEDGE + #undef MAGNET_WEDGE_BACKTRACK #undef MAGNET_QUAD_WEDGE + #undef MAGNET_QUAD_WEDGE_BACKTRACK } else if (model == 3) { // only ax ay cancellation (already done above) // do nothing