Merge pull request #404 from jgray-19/fix_sbend_face

Backtracking Improvements

Author: ldeniau

src/mad_dynmap.cpp

@@ -1193,18 +1193,32 @@ inline void bend_face (cflw<M> &m, num_t lw, const V &h)
     if (m.sdir == 1) p.px += k0hq*sqr(p.x);
 
     T dpp      =  1 + 2/m.beta*p.pt + sqr(p.pt);
+    T y2       = sqr(p.y);
     T _pt2     =  1/(dpp - sqr(p.px));
     T xi       =  2*k0hq * sqrt(dpp)*_pt2;
     T dxi_px   =  2*p.px*xi         *_pt2;
-    T dxi_ddel = -2     *xi*(1+p.pt)*_pt2;
-    T y2       = sqr(p.y);
-
-    p.x  /= 1-dxi_px*y2;
-    p.px -= xi*y2;
-    p.py -= 2*xi*p.x*p.y;
-    p.t  += dxi_ddel*p.x*y2;
+    
+    if (m.sdir == -1) {
+      T npx = p.px - xi*y2;
+      _pt2   = 1/(dpp - sqr(npx));
+      xi     = 2*k0hq * sqrt(dpp)*_pt2;
+      dxi_px = 2*npx*xi          *_pt2;
+    }
 
-    if (m.sdir == -1) p.px += k0hq*sqr(p.x);
+    T dxi_ddel = -2     *xi*(1+p.pt)*_pt2;
+    if (m.sdir == 1){
+      p.x  /= 1-dxi_px*y2;
+      p.px -= xi*y2;
+      p.py -= 2*xi*p.x*p.y;
+      p.t  += dxi_ddel*p.x*y2;
+    } else {
+      p.t  += dxi_ddel*p.x*y2;
+      p.py -= 2*xi*p.x*p.y;
+
+      p.x  /= 1-dxi_px*y2;
+      p.px -= xi*y2;
+      p.px += k0hq*sqr(p.x);
+    }
   }
   mdump(1);
 }
@@ -1306,12 +1320,13 @@ inline void bend_fringe (cflw<M> &m, num_t lw)
     T fi2 = -2*co1*xyp*_yp2 - co3*2*xp*xyp    *pz;
     T fi3 =                 - co3*(1 + xp2*(1+yp2));
 
-    T kx = fi1*(1+xp2)*_pz  + fi2*xyp*_pz      - fi3*xp;
     T ky = fi1*xyp*_pz      + fi2*yp2*_pz      - fi3*yp;
-    T kz = fi1*tfac*xp*_pz2 + fi2*tfac*yp*_pz2 - fi3*tfac*_pz;
+    T y  = 2*p.y / (1 + sqrt(1-2*ky*p.y));
 
     if (m.sdir == 1) {
-      T y  = 2*p.y / (1 + sqrt(1-2*ky*p.y));
+      // Only calculate kx and kz once
+      T kx = fi1*(1+xp2)*_pz  + fi2*xyp*_pz      - fi3*xp;
+      T kz = fi1*tfac*xp*_pz2 + fi2*tfac*yp*_pz2 - fi3*tfac*_pz;
       T y2 = sqr(y);
 
       p.x  += 0.5*kx*y2;
@@ -1320,11 +1335,34 @@ inline void bend_fringe (cflw<M> &m, num_t lw)
       p.y   = y;
     } else { // need to reverse y-dependence
       T y2 = sqr(p.y);
+      // recalculate kx, ky and kz with a new py (Work out what it probably was)
+      T   npy = p.py + (4*c3*y2 + b0*tan(fi0))*y;
+      
+        pz = sqrt(dpp - sqr(p.px) - sqr(npy));
+       _pz = 1/pz;
+      _pz2 = sqr(_pz);
+      
+      xp  = p.px/pz,  yp  = npy/pz;
+      xyp = xp*yp  ,  yp2 = 1+sqr(yp);
+      xp2 = sqr(xp), _yp2 = 1/yp2;
+
+      fi0 = atan((xp*_yp2)) - c2*(1 + xp2*(1+yp2))*pz;
+      co2 = b0/sqr(cos(fi0));
+      co1 = co2/(1 + sqr(xp*_yp2))*_yp2;
+      co3 = co2*c2;
+
+      fi1 =    co1          - co3*2*xp*(1+yp2)*pz;
+      fi2 = -2*co1*xyp*_yp2 - co3*2*xp*xyp    *pz;
+      fi3 =                 - co3*(1 + xp2*(1+yp2));
+
+      T kx = fi1*(1+xp2)*_pz  + fi2*xyp*_pz      - fi3*xp;
+        ky = fi1*xyp*_pz      + fi2*yp2*_pz      - fi3*yp;
+      T kz = fi1*tfac*xp*_pz2 + fi2*tfac*yp*_pz2 - fi3*tfac*_pz;
 
       p.x  -= 0.5*kx*y2;
       p.py += (4*c3*y2 + b0*tan(fi0))*p.y;
       p.t  -= (0.5*kz  + c3*y2*sqr(relp)*tfac)*y2;
-      p.y  *= 0.5*(1 + sqrt(1-2*ky*p.y));
+      p.y  -= y2*ky/2;
     }
   }
   mdump(1);

