Merge pull request #285 from jgray-19/combined_sol

Combine curvature and solenoid multipole maps

Author: ldeniau

src/madl_dynmap.mad

@@ -499,75 +499,64 @@ function M.strex_kick (elm, m, lw, no_k0l) -- [KICKEX]                       --
   m.atdebug(elm, m, lw, 'strex_kick:1', no_k0l)
 end
 
-function M.strex_kickh (elm, m, lw) -- [KICKT]                               -- checked
-  if m.nmul == 0 then return end
+function M.strex_kickhs (elm, m, lw) -- [KICKT]                               -- unchecked
+  local lrad, el, nmul, ksi in m
+  if lrad == 0 or nmul == 0 and ksi == 0 then return end
 
-  m.atdebug(elm, m, lw, 'strex_kickh:0')
+  m.atdebug(elm, m, lw, 'strex_kickhs:0')
 
-  local el, tdir, nmul, knl, ksl, beam in m
+  local tdir, knl, ksl, beam in m
   local  wchg = lw*tdir*beam.charge
   local _beta = 1/beam.beta
+  local  hss  = (ksi^2/lrad) -- 4 H(s_0)^2/ds
 
   for i=1,m.npar do
     local x, px, y, py, t, pt, beam in m[i]
     local  wchg = beam and lw*tdir*beam.charge or wchg
     local _beta = beam and 1/beam.beta or _beta
     local    pz = sqrt(1 + (2*_beta)*pt + pt^2)
-    local bx,by = bxby(nmul, knl, ksl, x, y)
-
-    m[i].px = px - wchg*(by - knl[1]*pz)
-    m[i].py = py + wchg*(bx - ksl[1]*pz)
-    m[i].t  = t  - wchg*(knl[1]*x - ksl[1]*y)*(_beta+pt)/pz
-
-    -- field curvature
-    if el ~= 0 then
-      m[i].px = m[i].px - wchg*(knl[1]^2/el)*x 
-      m[i].py = m[i].py - wchg*(ksl[1]^2/el)*y
-    end
-  end
-
-  m.atdebug(elm, m, lw, 'strex_kickh:1')
-end
-
-function M.strex_kicks (elm, m, lw) -- [KICKT]                               -- unchecked
-  local lrad, nmul, ksi in m
-  if lrad == 0 or nmul == 0 and ksi == 0 then return end
-
-  m.atdebug(elm, m, lw, 'strex_kicks:0')
-
-  local tdir, knl, ksl, beam in m
-  local  wchg = lw*lrad*tdir*beam.charge
-  local _beta = 1/beam.beta
-  local  bsol = 0.5*ksi
-
-  for i=1,m.npar do
-    local x, px, y, py, t, pt, beam in m[i]
-    local  wchg = beam and lw*lrad*tdir*beam.charge or wchg
-    local _beta = beam and 1/beam.beta or _beta
 
     -- multipole
     local bx,by = bxby(nmul, knl, ksl, x, y)
 
-    m[i].px = px - wchg*by
+    m[i].px = px - wchg*by -- With 0 el and ksi ~= 0, we get problems with the multipole
     m[i].py = py + wchg*bx
 
-    -- solenoid
-    local _dpp = 1/sqrt(1 + (2*_beta)*pt + pt^2)
-    local  ang = bsol*_dpp
-    local ca, sa = cos(ang), sin(ang)
+    if abs(knl[1]) + abs(ksl[1]) > 0 then
+      m[i].px = m[i].px + (wchg*knl[1])*pz
+      m[i].py = m[i].py - (wchg*ksl[1])*pz
+      m[i].t  = t  - wchg*(knl[1]*x - ksl[1]*y)*(_beta+pt)/pz    
 
-    local nx  = ca* x + sa* y
-    local npx = ca*px + sa*py
-    local ny  = ca* y - sa* x
-    local npy = ca*py - sa*px
-    local nt  = t + ang*(_beta+pt)*(y*px - x*py)*_dp
+      -- field curvature
+      if lrad ~= 0 then
+        m[i].px = m[i].px - wchg*(knl[1]^2/lrad)*x 
+        m[i].py = m[i].py - wchg*(ksl[1]^2/lrad)*y
+      end
+    end
 
