madl_track.mad

--[=[
 o-----------------------------------------------------------------------------o
 |
 | Track module
 |
 | Methodical Accelerator Design - Copyright CERN 2016+
 | Support: http://cern.ch/mad  - mad at cern.ch
 | Authors: L. Deniau, laurent.deniau at cern.ch
 | Contrib: I. Shreyber, irina.tecker at cern.ch
 |
 o-----------------------------------------------------------------------------o
 | You can redistribute this file and/or modify it under the terms of the GNU
 | General Public License GPLv3 (or later), as published by the Free Software
 | Foundation. This file is distributed in the hope that it will be useful, but
 | WITHOUT ANY WARRANTY OF ANY KIND. See http://gnu.org/licenses for details.
 o-----------------------------------------------------------------------------o

  Purpose:
  - TODO

 o-----------------------------------------------------------------------------o
]=]

-- help -----------------------------------------------------------------------o

local __help = {}
__help.track = [=[
NAME
  track -- track

SYNOPSIS
  TODO

DESCRIPTION
  The module track provides...

RETURN VALUES
  The track command

SEE ALSO
  element, sequence, survey
]=]

__help['track: references'] = [=[
[Wolski14]    A. Wolski, "Beam Dynamics in High Energy Particle Accelerators".
[Forest98]    E. Forest, "Beam Dynamics, A New Attitude and Framework".
[Forest02]    E. Forest, "Introduction to Polymorphic Tracking Code",
              KEK report 2002-3.
[Forest06]    E. Forest, "Geometric Integration for Particle Accelerators",
              J.Phys. A, Math.Gen. 39 (2006) 5321-5377.
[Ripken95]    K. Heinemann, G. Ripken, and F. Schmidt. "Construction of
              Nonlinear Symplectic Six-dimensional Thin-lens Maps by
              Exponentiation", Technical Report 95-189, DESY, 1995.
[SixTrack]    R. De. Maria et al., "SixTrack Physics Manual", CERN,
              http://cern.ch/SixTrack/doc/physics_manual/sixphys.pdf, Sep. 2016.
[Yoshida90]   H. Yoshida, "Construction of higher order symplectic integrators",
              Phys.Let. A, Volume 150, number 5,6,7.
[Burkhardt13] H. Burkhardt et al., "Improved Teapot Method and Tracking with
              Thick Quadrupoles for the LHC and its Upgrade", IPAC 2013.
[MADX-PTC]    MAD-X PTC Code, rev. 6099
]=]

__help['track: straight maps'] = [[
  [2] p.352, ch.12.1
      DKD: ch.12.1.3, p.353
      MKM: ch.12.1.4, p.354
  [3] DKD: p.72, Section K.4.2,
           PTC, Function INTER_STREX in file Sh_def_kind.f90)
      MKM: p.75, Section K.4.7,
           PTC, Function INTER_TKTF  in file Sh_def_kind.f90)
]]

__help['track: curved maps'] = [[
  [2] p.357, ch.12.2
  [3] PTC, p.76, Section K.4.9,Function INTER_TEAPOT in file Sh_def_kind.f90)
]]

__help['track: pushtkt'] = [[
[2] MKM: p.75, Section K.4.7,
[3] PTC, Function INTER_TKTF  in file Sh_def_kind.f90)
]]

__help['track: straight matrix-kick'] = [[
  [2] p.358, ch.12.2.2 Drift-Kick Split, drift in polar coordinates
  [3] p.76, PTC, Function SPROTR in file Sh_def_kind.f90)
]]

--[[TODO:
-Add charge, direction, mass(?)
-Deltap (recalculation of pt)
- KICKEX: how to read properly strength of thick element (i.e. convert
k1,k1s, k2 etc to knl,ksl)
-????KICKEX: add a check for thick element that there is only one strength,
i.e. quad only has k1 and not k0, k2 or???
]]
-- add patch, rotation + translation

--[[ METHOD:
2 - 2nd order, 1 kick per integration step
4 - Ruth-Neri-Yoshida 4th order, 3 kicks per step
6 - Yoshida 6th order, 7kicks per step
]]

--[[
  frame of integration:

 angle == 0 and k0 ~= 0 rbend on (??? rbend with the straight ref. frame ??)
 STRAIGHT REF. FRAME FOR RBEND : STREX

 angle == 0 and k0 == 0 rbend off (i.e. drift)
 DRIFT-TRACK

 angle ~= 0 and k0 ~= 0 sbend on
 SBEND, default RBEND:  TEAPOT

 angle ~= 0 and k0 == 0 sbend off( useless, what kind of element is this? broken sbend?)
       <=> rbend with l = l_arc and k0 = angle, what about errors?

 angle == 0 abd k0 == 0 but k[n] ~= 0: straight magnet
 DKD (model = 1 ): STREX ; MKM (model = 2): TKTF
]]

--[[
  Documentation....

  List of top level frame_track methods:
  --------------------------------------

  - patch_frame
  - drift_frame
    + l = 0 -> sub element only
  - thin_frame ignore l and angle but not patch and tilt
    + nmul == 0 and #elm > 0 -> drift_frame
  - thick_frames (all-in-one)
    + l = 0 -> thin_frame
    + nmul = 0 -> drift_frame
    + l > 0 and nmul > 0 and angle == 0 -> strait_frame
    + l > 0 and nmul > 0 and angle ~= 0 -> curved_frame

  - do we need curved and strait differenciation?
  - decompose entry and exit in more steps?
  - do we perform radiation and aperture check after each nst?
    + what about thin lens?

  - no more guards about lengths in thin, allow to store non-zero length
  - no more guards about lengths in thick, allow to be compatible with madx
  - check for angle and length to switch between integration schemes
]]

__help['track: classification of drift-kick models versus thick elements'] = [=[
              drift kind        kick kind        fringe kind
- angle=0
  rbend:      rbend_drift       rfcavity_kick    straight_fringe
  quadrupole: quadrupole_drift  quadrupole_kick  straight_fringe
  solenoid:   solenoid_drift    rfcavity_kick    solenoid_fringe
  rfcavity:   solenoid_drift    rfcavity_kick    rfcavity_fringe
  twcavity:   twcavity_drift    twcavity_kick    twcavity_fringe
- angle~=0
  sbend:      sbend_drift       polar_kick       sbend_fringe       sbend
  combined:   combined_drift    combined_kick    sbend_fringe       KTK, TKTF

- precedences (ordered):
  if ks ~=0   solenoid_drift    rfcavity_kick    solenoid_fringe
  if ang~=0   sbend_drift       polar_kick       sbend_fringe
  if k0 ~=0   rbend_drift       rfcavity_kick    straight_fringe
  if v  ~=0   rbend_drift       rfcavity_kick    rfcavity_fringe
  if k1 ~=0   quadrupole_drift  quadrupole_kick  straight_fringe
  default     straight_drift    straight_kick    straight_fringe

- fallbacks:
  straight_drift(l=0)      => (none)
  rbend_drift(k0=0)        => straight_drift
  solenoid_drift(ks=0)     => rbend_drift
  quadrupole_drift(k1=0)   => rbend_drift

  polar_drift(ang=0)       => rbend_drift (*)
  sbend_drift(k0=ang)      => polar_drift
  combined_drift(k1=0)     => sbend_drift

  straight_kick(nmul=0)    => (none)
  rfcavity_kick(v=0)       => straight_kick
  quadrupole_kick(k1=0)    => rfcavity_kick (*)

  polar_kick(ang=0)        => straight_kick (*)
  combined_kick(k1=0)      => polar_kick

  straight_fringe(nmul=0)  => (none)
  solenoid_fringe(ks=0)    => straight_fringe (*)
  rfcavity_fringe(v=0)     => straight_fringe (*)
  sbend_fringe(ang=0)      => straight_fringe (*)

  (*) fallbacks that should never occur according to precedences

- remarks:
  + fringe kind should be driven by element kind and not by precedence
  + quadrupole_drift = quadrupole matrix (linear) + drift correction (nonlinear)
    quadrupole_kick  = straight_kick - k1_kick
  + no solution for combined ang and/or k0 and ks ??
]=]

__help['track: drift-kick models and references'] = [=[
                   PTC KIND                         References
straight_drift:    DRIFT{R/P}
rbend_drift:       SPAR{R/P}
solenoid_drift:    KICK_SOL{R/P}
polar_drift:       SPROT{R/P}
sbend_drift:       SSEC{R/P}
quadrupole_kick:   PUSHTKT7{R/P}
combined_drift:    PUSHTKT7{R/P} ?                  KTK, TKTF

straight_kick:     KICKEX{R/P}
rfcavity_kick:     KICK_CAV{R/P}
polar_kick:        SKICK{R/P} (GETELECTRIC)
quadrupole_kick:   KICKPATH{R/P}+KICKTKT7{R/P}
combined_kick:     KICKPATH{R/P}+KICKTKT7{R/P}

twcavity(drift+kick): CAVE{R/P}_TRAV => RK{2,4,6}_CAV{R/P}+...

straight_fringe:   FRINGE_STRAIGHT{R/P}
sbend_fringe:      FRINGE_TEAPOT{R/P}
solenoid_fringe:   FRINGE_SOL{R/P}
rfcavity_fringe:   FRINGECAV{R/P}
twcavity_fringe:   FRINGE_CAV_TRAV{R/P}
]=]

