Skip to content

Combine curvature and solenoid multipole maps #285

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Mar 3, 2023
Merged

Combine curvature and solenoid multipole maps #285

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8842936
Combine curvature and solenoid multipole maps
jgray-19 Mar 3, 2023
31824c0
Simplify solenoid for multipole
Mar 3, 2023
958d0b0
cleanup
Mar 3, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 38 additions & 49 deletions src/madl_dynmap.mad
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
36 changes: 19 additions & 17 deletions src/madl_etrck.mad
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down