src/madl_dynmap.mad

@@ -1777,21 +1777,33 @@ function bend_face (elm, m, lw, h) -- [NEWFACE]                              --
       px = px + k0hq*x^2
     end
 
-    -- This is still not fully revesible in 3rd and 4th order, due to self dependance of px
-    local dpp      =  1 + 2/beta*pt + pt^2
-    local _pt2     =  1/(dpp - px^2)
-    local xi       =  2*k0hq*sqrt(dpp)*_pt2
-    local dxi_px   =  2*px*xi         *_pt2
-    local dxi_ddel = -2   *xi*(1+pt)  *_pt2
-    local y2       = y^2
-
-    x  = x / (1-dxi_px*y2) -- this affects reversibility also, for py
-    px = px - xi*y2
-    py = py - 2*xi*x*y
-    t  = t  + dxi_ddel*x*y2
-
-    if sdir == -1 then -- to insure reversal symmetry; horizontal wedge (only sdir) -> PTC uses edir, and if we do this, reversing edir does not give the same result
-      px = px + k0hq*x^2
+    -- This is still not fully revesible in 4th order, due to dependancy of x on px and vice versa.
+    local dpp    = 1 + 2/beta*pt + pt^2
+    local y2     = y^2
+    local _pt2   = 1/(dpp - px^2)
+    local xi     = 2*k0hq*sqrt(dpp)*_pt2
+    local dxi_px = 2*px*xi         *_pt2
+
+    if sdir == -1 then 
+      local npx = px - xi*y2 
+      _pt2   = 1/(dpp - npx^2)
+      xi     = 2*k0hq*sqrt(dpp)*_pt2
+      dxi_px = 2*npx*xi        *_pt2
+    end
+    
+    local dxi_ddel = -2*xi*(1+pt)*_pt2
+
+    if sdir == 1 then 
+      x  = x / (1-dxi_px*y2)
+      px = px - xi*y2
+      py = py - 2*xi*x*y
+      t  = t  + dxi_ddel*x*y2
+    else
+      t = t + dxi_ddel*x*y2
+      py = py - 2*xi*x*y
+
+      x = x / (1-dxi_px*y2)
+      px = px - xi*y2 + k0hq*x^2
     end
 
     m[i].x  = x
@@ -1863,6 +1875,34 @@ function bend_wedge (elm, m, lw_, e) -- [WEDGE] see also [sr]bend_thick      --
   m.atdebug(elm, m, 'bend_wedge:1')
 end
 
+local function bend_fringe_param (dpp, px, py, c2, b0, tfac, only_ky)
+  -- This map could do with some cleaning, use of pz, _pz and _pz2 is inconsistent
+  local   pz = sqrt(dpp - px^2 - py^2)
+  local  _pz = 1/pz
+  local _pz2 = _pz^2
+
+  local xp,  yp  = px/pz, py/pz
+  local xyp, yp2 = xp*yp, 1+yp^2
+  local xp2,_yp2 = xp^2 , 1/yp2
+
+  local fi0 = atan((xp*_yp2)) - c2*(1 + xp2*(1+yp2))*pz
+  local co2 = b0/cos(fi0)^2
+  local co1 = co2/(1 + (xp*_yp2)^2)*_yp2
+  local co3 = co2*c2
+
+  local fi1 =    co1          - co3*2*xp*(1+yp2)*pz
+  local fi2 = -2*co1*xyp*_yp2 - co3*2*xp*xyp    *pz
+  local fi3 =                 - co3*(1 + xp2*(1+yp2))
+
+  local ky = fi1*xyp*_pz      + fi2*yp2*_pz      - fi3*yp
+  if only_ky then return fi0, nil, ky, nil end -- only ky is needed, speed up calculation
+
+  local kx = fi1*(1+xp2)*_pz  + fi2*xyp*_pz      - fi3*xp
+  local kz = fi1*tfac*xp*_pz2 + fi2*tfac*yp*_pz2 - fi3*tfac*_pz
+  return fi0, kx, ky, kz
+end
+
+
 function bend_fringe (elm, m, lw) -- [FRINGE_DIPOLE]                         -- checked
   if abs(m.knl[1]) < minang then return end
   m.atdebug(elm, m, 'bend_fringe:0')