-- locals ---------------------------------------------------------------------o

local vector, matrix, element, mtable, Command                  in MAD
local is_nil, is_number, is_sequence, is_beam                   in MAD.typeid
local minlen, minang, clight                                    in MAD.constant
local observe, ismoved                                          in element.flags

local abs, sqrt, max, sin, cos, tan, asin, acos, atan, atan2,
      sinh, cosh, sinc, fact, twopi                             in MAD.gmath

local minvolt, volt_c = 1e-6, 1e-3

-- forward declarations -------------------------------------------------------o

-- patches
local entry_patch, exit_patch                                     -- patches
local xrotation, yrotation, srotation                             -- rotations
local translate, changedir, changeref                             -- specials

-- drifts
local straight_drift , rbend_drift, quadrupole_drift              -- straight
local polar_drift    , sbend_drift, combined_drift                -- curved
local solenoid_drift                                              -- solenoid
local twcavity_drift                                              -- cavities

-- kicks
local thin_kick                                                   -- thin
local straight_kick  , quadrupole_kick                            -- straight
local polar_kick     , combined_kick                              -- curved
local rfcavity_kick  , twcavity_kick                              -- cavities

-- fringes
local straight_fringe                                             -- straight
local polar_fringe   , sbend_fringe                               -- curved
local solenoid_fringe                                             -- solenoid
local rfcavity_fringe, twcavity_fringe                            -- cavities

-- special --------------------------------------------------------------------o

local function empty_track () end

local function invalid_track (elm)
  error(string.format("invalid element %s of kind %s with length %s",
                      elm.name, elm.kind, elm.l))
end

element.  extra_element :set_functions { track = invalid_track,
                                     backtrack = invalid_track } -- not valid
element.special_element :set_functions { track =   empty_track,
                                     backtrack =   empty_track } -- default

-- patches --------------------------------------------------------------------o

-- X-rotation (pitch)

function xrotation (_, m, angle)   -- ROT_YZR <-> Rx(angle)
  if abs(angle) < minang then return end -- XZrot with x,px <-> y,py)

  local x, px, y, py, t, pt, sdir in m
  local sa, ca, ta = sin(angle*sdir), cos(angle*sdir), tan(angle*sdir)
  local _beta = 1/m.beam.beta
  local    pz = sqrt(1 + 2*_beta*pt + pt^2 - px^2 - py^2)
  local   _pz = 1/pz
  local   ptt = 1 - ta*py*_pz
  local  _ptt = 1/ptt

  m.y  = y/(ca*ptt)
  m.py = ca*py + sa*pz
  m.x  = x + ta*y*px*_pz*_ptt
  m.t  = t - ta*y   *_pz*_ptt*(_beta+pt)
end

-- Y-rotation (yaw)

