Thin lense rf cavity

src/MatrixKick/MatrixKick.jl

@@ -15,6 +15,13 @@ Base.@kwdef struct Quadrupole{T}
   Bn1::T  # quadrupole gradient / (T·m^-1)
 end
 
+Base.@kwdef struct ThinLensRFCavity{T}
+  L::T    # RF-cavity length     / m
+  V::T    # Cavity voltage       / V
+  k::T    # RF wavenumber (2π/λ) / raduans·m^-1
+  phi0::T # RF phase offset      / radians
+end
+
 
 function track!(bunch::Bunch, ele::MatrixKick.Drift; work=get_work(bunch, Val{1}()))
 #=
@@ -67,7 +74,7 @@ function track!(bunch::Bunch, ele::MatrixKick.Quadrupole; work=get_work(bunch, V
   v_work = StructArray{Coord{eltype(work[1])}}((work[1], work[2], work[3], work[4], work[5], work[6]))
 
   trackQuadMx!(v_work, v, k2_num, L / 2)
-  trackQuadK!( v, v_work, bunch.beta_gamma_ref, L)
+  trackQuadK!( v, v_work, bunch.beta_gamma_ref, L) 
   trackQuadMx!(v_work, v, k2_num, L / 2)
 
   v .= v_work
@@ -150,5 +157,55 @@ function trackQuadK!(vf, vi, betgam_ref, s)
 end # function trackQ!M::Quadrupole()
 
 
+
+"""
+function track!(bunch::Bunch, ele::MatrixKick.ThinLensRFCavity; work=get_work(bunch, Val{1}()))
+
+Apply a thin  RF cavity map in-place, where momenta (p_x, p_y, p_z)
+and rest mass are expressed in eV (treating c=1 for momentum units),
+while the longitudinal coordinate z is in meters.
+
+**Algorithm**:
+1) Half-drift by L/2 using the current momentum (which determines velocity).
+2) Kick in pz by `Δpz = q*V (eV) × `cos(k*z + phi0)`
+3) Half-drift by L/2 using the updated momentum.
+
+"""
+function track!(bunch::Bunch, ele::MatrixKick.ThinLensRFCavity; work=get_work(bunch, Val{1}()))
+  v = bunch.v
+  q = chargeof(bunch.species)
+  m = massof(bunch.species)
+  k = ele.k
+  # pleccing minus according to 4.46 in manual
+  dE = -q * ele.V # energy gain/loss in eV 
+  L = ele.L
+
+  @FastGTPSA! begin
+     # total momentum p [eV]
+  #------------------------------------------------------
+  # Step 1: Half-drift using old momentum
+  #------------------------------------------------------
+  # @. work[1] = sqrt(v.px^2 + v.py^2 + v.pz^2 + m^2)
+  @. work[1] = sqrt((1.0 + v.pz)^2 - (v.px^2 + v.py^2))  # P_s
+  @. work[1] = ifelse(work[1] > 0.0, work[1], 1.0)
+  @. v.z += 0.5 * L * v.pz /  work[1] 
+
+  #------------------------------------------------------
+  # Step 2: RF momentum kick in pz
+  #------------------------------------------------------
+  @. v.pz += dE * cos(k*v.z + ele.phi0)
+
+  #------------------------------------------------------
+  # Step 3: Half-drift using updated momentum
+  #------------------------------------------------------
+  @. work[1] = sqrt((1.0 + v.pz)^2 - (v.px^2 + v.py^2))  # P_s
+  @. work[1] = ifelse(work[1] > 0.0, work[1], 1.0)
+  @. v.z += 0.5 * L * v.pz /  work[1] 
+
+  end
+  # Spin unchanged
+  return bunch
+end
+
 end # module MatrixKick
 

test/element_tests.jl

@@ -28,6 +28,16 @@ qf3 = MatrixKick.Quadrupole(L = lq3, Bn1 = +gr3)
 qd3 = MatrixKick.Quadrupole(L = lq3, Bn1 = -gr3)
 qf4 = MatrixKick.Quadrupole(L = lq4, Bn1 = +gr4)
 qd4 = MatrixKick.Quadrupole(L = lq4, Bn1 = -gr4)
+phi0_1 = 0.0
+phi0_2 = 1.2
+rfk = 2π / 0.1
+V1 = 1.e5
+V2 = 2.e5
+lrf1 = 0.0 
+lrf2 = 0.02 # m
+rf1 = MatrixKick.ThinLensRFCavity(L = lrf1, V = V1, k = rfk, phi0 = phi0_1) # 100 kV, zero length
+rf2 = MatrixKick.ThinLensRFCavity(L = lrf2, V = V2, k = rfk, phi0 = phi0_2) # 200 kV, non-zero length
+
 
 # define beams
 # -- species
@@ -155,6 +165,9 @@ pyf_qd4 = [ 0.,  0.,                      0.,                    -5.012687615023
 zf_qd4 =  [ 0.,  1.675151081511833e-8,   -1.680184012110547e-8,   1.6726365460219405e-8,  -3.78176641902771e-7, -3.78176641902771e-7 ]
 pzf_qd4 = [ 0.,  1.e-3,                  -1.e-3,                  1.e-3,                   1.e-3,                  1.e-3                 ]
 
+# beam6.qd4.final
+pze_rf2 = [72471.55089533473, 72471.5516661298, 72471.55012453966, 72471.5516661298,  72471.5516661298, 72471.5516661298]
+
 # test individual elements
 @testset "element_tests" begin
   # === drifts ===
@@ -285,6 +298,27 @@ pzf_qd4 = [ 0.,  1.e-3,                  -1.e-3,                  1.e-3,
   @test beam4.v.z  ≈  zf_qd4  (rtol=5.e-13)
   @test beam4.v.pz == pzf_qd4
 
+  beam5 = Bunch(species = e_minus, beta_gamma_ref = bg3,
+               x = copy(xi), px = copy(pxi), y = copy(yi), py = copy(pyi), z = copy(zi), pz = copy(pzi))
+  track!(beam5, rf1);
+  @test beam5.v.x  ==  xi  
+  @test beam5.v.px ==  pxi 
+  @test beam5.v.y  ==  yi  
+  @test beam5.v.py ==  pyi 
+  @test beam5.v.z  ==  zi  
+  expected_dpz = -chargeof(e_minus) * V1* cos.(rfk * zi .+ phi0_1)
+  @test isapprox(beam5.v.pz, pzi .+ expected_dpz; atol=1e-13)
+
+  beam6 = Bunch(species = e_minus, beta_gamma_ref = bg3,
+               x = copy(xi), px = copy(pxi), y = copy(yi), py = copy(pyi), z = copy(zi), pz = copy(pzi))
+  track!(beam6, rf2);
+  @test beam6.v.x  ==  xi  
+  @test beam6.v.px ==  pxi 
+  @test beam6.v.y  ==  yi  
+  @test beam6.v.py ==  pyi 
+  @test beam6.v.z  !=  zi  
+  @test isapprox(beam6.v.pz, pze_rf2; atol=1e-13)
+
 
 end