@@ -1884,46 +1924,29 @@ function bend_fringe (elm, m, lw) -- [FRINGE_DIPOLE]                         --
     local   c2 = b0*fh*2
 
     local  dpp = 1 + 2/beta*pt + pt^2
-    local   pz = sqrt(dpp - px^2 - py^2)
-    local  _pz = 1/pz
-    local _pz2 = _pz^2
     local relp = invsqrt(dpp)
     local tfac = -(1/beta + pt)
     local   c3 = b0^2*fsad*relp
 
-    local xp,  yp  = px/pz, py/pz
-    local xyp, yp2 = xp*yp, 1+yp^2
-    local xp2,_yp2 = xp^2 , 1/yp2
-
-    local fi0 = atan((xp*_yp2)) - c2*(1 + xp2*(1+yp2))*pz
-    local co2 = b0/cos(fi0)^2
-    local co1 = co2/(1 + (xp*_yp2)^2)*_yp2
-    local co3 = co2*c2
-
-    local fi1 =    co1          - co3*2*xp*(1+yp2)*pz
-    local fi2 = -2*co1*xyp*_yp2 - co3*2*xp*xyp    *pz
-    local fi3 =                 - co3*(1 + xp2*(1+yp2))
-
-    local kx = fi1*(1+xp2)*_pz  + fi2*xyp*_pz      - fi3*xp
-    local ky = fi1*xyp*_pz      + fi2*yp2*_pz      - fi3*yp
-    local kz = fi1*tfac*xp*_pz2 + fi2*tfac*yp*_pz2 - fi3*tfac*_pz
-
+    local fi0, kx, ky, kz = bend_fringe_param(dpp, px, py, c2, b0, tfac, sdir==-1) -- this method is not perfect.
+    local ny = 2*y / (1 + sqrt(1-2*ky*y))
     if sdir == 1 then
-      y  = 2*y / (1 + sqrt(1-2*ky*y))
-      local y2 = y^2
+      local ny2 = ny^2
 
-      m[i].x  = x  + 0.5*kx*y2
-      m[i].py = py - (4*c3*y2 + b0*tan(fi0))*y
-      m[i].t  = t  + (0.5*kz  + c3*y2*relp^2*tfac)*y2
-      m[i].y  = y
+      m[i].x  = x  + 0.5*kx*ny2
+      m[i].py = py - (4*c3*ny2 + b0*tan(fi0))*ny
+      m[i].t  = t  + (0.5*kz  + c3*ny2*relp^2*tfac)*ny2
+      m[i].y  = ny
     else -- need to reverse y-dependence
-      -- This is still not fully revesible in 3rd and 4th order, due to self dependance of y
+      -- This is still not fully revesible in 5thand 6th order, due to self dependance of py -> results in max error of 1e-11 and 2e-9 for 5th and 6th order respectively
+      local npy = py + (4*c3*ny^2 + b0*tan(fi0))*ny
+      fi0, kx, ky, kz = bend_fringe_param(dpp, px, npy, c2, b0, tfac) -- this method is not perfect.
       local y2 = y^2
 
       m[i].x  = x  - 0.5*kx*y2
       m[i].py = py + (4*c3*y2 + b0*tan(fi0))*y
       m[i].t  = t  - (0.5*kz + c3*y2*relp^2*tfac)*y2
-      m[i].y  = 0.5*y * (1 + sqrt(1-2*ky*y))
+      m[i].y  = y - y2*ky/2
     end
   end
 

tests/tests/test-ng-maps/test-curved-maps.mad

@@ -114,10 +114,10 @@ function TestCurved:testBENDFACE()
     nslice = {1},
     energy = {1, 6500},  -- {1, 450, 6500}
 
-    tol = 100, 
+    tol = {5, 200, 200, 400, 4e3},--Non zero h2 causes failure at 4th order (h1 is fine for these tolerances)
 
-    h1 = {0,-0.04, 0.05},
-    h2 = {0, 0.04,-0.05},
+    h1 = {0,-0.04, 0.05}, -- This works with backtracking for tolerances above.
+    h2 = {0, 0.04,-0.05}, -- This does not work with backtracking at all
     alist = tblcat(ref_cfg.alist, {"h1", "h2"}),
     plot_info = {
       title    = "Bend Face NG v NG Maps",

tests/tests/test-ng-maps/test-fringe-maps.mad

@@ -69,7 +69,7 @@ function TestFringes:testQUADf() -- Test the straight fringes
     nslice = {1},
     energy = {1, 6500},  -- {1, 450, 6500}
 
-    tol = 1e3,
+    tol = 500, -- bend_fringe will pass backtracking at 480 eps
 
     k1     = {0, -0.15, 0.2},
     k1s    = {0, -0.15, 0.2},