Improve backtracking for bend_face

Author: jgray-19

src/madl_dynmap.mad

@@ -1649,7 +1649,7 @@ function M.curex_fringe (elm, m, lw_) -- [FRINGE_TEAPOT]                     --
   if sdir*lw == 1 then
 -- print(lw == 1 and 'forward entry' or 'backward exit')
     yrotation (-e, m, 1)
-    bend_face ( h, m   )
+    bend_face ( h, m, 1)
     if m.frng ~= 0 then
       if ftst(m.frng, fringe.bend) then bend_fringe(elm, m, 1) end
       if ftst(m.frng, fringe.mult) then mult_fringe(elm, m, 1) end
@@ -1666,7 +1666,7 @@ function M.curex_fringe (elm, m, lw_) -- [FRINGE_TEAPOT]                     --
       if ftst(m.frng, fringe.mult) then mult_fringe(elm, m, -1) end
       if ftst(m.frng, fringe.bend) then bend_fringe(elm, m, -1) end
     end
-    bend_face ( h, m    )
+    bend_face ( h, m, -1)
     yrotation ( e, m, -1)
   end
 
@@ -1735,40 +1735,42 @@ end
 
 -- fringes helpers ------------------------------------------------------------o
 
-function bend_face (h, m) -- [NEWFACE]                                   -- checked
+function bend_face (h, m, lw) -- [NEWFACE]                                   -- checked
   if h == 0 then return end
   local el, eld, knl in m
   if abs(el) < minlen or abs(knl[1]) < minstr then return end
 
-  m.atdebug(h, m, 'bend_face:0')
+  m.atdebug(h, m, lw, 'bend_face:0')
 
-  local edir, tdir, beam in m
-  local   k0h = 0.5*h * knl[1]/(eld or el)*edir
+  local edir, sdir, beam in m
+  ! h is weighted by tdir, el is weighted by sdir, so to get tdir back, we need to multiply by sdir
+  local   k0h = 0.5*h * knl[1]/(eld or el) * sdir
   local   chg = beam.charge
   local _beta = 1/beam.beta
 
   for i=1,m.npar do
     local x, px, y, py, t, pt, beam in m[i]
     local _beta = beam and 1/beam.beta or _beta
     local chg   = beam and beam.charge or chg
+    local k0hq  = k0h*chg -- k0hq must be weighted by bdir
 
-    if tdir == 1 then -- to insure reversal symmetry ; horizontal wedge (only dir)
-      px = px + chg*k0h*x^2
+    if sdir == 1 then -- to insure reversal symmetry ; horizontal wedge (only sdir)
+      px = px + k0hq*x^2
     end
 
     local dpp      =  1 + 2*pt*_beta + pt^2
     local _pt2     =  1/(dpp - px^2)
-    local xi       =  2*k0h*tdir*sqrt(dpp)*_pt2
-    local dxi_px   =  2*px*xi           *_pt2
-    local dxi_ddel = -2   *xi*(1+pt)    *_pt2
+    local xi       =  2*k0hq*sqrt(dpp)*_pt2
+    local dxi_px   =  2*px*xi         *_pt2
+    local dxi_ddel = -2   *xi*(1+pt)  *_pt2
 
     x  = x / (1-dxi_px*y^2)
     px = px - xi*y^2
     py = py - 2*xi*x*y
     t  = t  - dxi_ddel*x*y^2
 
-    if tdir == -1 then -- to insure reversal symmetry; horizontal wedge (only dir)
-      px = px - chg*k0h*x^2
+    if sdir == -1 then -- to insure reversal symmetry; horizontal wedge (only sdir)
+      px = px + k0hq*x^2
     end
 
     m[i].x  = x
@@ -1777,7 +1779,7 @@ function bend_face (h, m) -- [NEWFACE]                                   -- chec
     m[i].t  = t
   end
 
-  m.atdebug(h, m, 'bend_face:1')
+  m.atdebug(h, m, lw, 'bend_face:1')
 end
 
 function mad8_wedge (e, m, lw_, wc) -- [FRINGE_TEAPOT:MAD8_WEDGE]            -- checked