function yrotation (_, m, angle)   -- ROT_XZR (Dragt's PROT) <-> Ry(angle)
  if abs(angle) < minang then return end

  local x, px, y, py, t, pt, sdir in m
  local sa, ca, ta = sin(angle*sdir), cos(angle*sdir), tan(angle*sdir)
  local _beta = 1/m.beam.beta
  local    pz = sqrt(1 + 2*_beta*pt + pt^2 - px^2 - py^2)
  local   _pz = 1/pz
  local   ptt = 1 - ta*px*_pz
  local  _ptt = 1/ptt

  m.x  = x/(ca*ptt)
  m.px = ca*px + sa*pz
  m.y  = y + ta*x*py*_pz*_ptt
  m.t  = t - ta*x   *_pz*_ptt*(_beta+pt)
end

-- S-rotation (roll, tilt)

function srotation (_, m, angle)   -- XYrot <-> Rz(-angle)?
  if abs(angle) < minang then return end

  local x, px, y, py, bdir in m
  local sa, ca = sin(angle*bdir), cos(angle*bdir)

  m.x  = ca*x  + sa*y
  m.y  = ca*y  - sa*x
  m.px = ca*px + sa*py
  m.py = ca*py - sa*px
end

-- Translation (dz treated as drift length)

function translate (_, m, dx, dy, dz)
  local bdir in m

  straight_drift(nil, m, dz)
  m.x = m.x - dx*bdir
  m.y = m.y - dy*bdir
end

-- Changedir (reverse x and angles)

function changedir(_, m)
  m.sdir = -m.sdir
  m.bdir = -m.bdir
end

-- Changeref (change reference frame)

function changeref (elm, m)
  local x, y, z, theta, phi, psi in elm
  xrotation(nil, m, phi   or 0)                    -- phi  : Elevation angle.
  yrotation(nil, m, theta or 0)                    -- theta: Azimuthal angle.
  srotation(nil, m, psi   or 0)                    -- psi  : Roll      angle.
  translate(nil, m, x or 0, y or 0, z or 0)
end

-- Entry/Exit patches

 -- not correct, need local survey patch
function entry_patch (elm, m)
  local theta, phi, psi, dtheta, dphi, dpsi, tilt in elm
  local x, y, z, dx, dy, dz in elm
  local mdir in m
  -- TODO: survey m.S (mdir = -1)
  xrotation(nil, m, (phi   or 0)+(dphi   or 0))    -- phi  : Elevation angle.
  yrotation(nil, m, (theta or 0)+(dtheta or 0))    -- theta: Azimuthal angle.
  srotation(nil, m, (psi   or 0)+(dpsi   or 0))    -- psi  : Roll      angle.
  translate(nil, m, (x or 0)+(dx or 0),
                    (y or 0)+(dy or 0),
                   ((z or 0)+(dz or 0))*mdir )
  srotation(nil, m, tilt)
end

function exit_patch (elm, m)
  local theta, phi, psi, dtheta, dphi, dpsi, tilt in elm
  local x, y, z, dx, dy, dz in elm
  local mdir in m
  -- TODO: survey m.S (mdir = 1)
  srotation(nil, m, -tilt)
  translate(nil, m, (x or 0)+(dx or 0),
                    (y or 0)+(dy or 0),
                   ((z or 0)+(dz or 0))*mdir )
  srotation(nil, m, (psi   or 0)+(dpsi   or 0))    -- psi  : Roll      angle.
  yrotation(nil, m, (theta or 0)+(dtheta or 0))    -- theta: Azimuthal angle.
  xrotation(nil, m, (phi   or 0)+(dphi   or 0))    -- phi  : Elevation angle.
end

-- strengths and phases -------------------------------------------------------o

-- Load element phases in mflow

local function get_phas (elm, m, nmul) -- l=0 (or ignored)
  local pnl, dpnl = elm.pnl or {}, elm.dpnl or {}
  local psl, dpsl = elm.psl or {}, elm.dpsl or {}
  for i=1,nmul do
    m.pnl[i] = (pnl[i] or 0) + (dpnl[i] or 0)
    m.psl[i] = (psl[i] or 0) + (dpsl[i] or 0)
  end
end

-- Load element multipoles and volts in mflow

local function get_mult (elm, m) -- l=0 (or ignored)
  local knl, dknl = elm.knl or {}, elm.dknl or {}
  local ksl, dksl = elm.ksl or {}, elm.dksl or {}
  local nmul = max(#knl, #ksl, #dknl, #dksl)

  m.ksi  = (elm.ksi  or 0) + (elm.dksi or 0) -- solenoid
  m.krf  = (elm.krf  or 0) + (elm.dkrf or 0) -- rf cavity (krf*z = omega_rf*t)
  m.vrf  = (elm.volt or 0) + (elm.dvrf or 0) -- rf volts
  m.ktot = nmul + abs(m.ksi) + abs(m.krf) + abs(m.vrf)

  if nmul >= 1 then
    for i=1,nmul do
      m.knl[i] = (knl[i] or 0) + (dknl[i] or 0)
      m.ksl[i] = (ksl[i] or 0) + (dksl[i] or 0)
    end
    if abs(m.krf) >= minang then
      get_phas(elm, m, nmul)
    end
    m.nmul = nmul
  else
    m.nmul, m.knl[1], m.ksl[1], m.knl[2] = 0, 0, 0, 0
  end
end

-- Load element strengths, multipoles and volts in mflow

local function get_kmult (elm, m) -- l~=0
  local l in elm
  if abs(l) < minlen then return get_mult(elm, m) end

  local k0 , k1 , k2 , k3  in elm
  local k0s, k1s, k2s, k3s in elm
  local nkn = k3  and 4 or k2  and 3 or k1  and 2 or k0  and 1 or 0
  local nks = k3s and 4 or k2s and 3 or k1s and 2 or k0s and 1 or 0
  local knl, dknl = elm.knl or {}, elm.dknl or {}
  local ksl, dksl = elm.ksl or {}, elm.dksl or {}
  local nmul = max(#knl, #ksl, #dknl, #dksl, nkn, nks)

  m.ksi  = (elm.ks or 0)*l + (elm.ksi  or 0) + (elm.dksi or 0)
  m.krf  =                   (elm.krf  or 0) + (elm.dkrf or 0)
  m.vrf  =                   (elm.volt or 0) + (elm.dvrf or 0)
  m.ktot = nmul + abs(m.ksi) + abs(m.krf) + abs(m.vrf)

  if nmul >= 1 then
    local kn = { k0 , k1 , k2 , k3  }
    local ks = { k0s, k1s, k2s, k3s }
    for i=1,nmul do
      m.knl[i] = (kn[i] or 0)*l + (knl[i] or 0) + (dknl[i] or 0)
      m.ksl[i] = (ks[i] or 0)*l + (ksl[i] or 0) + (dksl[i] or 0)
    end
    if abs(m.krf) >= minang then
      get_phas(elm, m, nmul)
    end
    m.nmul = nmul
  else
    m.nmul, m.knl[1], m.ksl[1], m.knl[2] = 0, 0, 0, 0
  end
end

-- Compute weighted multipoles

local function get_bxby (nmul, knl, ksl, x, y, lw)
  if nmul == 0 then return 0, 0 end

  local by = lw*knl[nmul] / fact(nmul-1)
  local bx = lw*ksl[nmul] / fact(nmul-1)
  local byt

  for i=nmul-1,1,-1 do
    byt = x*by - y*bx + lw*knl[i] / fact(i-1)
    bx  = y*by + x*bx + lw*ksl[i] / fact(i-1)
    by  = byt
  end
  return bx, by
end

-- drifts ---------------------------------------------------------------------o

__help['track: exact straight drift'] = [=[
  [Wolski14] p.86, ch. 3.1, eq. 3.13-3.18
  [Forest98] p.306, ch. 10.4.3.1, eq. 10.23a-10.23c
  [MADX-PTC] Function DRIFTR in file Sh_def_kind.f90
]=]

function straight_drift (_, m, l) -- exact, l~=0
  if abs(l) < minlen then return end

  local x, px, y, py, t, pt, T in m
  local _beta = 1/m.beam.beta
  local l_pz = l/sqrt(1 + 2*_beta*pt + pt^2 - px^2 - py^2)

  m.x = x + px*l_pz
  m.y = y + py*l_pz
  m.t = t - (_beta+pt)*l_pz + (1-T)*l*_beta
end

__help['track: exact rbend drift'] = [=[
  [Forest98] p.367, ch. 12.3.1, eq. 12.39a,12.39c, 12.39f
  [MADX-PTC] Function SPARR in file Sh_def_kind.f90
]=]

function rbend_drift (elm, m, l) -- exact, l~=0
  local knl in m
  if abs(knl[1]) < minang then
    return straight_drift(nil,m,l)
  end

  local x, px, y, py, t, pt, bdir, T in m
  local _beta = 1/m.beam.beta
  local lw = l/elm.l
  local k0 = bdir*lw*knl[1]

  local pz = sqrt(1 + 2*_beta*pt + pt^2 - px^2 - py^2)
  local pt = sqrt(1 + 2*_beta*pt + pt^2 - py^2)
  local xp = px - k0
  local ps = sqrt(1 + 2*_beta*pt + pt^2 - xp^2 - py^2)
  local xs = (asin(px/pt) - asin(xp/pt))/k0

  m.x  = x + (ps-pz)/k0
  m.px = xp
  m.y  = y + py*xs
  m.t  = t - (_beta+pt)*xs + (1-T)*l*_beta
end

__help['track: exact solenoid drift'] = [=[
  [Forest98] p.377, ch. 12.5.2, eq. 12.50 and solution
  [MADX-PTC] Function KICK_SOL in file Sh_def_kind.f90
]=]

function solenoid_drift (elm, m, l) -- exact, l~=0          TO BE REVIEWED
  if abs(m.ksi) < minang then
    return rbend_drift(elm,m,l)
  end

  local x, px, y, py, t, pt, T, bdir in m
  local _beta = 1/m.beam.beta
  local _el = 1/elm.l

  local bsol = m.ksi*bdir*_el/2
  local xp = px + bsol*y
  local yp = py - bsol*x

  local l_pz  = l/sqrt(1 + 2*_beta*pt + pt^2 - xp^2 - yp^2)
  local angle = l_pz*bsol

  local ca, sa, sc = cos(angle), sin(angle), sinc(angle)

  local xt  = ca*x  + l_pz*sc*px
  local pxt = ca*px - l_pz*sc*x *bsol^2
  local yt  = ca*y  + l_pz*sc*py
  local pyt = ca*py - l_pz*sc*y *bsol^2

  m.x  = ca*xt  + sa*yt
  m.px = ca*pxt + sa*pyt
  m.y  = ca*yt  - sa*xt
  m.py = ca*pyt - sa*pxt

  m.t  = t - (_beta+pt)*l_pz + (1-T)*l*_beta
end

__help['track: drift to correct pseudo-exact combined kick'] = [=[
  [MADX-PTC] Function KICKPATHR in file Sh_def_kind.f90
]=]

local function combadj_drift (_, m, l) -- l~=0, adjust combined_drift
  local x, px, y, py, t, pt, T in m
  local _beta = 1/m.beam.beta

  local l_pz = l/sqrt(1 + 2*_beta*pt + pt^2 - px^2 - py^2)
  m.x  = x + px*(l_pz-l)
  m.y  = y + py*(l_pz-l)
  m.t  = t - (_beta+pt)*l_pz + (1-T)*l*_beta
end

__help['track: pseudo-exact quad drift'] = [=[
  [MADX-PTC] Function PUSHTKT7 in file Sh_def_kind.f90
]=]

function quadrupole_drift (elm, m, l) -- comb_drift with angle == 0
  local _el   = 1/elm.l
  local _beta = 1/m.beam.beta
  local l_bet2_gam2 = l/(m.beam.gamma^2-1)

  local knl, bdir in m
  local k0, k1 = bdir*knl[1]*_el, bdir*knl[2]*_el
  local w2 = abs(k1)
  local w  = sqrt(w2)

  -- TODO: case w=0...
  local cw, sw = cos (w*l), sin (w*l)/w
  local ch, sh = cosh(w*l), sinh(w*l)/w

  if k1 < 0 then
    cw, sw, ch, sh = ch, sh, cw, sw
  end

  local npos in m
  if is_nil(npos) or npos == 'last' then
    local dl = is_nil(pos) and l/2 or l
    combadj_drift(_, m, dl)
  end

  local x, px, y, py, t, pt in m
  local m16 = (1-cw)/w2

  m.x  =    cw*x + sw*px -k0*m16
  m.px = -w*sw*x + cw*px -k0*sw
  m.y  =    ch*y + sh*py
  m.py =  w*sh*y + ch*py
  m.t  = t + l_bet2_gam2*pt

  if is_nil(npos) or npos == 'first' then
    local dl = is_nil(pos) and l/2 or l
    combadj_drift(_, m, l)
  end
  m.npos = nil
end

__help['track: pseudo-exact combined drift'] = [=[
  [MADX-PTC] Function PUSHTKT7 in file Sh_def_kind.f90
]=]

function combined_drift (elm, m, l)
  local _el   = 1/elm.l
  local _beta = 1/m.beam.beta
  local l_bet2_gam2 = l/(m.beam.gamma^2-1)

  local angle = elm.angle
  local lw    = l*_el
  local h     = lw*angle
  local hbar  = h*_beta

  local knl, bdir in m
  local k0, k1 = bdir*knl[1]*_el, bdir*knl[2]*_el
  local hk0pk1 = h*k0+k1
  local wx2, wy2 = abs(hk0pk1), abs(k1)
  local wx , wy  = sqrt(wx2), sqrt(wy2)

  local cx, sx, cy, sy
  if hk0pk1 < 0
  then cx, sx = cosh(wx*l), sinh(wx*l)/wx -- check wx == 0
  else cx, sx = cos (wx*l), sin (wx*l)/wx -- check wx == 0
  end
  if k1 < 0
  then cy, sy = cosh(wy*l), sinh(wy*l)/wy -- check wy == 0
  else cy, sy = cos (wy*l), sin (wy*l)/wy -- check wy == 0
  end

  local npos in m
  if is_nil(npos) or npos == 'last' then
    local dl = is_nil(pos) and l/2 or l
    combadj_drift(_, m, dl)
  end

  local x, px, y, py, t, pt in m
  local m16, m56 = (1-cx)/wx2, hbar^2*(l-sx)/wx2

  m.x  =      cx*x + sx*px + hbar*m16*pt + (h-k0)*m16
  m.px = -wx2*sx*x + cx*px + hbar*sx*pt  + (h-k0)*sx
  m.y  =      cy*y + sy*py
  m.py =  wy2*sy*y + cy*py
  m.t  = t - hbar*(sx*x+m16*px) + (l_bet2_gam2-m56)*pt

  if is_nil(npos) or npos == 'first' then
    local dl = is_nil(pos) and l/2 or l
    combadj_drift(_, m, l)
  end
  m.npos = nil
end

__help['track: exact polar drift'] = [=[
  [Forest02] p.76, ch. K.4.9, (general explaination)
  [Forest06] p.5365, eq. 127, exact drift -- WRONG: ref is an yrotation
  [MADX-PTC] Function SPROTR in file Sh_def_kind.f90
]=]

function polar_drift (elm, m, l) -- exact, l~=0, k0=angle
  local rho = elm.l/(elm.angle*m.sdir)
  local ang = l/rho
  local sa, ca, ta, sa2 = sin(ang), cos(ang), tan(ang), sin(ang/2)

  local x, px, y, py, t, pt, T in m
  local _beta = 1/m.beam.beta
  local pz  = sqrt(1 + 2*_beta*pt + pt^2 - px^2 - py^2)
  local _pz = 1/pz
  local ptt = 1 - ta*px*_pz
  local xr  = (x+rho)*ta*_pz

  m.x  = (x + rho*(2*sa2^2 + sa*px*_pz))/(ca*ptt)
  m.px = ca*px + sa*pz
  m.y  = y + xr*py/ptt
  m.t  = t - xr*(_beta+pt)/ptt + (1-T)*l*_beta
end

__help['track: exact sbend drift'] = [=[
  [Forest06] p.5365, eq. 126, exact drift
]=]

function sbend_drift (elm, m, l) -- exact, l~=0, k0~=angle
  local angle in elm
  local knl in m
  if abs(angle - knl[1]) < minang then
    return polar_drift(elm,m,l)
  end
  error("NYI")
end

__help['track: pseudo-exact combined drift'] = [=[
  [MADX-PTC] Function PUSHTKT7 in file Sh_def_kind.f90
]=]

-- PRECEDENCE SHOULD BE OUTSIDE INTEGRATORS
--local function combined_drift (elm, m, l) -- angle, k0, k1 (tunes)
--  if abs(m.knl[2]) < minang then -- Put these checks outside the integrator
--    if abs(angle) < minang then
--      return rbend_drift(elm,m,l)
--    else
--      return sbend_drift(elm,m,l)
--    end
--  else
--    if abs(angle) < minang then
--      return quadrupole_drift(elm,m,l)
--    else
--      return combined_drift  (elm,m,l)
--    end
--  end
--end

-- kicks ----------------------------------------------------------------------o

__help['track: generic thin kick'] = [[
  [MADX-PTC] Function KICKTR in file Sh_def_kind.f90
]]

-- Special case TO BE REVIEWED !!!
function thin_kick (elm, m) -- l=0 (ignored)
  local nmul in m
  if nmul == 0 then return end

--  local knlt, kslt = elm.knl or {}, elm.ksl or {} -- BUG: not retrieved by getmul
--  local knl1, ksl1 = (knlt[1] or 0), (kslt[1] or 0)

  local x, px, y, py, bdir, knl, ksl in m
  local bx, by = get_bxby(nmul, knl, ksl, x, y, 1)

  m.px = px - bdir*(by + knl[1])
  m.py = py + bdir*(bx - ksl[1])

  if knl1 ~= 0 or ksl1 ~= 0 then
    local t, pt in m
    local _beta = 1/m.beam.beta
    local pz = sqrt(1 + 2*_beta*pt + pt^2)

    local l = elm.lrad or 0 -- TODO: check if lrad or l?
    if l ~= 0 then -- dipole focusing and deltap
      m.px = m.px - knl[1]^2*x/l + bdir*knl[1]*(pz-1)
      m.py = m.py - ksl[1]^2*y/l + bdir*ksl[1]*(pz-1)
    end
--    m.t = t + bdir*(knl1*x - ksl1*y) * (_beta+pt)/pz
    m.t = t - bdir*(knl[1]*x - ksl[1]*y) * (_beta+pt)/pz
  end
end

__help['track: straight kick'] = [[
  [Forest98] p.354, ch. 12.1.3, eq. 12.6
  [MADX-PTC] Function KICKEXR in file Sh_def_kind.f90
  [SixTrack] p.12, ch. 3.5 (thin multipole), eq. 91-94
  [Ripken95] p.26, ch. 4.6.2, eq. 4.23 (and solution)
]]

function straight_kick (elm, m, l) -- l~=0
  local nmul in m
  if nmul == 0 then return end

  local x, px, y, py, knl, ksl, bdir in m
  local lw = abs(elm.l)>minlen and l/elm.l or 1
  local bx, by = get_bxby(nmul, knl, ksl, x, y, lw)

  m.px = px - bdir*by
  m.py = py + bdir*bx
end

__help['track: exact quad kick'] = [=[
  [Forest98] p.354, ch. 12.1.3, eq. 12.6
  [MADX-PTC] Function KICKTKT7R in file Sh_def_kind.f90
]=]

function quadrupole_kick (elm, m, l)
  if abs(m.knl[2]) < minang then -- should never occur due to precedences
    return rfcavity_kick(elm,m,l)
  end

  local x, px, y, py, t, pt, bdir, T, nmul, knl, ksl in m
  local _beta = 1/m.beam.beta
  local pz = sqrt(1 + 2*_beta*pt + pt^2) - 1

  local lw = l/elm.l
  local bx, by = get_bxby(nmul, knl, ksl, x, y, lw)

   m.px = px + lw*knl[1]*(pz-pt)*_beta
   m.t  = t  - lw*knl[1]*((pt+_beta)/(1+pz) - _beta)*x

   m.px = m.px - bdir*(by - lw*(knl[1] - knl[2]*x))
   m.py = m.py + bdir*(bx - lw*          knl[2]*y )
end

__help['track: polar kick'] = [[
  [Forest98] p. 358-361, ch. 12.2.2, eq. 12.18a-12.18f
  [Forest02] p. 76
  [MADX-PTC] Function SKICKR in file Sh_def_kind.f90
  [Ripken95] p. 22, ch. 4.3 (bending magnet) eq. 4.12 (and solution)
  [SixTrack] Thin Dipole, ch. 3.2.2, eq. 48-50
]]

function polar_kick (elm, m, l) -- [PTC] SKICK
  local angle in elm
  if abs(angle) < minang then -- should never occur due to precedences
    return straight_kick(elm,m,l)
  end

  local x, px, y, py, nmul, knl, ksl, bdir in m
  local lw = l/elm.l
  local bx, by = get_bxby(nmul, knl, ksl, x, y, lw) -- TODO: GETELECTRIC

  local h = lw*angle
  m.px = px - bdir*by*(1 + h*x)
  m.py = py + bdir*bx*(1 + h*x)
end

__help['track: pseudo-exact combined kick'] = [=[
Ref??
]=]

function combined_kick (elm, m, l) -- angle, k0, k1 (tunes)
--  local knl in m
--  if abs(knl[2]) < minang then
    return polar_kick(elm,m,l)
--  end
  -- error("NYI") -- TODO
end

__help['track: exact cavity kick'] = [=[
[MADX-PTC] Function KICCAVR in file Sh_def_kind.f90
]=]

function rfcavity_kick (elm, m, l) --RF CAV, CRAB, RFMULT with l>0 kick (kickcavr)
  --local vrf in m
  --if abs(vrf) < minvolt then
    return straight_kick(elm,m,l)
  --end

-- missing input: tt, phase, phase0, ff, ph, a, r, nbessel, freq, volt, lag

--  local x, px, y, py, t, pt, bdir, nmul, knl, ksl in m
--  local volt, freq, lag, nbessel in elm
--
--  if freq == 0 then
--    error("Frequency number has to be defined for RF Cavity (i.e. harmon NYI) "
--          .. elm.name)
--  end
--
--  local pc = m.beam.pc
--  local omega = twopi*freq / clight
--  local vl    = bdir*l*volt*volt_c/pc
--  local df, f, r2 = 0, 1, 1
--
--  do i=1,nbessel
--    r2  = -r2*(ko*omega)^2/(4*(i+1)^2)
--    dr2 =  r2*i
--    df  =  df + 2*dr2
--    r2  =  r2*(x^2+y^2)
--    f   =  f + r2
--  end
--
--  local ff, ph = elm.f or {}, elm.ph or {}
--  local a , r  = elm.a or  0, elm.r or 0
--  local arg = omega*(t + tt) + phasе + ph + phase0
--
--  if nbessel > 0 then
--    m.px = px-x*ff[ko]* df * vl * cos(arg)/(omega)
--    m.py = py-y*ff[ko]* df * vl * cos(arg)/(omega)
--  end
--
--  m.pt = m.pt -ff*f*vl*sin(arg)
--
--  local lw = elm.l>0 and l/elm.l or 1
--  local bx, by = get_bxby(nmul, knl, ksl, x, y, lw)
--
--  m.px = m.px - ff[ko]*           bdir*by/pc*(a + r*cos(arg))
--  m.py = m.py + ff[ko]*           bdir*bx/pc*(a + r*cos(arg))
--  m.pt = m.pt + ff[ko]*(ko*omega)*bdir*by/pc*     r*sin(arg)

--[[ useful???
  by = -lw*knl[nmul] / fact(nmul-1) / nmul
  bx = -lw*ksl[nmul] / fact(nmul-1) / nmul

  for i=nmul-1,1,-1 do
    byt = x*by - y*bx - lw*knl[i] / fact(i-1) / i
    bx  = y*by + x*bx - lw*ksl[i] / fact(i-1) / i
    by  = byt
  end

  byt = x*by - y*bx
  bx  = y*by + x*bx
  by  = byt
]]
end

-- fringes --------------------------------------------------------------------o

__help['track: straight fringe'] = [=[
  [Forest98] p. 389, ch. 13.2.2, eq. 13.31a-13.31f
]=]

local function mult_fringe_field (elm, m, pos)
	if elm.kill_fringe                        or
	   elm.kill_ent_fringe and pos == 'entry' or
	   elm.kill_exi_fringe and pos == 'exit'  then return end

  local x, px, y, py, bdir, nmul, knl, ksl in m
  local _beta = 1/m.beam.beta
  local _l = 1/elm.l

  local rx, ix, fx, fy, fxx, fxy, fyx, fyy = 1, 0, 0, 0, 0, 0, 0, 0

	for i=1, nmul do
	  local drx, dix = rx, ix
    rx  = drx*x - dix*y
    ix  = drx*y + dix*x

    local ni, nf = -bdir/(4*(i+1)), i+2/i
    local lkn, lks = knl[i]*_l, ksl[i]*_l

    local u   =  ni*(lkn*rx  - lks*ix )
    local v   =  ni*(lkn*ix  + lks*rx )
    local du  =  ni*(lkn*drx - lks*dix)
    local dv  =  ni*(lkn*dix + lks*drx)
	  local dux =   i*du
    local dvx =   i*dv
    local duy =  -i*dv
    local dvy =   i*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*(v + dvy*y)
    fyx = fyx + dux*y - nf*(v - dvx*x)
	end

  local _pz = 1/sqrt(1 + 2*pt*_beta + pt^2)
  local a = 1 - fxx*_pz
  local d = 1 - fyy*_pz
  local b =   - fyx*_pz
  local c =   - fxy*_pz
  local _det = 1/(a*d - b*c)

  m.x  = x - fx*_pz
  m.y  = y - fy*_pz
  m.px = (d*px - b*py)*_det
  m.py = (a*py - c*px)*_det
  m.t  = t + (_beta+pt)*(m.px*fx + m.py*fy)*_pz^3
end

local function straight_fringe_field (elm, m, pos)
  if elm.kill_fringe                        or
     elm.kill_ent_fringe and pos == 'entry' or
     elm.kill_exi_fringe and pos == 'exit'  then return end

-- Lee-Whiting formula, E. Forest ch 13.2.3, eq 13.33
	local x, px, y, py, t, pt, T, bdir in m
  local k0   = elm.k0 or 0
  local _beta = 1/m.beam.beta
-- soft edge
  local fint = elm.fint or 0  -- fringe field integral at entrance and exit of the bend
  local hgap = elm.hgap or 0  -- half gap of the magnet

  local _pz  = 1/sqrt(1 + 2*_beta*pt + pt^2 )
  local time_fac = (_beta + pt)*_pz

  local a2, b2 = elm.ksl[2]/elm.l or 0, elm.knl[2]/elm.l or 0
  local b  = sqrt(b2^2 + a2^2)

  local f1 = -bdir*fint*abs(fint)*b*_pz/24.0
  local f2 =  hgap*b*_pz

  local ang=-atan2(a2,b2)/2

  srotation(nil, m, -ang or 0)

  m.t = t - time_fac*(f1*x + f2*(1+f1/2)*px*_pz*exp(-f1))*px*_pz
          + time_fac*(f1*y + f2*(1-f1/2)*py*_pz*exp( f1))*py*_pz

  m.x = x*exp( f1) + px*f2*_pz
  m.y = y*exp(-f1) - py*f2*_pz

  m.px = px*exp(-f1)
  m.py = py*exp( f1)

  srotation(nil, m, -ang or 0)
end

function straight_fringe (elm, m, l) -- ADD DIRECTION
  local nmul, sdir in m
  if nmul == 0 then return end

  if pos == 'entry'  then
	  local e, h = elm.e1 or 0, elm.h1 or 0
    if sdir == 1 then
    	yrotation          (elm, m,  e)
    	sbend_face         (elm, m,  h)
    	sbend_fringe_field (elm, m,  'entry')
	  else
      sbend_fringe_field (elm, m,  'entry')
      sbend_face         (elm, m,  h)
      yrotation          (elm, m,  e)
    end
    mult_fringe_field    (elm, m,  'entry')
    straight_fringe_field(elm, m,  'entry')
    return
  end

  if pos == 'exit'  then
    local e, h = elm.e2 or 0, elm.h2 or 0
    straight_fringe_field(elm, m,  'exit')
    mult_fringe_field    (elm, m,  'exit')
  	if sdir == 1 then
      sbend_fringe_field (elm, m,  'exit')
      sbend_face         (elm, m,  h)
      yrotation          (elm, m,  e)
	  else
  	  yrotation          (elm, m,  e)
  	  sbend_face         (elm, m,  h)
  	  sbend_fringe_field (elm, m,  'exit')
    end
    return
  end
end

__help['track: solenoid fringe'] = [=[

]=]

function solenoid_fringe (elm, m, l)
--  local ksi in m
--  if abs(ksi) < minang then -- should never occur due to precedences
    return straight_fringe(elm,m,l)
--  end
  -- error("NYI") -- TODO
end

__help['track: rfcavity fringe'] = [=[
Ref??
]=]

local rfcavity_fringe_ptc -- forward

function rfcavity_fringe (elm, m, l)
--  local v in m
--  if abs(v) < minvolt then -- should never occur due to precedences
    return straight_fringe(elm,m,l)
--  end
  -- rfcavity_fringe_ptc(elm, m, l) -- TODO
end

__help['track: sbend fringe'] = [=[
Ref??
]=]

local sbend_fringe_ptc -- forward

function sbend_fringe (elm, m, l)
  local angle in elm
  if abs(angle) < minang then -- should never occur due to precedences
    return straight_fringe(elm,m,l)
  end
  sbend_fringe_ptc(elm, m, l)
end

-- sbend ----------------------------------------------------------------------o

__help['track: exact sbend'] = [[
  [Forest98] p.360, ch. 12.2.3
  [Forest02] p.76 (general explaination)
  [Forest06] p.5365, eq. 126, exact drift
  [MADX-PTC] Function SPROTR in file Sh_def_kind.f90
]]

local function sbend_face (elm, m, h)
  local x, px, y, py, t, pt, bdir in m
  local k0 = elm.k0 or 0
  local _beta = 1/m.beam.beta

  if bdir == 1 then -- to insure reversal symmetry ; horizontal wedge (only dir)
    px = px +(bdir*k0*h/2)*x^2
  end

  local ptt =        sqrt(1 + 2*pt*_beta + pt^2 - px^2)
  local xi  = bdir*sqrt(1 + 2*pt*_beta + pt^2)*k0*h/ptt^2
  local dxi_px   =  2*px*xi/ptt^2
  local dxi_ddel = -2*xi*(1 + pt)/ptt^2

  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(bdir == -1) then -- to insure reversal symmetry; horizontal wedge (only dir)
    px = px +(bdir*k0*h/2)*x^2
  end

  m.x  = x
  m.px = px
  m.py = py
  m.t  = t
end

local function sbend_wedge (elm, m, e)
  local x, px, y, py, t, pt, T, bdir in m
  local k0 = elm.k0 or 0
  local _beta = 1/m.beam.beta
  local b1 = bdir*k0
  local sa, ca, sa2 = sin(e), cos(e), sin(2*e)

  if b1 == 0 then
    yrot_track (elm, m, e) return
  end

  local pz  = sqrt(1 + 2*_beta*pt + pt^2 - px^2 - py^2)
  local ptt = sqrt(1 + 2*_beta*pt + pt^2        - py^2)
  local pxt = px*ca + (pz - b1*x)*sa
  local pzs = sqrt(1 + 2*_beta*pt + pt^2 - py^2 - pxt^2)
  local yt  = (e + asin(px/ptt) - asin(pxt/ptt)) / b1

  m.x  = x*ca + (x*px*sa2 + sa^2*(2*x*pz-b1*x^2)) / (pzs+pz*ca-px*sa)
  m.px = pxt
  m.y  = y + py * yt
  m.t  = t - yt*(_beta + pt)
end

local function sbend_fringe_field (elm, m, pos)
--  TODO: FINT and FINTX to have different FINT at entrance /exit
  local x, px, y, py, t, pt, T in m
  local k0   = elm.k0 or 0
  local _beta = 1/m.beam.beta
-- soft edge
  local fint = elm.fint or 0  -- fringe field integral at entrance and exit of the bend
  local hgap = elm.hgap or 0  -- half gap of the magnet
  local fsad = 0

  if fint*hgap ~= 0.0 then  -- soft edge
      fsad = 1/(fint*hgap*2)/36
  end

  local b = m.beam.charge*k0
  if pos == 'exit' then
     b = -b
  end

  local pz = sqrt(1 + 2*_beta*pt + pt^2 - px^2 - py^2)
  local time_fac = _beta + pt
  local rel_p    = sqrt(1 + 2*_beta*pt + pt^2)

  local c3=b^2*fsad/rel_p

  local xp, yp = px/pz or 0, py/pz or 0
  local d  = matrix {
    { (1+xp^2) / pz, xp*yp    / pz, -time_fac*xp /pz^2 },
    {  xp*yp   / pz, (1+yp^2) / pz, -time_fac*yp /pz^2 },
    { -xp          ,   -yp        ,  time_fac    /pz   },
  }

  local fi0 = atan((xp/(1 + yp^2)))-2*b*fint*hgap*(1 + xp^2*(2 + yp^2))*pz
  local co2 = b/cos(fi0)^2
  local co1 = co2/(1 + (xp/(1 + yp^2))^2)

  local fi = {0, 0, 0}
  fi[1] =    co1 /       (1 + yp^2) - 2*co2*b*fint*hgap* (2*xp*(2 + yp^2)*pz)
  fi[2] =-2*co1*xp*yp/(1 + yp^2)^2  - 2*co2*b*fint*hgap* (2*xp^2*yp)*pz
  fi[3] =                           - 2*co2*b*fint*hgap* (1 + xp^2*(2 + yp^2))

  fi0=b*tan(fi0)

  b=0
  for i=1,3 do -- y column
    b=fi[i]*d:get(i,2) + b
  end
  m.y  = 2*m.y/(1 + sqrt(1 - 2*b*m.y))
  m.py = py - fi0*m.y

  b=0
  for i=1,3 do -- x column
    b=fi[i]*d:get(i,1) + b
  end
  m.x = m.x+0.5*b*m.y^2

  b=0
  for i=1,3 do -- z column
    b=fi[i]*d:get(i,3) + b
  end
  m.t = t + 0.5*b*m.y^2

  -- soft edge
  m.py = m.py - 4*c3*m.y^3
  m.t  = m.t  +   c3*m.y^4/rel_p^2*time_fac
end

function sbend_fringe_ptc (elm, m, pos) -- has forward decl
  if elm.kill_fringe then return end

  if pos == 'entry' and not elm.kill_ent_fringe then
    local e, h = elm.e1 or 0, elm.h1 or 0
    if abs(e) < minang then
      sbend_face         (elm, m,  h)
      sbend_fringe_field (elm, m,  'entry')
    else
      yrotation          (elm, m,  e)
      sbend_face         (elm, m,  h)
      sbend_fringe_field (elm, m,  'entry')
      sbend_wedge        (elm, m, -e)
    end
    return
  end

  if pos == 'exit' and not elm.kill_exi_fringe then
    local e, h = elm.e2 or 0, elm.h2 or 0
    if abs(e) < minang then
      sbend_fringe_field (elm, m, 'exit')
      sbend_face         (elm, m,  h)
    else
      sbend_wedge        (elm, m, -e)
      sbend_fringe_field (elm, m, 'exit')
      sbend_face         (elm, m,  h)
      yrotation          (elm, m,  e)
    end
    return
  end
end

-- cavities -------------------------------------------------------------------o

--[=[
local function rfcavity_kick (elm, m, l) -- RF CAV, CRAB, RFMULT with l>0 kick (kickcavr)
  local x, px, y, py, t, pt, nmul, bdir in m
  local tt, phase, phase0 = elm.tt or 0, elm.phase or 0, elm.phase0 or 0
  local ff, ph = elm.f or {}, elm.ph or {}
  local a, r = elm.a or  0, elm.r or 0

  -- nocavity.and.(.not.always_on)  return
  if 'nocavity' and (not 'always_on') then return end -- non-sense!!!

  local pc = m.beam.pc
  local volt_c = 1.0e-3
  local freq, harm, volt, lag = elm.freq or 0, elm.harm or 0, elm.volt or 0, elm.lag or 0
  local nbessel, nf = elm.nbessel or  0, elm.nf or 1
  if freq == 0 and harm == 0 then
    error("either frequency or harmonic number has to be defined for RF Cavity "
          .. elm.name)
  end

  local omega = twopi*freq / clight
  local vl    = m.bdir*l*volt*volt_c/pc
  local df, f, r2 = 0, 1, 1

  for ko=1,nf do -- over modes
    for i=1,nbessel do
      r2  = -r2*(ko*omega)^2/(4*(i+1)^2)
      dr2 =  r2*i
      df  =  df + 2*dr2
      r2  =  r2*(x^2+y^2)
      f   =  f + r2
    end

    if nbessel > 0 then
      m.px = px-x*ff[ko]* df * vl * cos(arg)/(ko*omega)
      m.py = py-y*ff[ko]* df * vl * cos(arg)/(ko*omega)
    end

    local arg = ko*omega*(t + tt) + phas + ph[ko] + phase0
    m.pt = m.pt -ff[ko]*f*vl*sin(arg)

    local lw, byt = elm.l>0 and l/elm.l or 1
    local by = lw*knl[nmul] / fact(nmul-1)
    local bx = lw*ksl[nmul] / fact(nmul-1)

    for i=nmul-1,1,-1 do
      byt = x*by - y*bx + lw*knl[i] / fact(i-1)
      bx  = y*by + x*bx + lw*ksl[i] / fact(i-1)
      by  = byt
    end

    m.px = m.px - ff[ko]*bdir*by/pc*(a + r*cos(arg))
    m.py = m.py + ff[ko]*bdir*bx/pc*(a + r*cos(arg))

    by = -lw*knl[nmul] / fact(nmul-1) / nmul
    bx = -lw*ksl[nmul] / fact(nmul-1) / nmul

    for i=nmul-1,1,-1 do
      byt = x*by - y*bx - lw*knl[i] / fact(i-1) / i
      bx  = y*by + x*bx - lw*ksl[i] / fact(i-1) / i
      by  = byt
    end

    byt = x*by - y*bx
    bx  = y*by + x*bx
    by  = byt

    m.pt = m.pt + ff[ko]*(ko*omega)*bdir*by/pc * r *sin(arg)
  end --over modes
end

local function rfcavity_thin_kick (elm, m, l) -- for tracking thin rf mult
  local x, px, y, py, t, pt, nmul, bdir in m
  local tt, phase, phase0 = elm.tt or 0, elm.phase or 0, elm.phase0 or 0
  local ff, ph = elm.f or {}, elm.ph or {}
  local a, r = elm.a or  0, elm.r or 0
  -- nocavity.and.(.not.always_on)  return
  if 'nocavity' and (not 'always_on') then return end -- non-sense!!!

  local pc = beam.pc
  local volt =elm.volt
  local volt_c = 1.0e-3
  local omega = twopi*freq / clight

  for ko=1,nf do -- over modes
    local arg = ko*omega*(t + tt) + phas + ph[ko] + phase0
    m.pt = m.pt - bdir*volt_c*ff[ko]*volt*sin(arg)/pc

    -- BUG: knl, ksl never retrieved
    local lw, byt = elm.l>0 and l/elm.l or 1
    local by = lw*knl[nmul] / fact(nmul-1)
    local bx = lw*ksl[nmul] / fact(nmul-1)

    for i=nmul-1,1,-1 do
      byt = x*by - y*bx + lw*knl[i] / fact(i-1)
      bx  = y*by + x*bx + lw*ksl[i] / fact(i-1)
      by  = byt
    end

    m.px = m.px - ff[ko]*bdir*by/pc*(a + r*cos(arg))
    m.py = m.py + ff[ko]*bdir*bx/pc*(a + r*cos(arg))

    by = -lw*knl[nmul] / fact(nmul-1) / nmul
    bx = -lw*ksl[nmul] / fact(nmul-1) / nmul

    for i=nmul,2,-1 do
      byt = x*by - y*bx - lw*knl[i-i] / fact(i-2) / (i-1)
      bx  = y*by + x*bx - lw*ksl[i-i] / fact(i-2) / (i-1)
      by  = byt
    end

    byt = x*by - y*bx
    bx  = y*by + x*bx
    by  = byt

    m.t = t - ff[ko]*(ko*omega)*bdir*by/pc * r *sin(arg)
  end --over modes
end

local function adjust_time_cav(elm, m, pos)
  local x, px, y, py, t, pt, T, dir in m
  local _beta = 1/m.beam.beta

  local delta_e = elm.delta_e or 0
  if pos == 'entry' then
    delta_e = pt
    if 'nocavity' and (not 'always_on') then return end -- non-sense!!!

    if elm.l < minlen then
      rfcavity_thin_kick(elm, m)
      delta_e=(m.pt - delta_e)*m.beam.pc -- side effects?
    end
  else
    if elm.l < minlen then return end
    t = t - elm.l*_beta -- TO CHECK      -- side effects?
    delta_e = (pt - delta_e)*m.beam.pc   -- side effects?
  end
end

local function rfcavity_fringe_ptc (elm, m, pos)
  local x, px, y, py, t, pt, T, dir in m

  if ('nocavity' and (not 'always_on')) then return end -- non-sense!!!

  if elm.kill_fringe     == 1                     or
     elm.kill_ent_fringe == 1 and pos ==  'entry' or
     elm.kill_ext_fringe == 1 and pos ==  'exit' then return end

  adjust_time_cav(elm, m, 'entry')

  if not ('fringe' or elm.nbessel == -1) then return end
  if elm.l < minlen then return end

  local kbmad = 0
  if elm.nbessel == -1 then kbmad = 1 end

  local z = 0
  if pos == 'entry' then z = elm.l end
  local freq, volt = elm.freg or 0, elm.volt or 0
  local volt_c = 1.0e-3
  local pc = m.beam.pc
  local omega = twopi*freq / clight
  local vl = m.bdir*volt*volt_c/pc
  if pos == 'exit' then vl = -vl end

  local tt, phase, phase0 = elm.tt or 0, elm.phase or 0, elm.phase0 or 0
  local ph = elm.f or {}, elm.ph or {}
  local nf = elm.nf or 1

  for ko=1,nf do   -- over modes
   s1 = cos(kbmad*ko*omega*z)*sin(ko*omega*(t+tt) + phas + phase0 + ph[ko])
   c1 = cos(kbmad*ko*omega*z)*cos(ko*omega*(t+tt) + phas + phase0 + ph[ko])
   m.px = m.px + 0.5*vl*s1*m.x
   m.py = m.py + 0.5*vl*s1*m.y
   m.pt = m.pt - 0.25*(m.x^2 + m.y^2)*vl*c1*omega*ko
  end

  adjust_time_cav(elm, m, 'exit')
end
--]=]

-- integrator schemes ---------------------------------------------------------o

local integrators = require 'madl_sympint'

-- sub element tracking -------------------------------------------------------o

local function subelm_track (elm, m, integrator)
  local drift_track, kick_track in elm
  local selm, selm_s in m

  -- save current sub-element and sub-position
  local nselm = #selm
  selm[nselm], selm_s[nselm] = elm, 0

  -- sub-tracking
  local n = #elm
  for i=1,n do
    -- implicit drift
    local dl = elm:spos(i) - selm_s[nselm]
    if dl >= minlen then
      integrator(elm, m, dl, drift_track, kick_track)
      selm_s[nelm] = selm_s[nelm] + dl
    elseif dl <= -minlen then
      error(string.format("invalid sub-element %s['%s'] relative position at %sm",
            m.sequence.name, elm[i].name, dl))
    end

    -- sub-element
    elm[i]:track(m)
    selm_s[nelm] = selm_s[nelm] + elm[i].l
  end

  -- last implicit drift (if any)
  local dl = elm.l - elm:spos(n)
  if dl >= minlen then
    integrator(elm, m, dl, drift_track, kick_track)
  elseif dl <= -minlen then
    error(string.format("invalid sub-element %s['%s'] relative position at %sm",
          m.sequence.name, elm[n].name, dl))
  end

  -- remove current sub-element
  selm[nselm] = nil
end

-- sub element backtracking ---------------------------------------------------o

local function subelm_backtrack (elm, m, integrator)
  local drift_track, kick_track in elm
  local selm, selm_s in m

  -- save current sub-element and sub-position
  local nselm = #selm
  selm[nselm], selm_s[nselm] = elm, 0

  -- sub-tracking
  local n = #elm
  for i=n,1,-1 do
    -- implicit drift
    local dl = selm_s[nselm] - (elm:spos(i) + elm[i].l)
    if dl <= -minlen then
      integrator(elm, m, dl, drift_track, kick_track)
      selm_s[nelm] = selm_s[nelm] + dl
    elseif dl >= minlen then
      error(string.format("invalid sub-element %s['%s'] relative position at %sm",
            m.sequence.name, elm[i].name, dl))
    end

    -- sub-element
    elm[i]:backtrack(m)
    selm_s[nelm] = selm_s[nelm] - elm[i].l
  end

  -- last implicit drift (if any)
  local dl = -elm:spos(1)
  if dl <= -minlen then
    integrator(elm, m, dl, drift_track, kick_track)
  elseif dl <= -minlen then
    error(string.format("invalid sub-element %s['%s'] relative position at %sm",
          m.sequence.name, elm[n].name, dl))
  end

  -- remove current sub-element
  selm[nselm] = nil
end

-- element tracking -----------------------------------------------------------o

local function patch_track (elm, m, pos)
  if pos == 'entry'
  then entry_patch(elm, m)
  else  exit_patch(elm, m)
  end
end

local function thin_track (elm, m) -- e.g. multipole
  get_mult(elm, m) -- load strengths in map

  if m.ktot == 0 then return end

  patch_track(elm, m, 'entry')
  thin_kick  (elm, m)
  patch_track(elm, m, 'exit' )
end

local function drift_track (elm, m) -- e.g. monitor
  local no_selm = #m.selm == 0

  patch_track(elm, m, 'entry')

  if no_selm and #elm == 0 then -- any surrounding or sub-elements?
    integrators.thick(elm, m, m.el, straight_drift)
  elseif no_selm then
    subelm_track(elm, m, straight_drift)
  else
    subelm_track(elm, m, integrators[elm.method or m.mth])
  end

  patch_track(elm, m, 'exit')
end

local function thick_track (elm, m) -- e.g. sbend, rbend, quadrupole, ...
  get_kmult(elm, m)   -- load strengths in map

  if m.ktot == 0 then -- no strengths in elm
    return drift_track(elm, m)
  end

  patch_track(elm, m, 'entry')
  elm:fringe_track(m, 'entry')

  local integrator = integrators[elm.method or m.mth]
  if #elm == 0 then   -- any sub-elements?
    integrator  (elm, m, m.el, elm.drift_track, elm.kick_track)
  else
    subelm_track(elm, m, integrator)
  end

  elm:fringe_track(m, 'exit')
  patch_track(elm, m, 'exit')
end

-- element back tracking ------------------------------------------------------o

local patch_backtrack = patch_track

local function thin_backtrack (elm, m) -- e.g. multipole
  get_mult(elm, m) -- load strengths in map

  if m.ktot == 0 then return end

  patch_backtrack(elm, m, 'exit' )
  thin_kick      (elm, m)
  patch_backtrack(elm, m, 'entry')
end

local function drift_backtrack (elm, m) -- e.g. monitor
  local no_selm = #m.selm == 0

  patch_backtrack(elm, m, 'exit')

  if no_selm and #elm == 0 then -- any surrounding or sub-elements?
    integrators.thick(elm, m, m.el, straight_drift)
  elseif no_selm then
    subelm_backtrack(elm, m, straight_drift)
  else
    subelm_backtrack(elm, m, integrators[elm.method or m.mth])
  end

  patch_backtrack (elm, m, 'entry')
end


local function thick_backtrack (elm, m) -- e.g. sbend, rbend, quadrupole, ...
  get_kmult(elm, m)   -- load strengths in map

  if m.ktot == 0 then -- no strengths in elm
    return drift_backtrack(elm, m)
  end

  patch_backtrack(elm, m, 'exit')
  elm:fringe_track    (m, 'exit')

  local integrator = integrators[elm.method or m.mth]
  if #elm == 0 then   -- any sub-elements ?
    integrator      (elm, m, m.el, elm.drift_track, elm.kick_track)
  else
    subelm_backtrack(elm, m, integrator)
  end

  elm:fringe_track    (m, 'entry')
  patch_backtrack(elm, m, 'entry')
end

-- sequence tracking ----------------------------------------------------------o

local dft = element.drift {}

local function sequ_track (seq, map, rng_, ntrn_)
  local tbl = map.mtable

  -- retrieve information
  local save, drift, nlen in map

  -- dynamic tracking
  for i,elm in seq:iter(rng_, ntrn_, 1) do

    -- implicit drift
    local dl in map
    if dl >= minlen then
      straight_drift(nil, map, dl)
      if drift and elm:is_selected(observe) then
        map.ndft = map.ndft+1
        dft.name = 'DRIFT_'..map.ndft
        tbl:fill(dft, map, dl)
      end
      map.s, map.dl = map.s + dl, 0 -- update s
    elseif dl <= -minlen then
      error(string.format("invalid implicit drift of %sm in %s['%s'] at %sm",
            dl, seq.name, elm.name, map.s))
    end

    -- sequence element
    map.el = elm.l
    local stop = elm:track(map)
    if stop then return tbl, map end -- stop (e.g. slink)

    local el in map
    if save and elm:is_selected(observe) then
      tbl:fill(elm, map, el)
    end
    map.s = map.s + el -- update s

    -- update dl
    if i == nlen
    then map.dl = 0
    else map.dl = seq:spos(i+1) - (seq:spos(i) + el) -- patch and dl?
    end
  end

  return tbl, map
end

-- sequence backtracking ------------------------------------------------------o

local function sequ_backtrack (seq, map, rng_, ntrn_)
  local tbl = map.mtable

  -- retrieve information
  local save, drift, nlen in map

  -- dynamic tracking
  for i,elm in seq:iter(rng_, ntrn_, -1) do

    -- implicit drift
    local dl in map
    if dl <= -minlen then
      straight_drift(nil, map, dl)
      if drift and elm:is_selected(observe) then
        map.ndft = map.ndft+1
        dft.name = 'DRIFT_'..map.ndft
        tbl:fill(dft, map, dl)
      end
      map.s, map.dl = map.s + dl, 0 -- update s
    elseif dl >= minlen then
      error(string.format("invalid implicit drift of %sm in %s['%s'] at %sm",
            dl, seq.name, elm.name, map.s))
    end

    -- sequence element
    map.el = -elm.l
    local stop = elm:backtrack(map)
    if stop then return tbl, map end -- stop (e.g. slink)

    local el in map
    if save and elm:is_selected(observe) then
      tbl:fill(elm, map, el)
    end
    map.s = map.s + el -- update s

    -- update dl
    if i == 1
    then map.dl = 0
    else map.dl = seq:spos(i-1) + seq[i-1].l - seq:spos(i) -- patch and dl?
    end
  end

  return tbl, map
end

-- load maps into elements ----------------------------------------------------o

local E = element

E.patch_element :set_functions {    track = patch_track,
                                backtrack = patch_backtrack }

E. thin_element :set_functions {    track = thin_track,
                                backtrack = thin_backtrack }

E.drift_element :set_functions {    track = drift_track,
                                backtrack = drift_backtrack }

E.thick_element :set_functions {    track = thick_track,
                                backtrack = thick_backtrack }

E.sequence :set_readonly(false)
                :set_functions {    track = sequ_track,
                                backtrack = sequ_backtrack } :set_readonly()

-- specialized maps

E.thick_element :set_functions { -- default specializations
                              drift_track = straight_drift,
                               kick_track = straight_kick,
                             fringe_track = straight_fringe }

E.sbend :set_functions {      drift_track = sbend_drift,
                               kick_track = polar_kick,
                             fringe_track = sbend_fringe }

E.combined :set_functions   { drift_track = combined_drift,
                               kick_track = combined_kick,
                             fringe_track = sbend_fringe }

E.rbend :set_functions      { drift_track = rbend_drift,
                               kick_track = rfcavity_kick,
                             fringe_track = straight_fringe }

E.quadrupole :set_functions { drift_track = quadrupole_drift,
                               kick_track = quadrupole_kick,
                             fringe_track = straight_fringe }

E.solenoid :set_functions   { drift_track = solenoid_drift,
                               kick_track = rfcavity_kick,
                             fringe_track = solenoid_fringe }

E.rfcavity :set_functions   { drift_track = solenoid_drift,
                               kick_track = rfcavity_kick,
                             fringe_track = rfcavity_fringe }

E.twcavity :set_functions   { drift_track = twcavity_drift,
                               kick_track = twcavity_kick,
                             fringe_track = twcavity_fringe }

E.changedir   :set_functions { track = \_,m -> changedir(nil, m) }
E.srotation   :set_functions { track = \e,m -> srotation(nil, m, e.angle) }
E.yrotation   :set_functions { track = \e,m -> yrotation(nil, m, e.angle) }
E.xrotation   :set_functions { track = \e,m -> xrotation(nil, m, e.angle) }
E.translation :set_functions { track = \e,m -> translate(nil, m, e.dx or 0,
                                                      e.dy or 0, e.dz or 0) }

-- track helper ---------------------------------------------------------------o

local _trk = {}

local function fill_mtable (self, elm, m, l)
  -- keep order!
  self = self+{ elm.name, elm.kind, m.s+l, l, m.x, m.px, m.y, m.py, m.t, m.pt }
end

local function make_mtable (self)
  if self.mtable then
    assert(self.mtable[_trk] == _trk, "corrupted tracked mtable")
    return self.mtable
  end

  local sequence, drift, range, nturn, direction in self
  local nrow = sequence:length_of(range, nturn) * (drift and 2 or 1)

  return mtable 'track' {
    type='track', title=sequence.name, direction=direction,
    {'name'}, 'kind', 's', 'l', 'x', 'px', 'y', 'py', 't', 'pt',
    [_trk]=_trk,
  } : reserve(nrow) : set_functions { fill=fill_mtable }
end

local function update_mflow (self)
  local sequence, beam, direction, mflow in self

  -- sanity checks
  assert(mflow[_trk] == _trk, "corrupted tracked mflow" )
  if mflow.save then
    assert(mflow.mtable[_trk] == _trk, "corrupted tracked mtable")
  end
  assert(is_beam    (mflow.beam    ), "corrupted tracked beam"    )
  assert(is_sequence(mflow.sequence), "corrupted tracked sequence")

  -- update with proper precedence
  if sequence and sequence ~= mflow.sequence then -- new sequence -> update
    assert(is_sequence(sequence), "invalid sequence")
    local len, sdir = #sequence, sequence.direction
    mflow.dl, mflow.sdir, mflow.nlen = 0, sdir, len
    mflow.sequence, mflow.bdir = sequence, sdir*mflow.beam.charge
  end
  if direction and direction ~= mflow.mdir then -- new direction -> update
    mflow.dl, mflow.mdir = 0, direction
  end
  if beam and beam ~= mflow.beam then -- new beam -> update
    assert(is_beam(beam), "invalid beam")
    mflow.beam, mflow.bdir = beam, mflow.sdir*mflow.mdir*beam.charge
  end

  return mflow
end

local function make_mflow (self)
  if self.mflow then
    return update_mflow(self)
  end

  local beam, sequence in self -- mandatory
  assert(is_beam(beam)        , "invalid beam")
  assert(is_sequence(sequence), "invalid sequence")

  local s, x, px, y, py, t, pt, X0 in self
  s  = s  or X0.s  or          0
  x  = x  or X0.x  or X0[1] or 0
  px = px or X0.px or X0[2] or 0
  y  = y  or X0.y  or X0[3] or 0
  py = py or X0.py or X0[4] or 0
  t  = t  or X0.t  or X0[5] or 0
  pt = pt or X0.pt or X0[6] or 0

  local save, drift, nstep, nslice, naction, method, total_path in self
  local nlen, sdir, mdir = #sequence, sequence.direction, self.direction
  local T = total_path == true and 1 or 0

  local mtable = save and make_mtable(self)

  return { x=x, px=px, y=y, py=py, t=t, pt=pt, T=T,          -- initial vars
           nmul=0, knl={}, ksl={}, pnl={}, psl={},           -- internal vars
           s=s, dl=0, el=0, nlen=nlen, ndft=0,               -- internal states
           selm={}, selm_s={},                               -- subelement stack
           nst=nstep, nsl=nslice, nact=naction, mth=method,  -- user setup
           save=save, drift=save and drift,                  -- user status
           sdir=sdir, mdir=mdir, bdir=sdir*mdir*beam.charge, -- directions
           sequence=sequence, beam=beam, mtable=mtable,      -- attached objects
           [_trk]=_trk }                                     -- tag
end

-- track command --------------------------------------------------------------o

local function exec (self)
  assert(self[_trk], "invalid argument #1 (track expected)")
  local map = make_mflow(self)                          -- build mflow
  local range, nturn in self                            -- range of tracking

  if map.mdir == 1
  then return map.sequence:    track(map, range, nturn) --     track sequence
  else return map.sequence:backtrack(map, range, nturn) -- backtrack sequence
  end
end

-- track command template

local track = Command 'track' {
  nturn=0, direction=1, total_path=false, drift=true, save=true,
  nstep=1, nslice=2, method='simple',
  exec=exec, [_trk]=true,
} :set_functions { naction=empty_track }
  :set_readonly()

-- end ------------------------------------------------------------------------o
return {
  track  = track,
  __help = __help,
}