-    m[i].px = npx - 0.25*wchg*nx*_dpp
-    m[i].py = npy - 0.25*wchg*ny*_dpp
-    m[i].t  = nt  + 0.125*(_beta+pt)*wchg*(nx^2+ny^2)*_dpp^3
+    -- solenoid
+    if ksi ~= 0 and lrad ~= 0 then
+      local x, px, y, py, t, pt, beam in m[i]
+      local _dp = 1/(1 + (2*_beta)*pt + pt^2)
+      local _dpp = 1/pz
+      local  ang = (0.5*ksi*lw)*_dpp
+      local ca, sa = cos(ang), sin(ang)
+
+      local nx  = ca* x + sa* y
+      local npx = ca*px + sa*py
+      local ny  = ca* y - sa* x
+      local npy = ca*py - sa*px
+      local nt  = t - ang*(_beta+pt)*(y*px - x*py)*_dp
+
+      m[i].x  = nx
+      m[i].px = npx - 0.25*(wchg*hss)*nx*_dpp
+      m[i].y  = ny
+      m[i].py = npy - 0.25*(wchg*hss)*ny*_dpp
+      m[i].t  = nt  - 0.125*(_beta+pt)*(wchg*hss)*(nx^2+ny^2)*_dpp^3      
+    end
   end
 
-  m.atdebug(elm, m, lw, 'strex_kicks:1', no_k0l)
+  m.atdebug(elm, m, lw, 'strex_kickhs:1')
 end
 
 -- DKD [INTER_TEAPOT] ---------------------------------------------------------o

src/madl_etrck.mad

@@ -57,8 +57,8 @@ local thinonly, thickonly, driftonly, DKD, MKM                     in MAD.symint
 
 -- straight elements (DKD)
 local strex_drift, strex_kick             , strex_fringe           in MAD.dynmap
-local              strex_kickh --[[ h   ]]                         in MAD.dynmap
-local              strex_kicks --[[ ksi ]]                         in MAD.dynmap
+!local              strex_kickh --[[ h   ]]                         in MAD.dynmap
+local              strex_kickhs --[[ h, ksi ]]                     in MAD.dynmap
 
 -- curved elements (DKD)
 local curex_drift, curex_kick             , curex_fringe           in MAD.dynmap
@@ -247,20 +247,18 @@ local function track_multipole (elm, m)
   if m.nmul == 0 then return track_marker(elm, m) end
 
   local angle, ksi, lrad in elm
-  local knl, ksl, edir in m
+  local knl, ksl in m
 
-  m.el = lrad
+  m.lrad = lrad
 
   local kick = strex_kick
-  if abs(knl[1]) + abs(ksl[1]) >= minang then
-    kick = strex_kickh
-  elseif abs(ksi) >= minstr then
+  if abs(knl[1]) + abs(ksl[1]) + abs(ksi) >= minang then
     m.ksi = ksi
-    kick = strex_kicks
+    kick = strex_kickhs
   end
 
   trackelm(elm, m, thinonly, nil, kick, fnil)
-  m.ksi = nil
+  m.ksi, m.lrad = nil, nil
 end
 
 local function track_bbeam (elm, m)
@@ -548,34 +546,38 @@ local function track_solenoid (elm, m)
   get_mult(elm, m)
 
   local ds, tdir, edir in m
-  local ks, ksi in elm
+  local ks, ksi, lrad in elm
   local l = abs(ds)
+  local model  = elm.model or m.model
 
   if l < minlen then
-    errorf("invalid solenoid '%s' length=%.4e [m] (>0 expected)", elm.name, l)
+    if (model == "DKD" and lrad ~=0 ) then 
+      return track_multipole(elm, m)
+    else
+      errorf("invalid solenoid '%s' length=%.4e [m] (>0 expected)", elm.name, l)
+    end
   end
 
   m.el, m.eh = ds, 0
 
   local ksi = ksi + ks*l
   local no_ksi = abs(ksi) < minstr
-  local model  = elm.model  or m.model
   local method = elm.method or m.method
   local thick, kick, fringe
 
   if model == 'DKD' and no_ksi then
     thick, kick, fringe = strex_drift, strex_kick , strex_fringe
   elseif model == 'DKD' then
-    m.ksi = ksi
-    thick, kick, fringe = strex_drift, strex_kicks, strex_fringe
+    m.ksi, m.lrad = ksi, lrad
+    thick, kick, fringe = strex_drift, strex_kickhs, strex_fringe
   else
     m.ks = ksi/l*edir
-    thick, kick, fringe = solen_thick, strex_kick , solen_fringe
+    thick, kick, fringe = solen_thick, strex_kick , m.ptcmodel and curex_fringe or solen_fringe
   end
 
   local track  = #elm == 0 and trackelm or tracksub
-  track(elm, m, DKD[method], thick, strex_kick, fringe)
-  m.ks, m.ksi = nil, nil
+  track(elm, m, DKD[method], thick, kick, fringe)
+  m.ks, m.ksi, m.lrad = nil, nil, nil
 end
 
 -- eseptum