diff --git a/benchmark/Manifest.toml b/benchmark/Manifest.toml
new file mode 100644
index 00000000..97bcd14e
--- /dev/null
+++ b/benchmark/Manifest.toml
@@ -0,0 +1,557 @@
+# This file is machine-generated - editing it directly is not advised
+
+julia_version = "1.11.7"
+manifest_format = "2.0"
+project_hash = "de75b5803ad1cba7c19770d761c418a0aaedf835"
+
+[[deps.Accessors]]
+deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "MacroTools"]
+git-tree-sha1 = "856ecd7cebb68e5fc87abecd2326ad59f0f911f3"
+uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
+version = "0.1.43"
+
+    [deps.Accessors.extensions]
+    AxisKeysExt = "AxisKeys"
+    IntervalSetsExt = "IntervalSets"
+    LinearAlgebraExt = "LinearAlgebra"
+    StaticArraysExt = "StaticArrays"
+    StructArraysExt = "StructArrays"
+    TestExt = "Test"
+    UnitfulExt = "Unitful"
+
+    [deps.Accessors.weakdeps]
+    AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"
+    IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+    LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+    StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
+    Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+    Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+
+[[deps.Adapt]]
+deps = ["LinearAlgebra", "Requires"]
+git-tree-sha1 = "7e35fca2bdfba44d797c53dfe63a51fabf39bfc0"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "4.4.0"
+
+    [deps.Adapt.extensions]
+    AdaptSparseArraysExt = "SparseArrays"
+    AdaptStaticArraysExt = "StaticArrays"
+
+    [deps.Adapt.weakdeps]
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[[deps.ArrayInterface]]
+deps = ["Adapt", "LinearAlgebra"]
+git-tree-sha1 = "d81ae5489e13bc03567d4fbbb06c546a5e53c857"
+uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+version = "7.22.0"
+
+    [deps.ArrayInterface.extensions]
+    ArrayInterfaceBandedMatricesExt = "BandedMatrices"
+    ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
+    ArrayInterfaceCUDAExt = "CUDA"
+    ArrayInterfaceCUDSSExt = ["CUDSS", "CUDA"]
+    ArrayInterfaceChainRulesCoreExt = "ChainRulesCore"
+    ArrayInterfaceChainRulesExt = "ChainRules"
+    ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    ArrayInterfaceMetalExt = "Metal"
+    ArrayInterfaceReverseDiffExt = "ReverseDiff"
+    ArrayInterfaceSparseArraysExt = "SparseArrays"
+    ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
+    ArrayInterfaceTrackerExt = "Tracker"
+
+    [deps.ArrayInterface.weakdeps]
+    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
+    ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+    Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
+    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
+
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+version = "1.11.0"
+
+[[deps.AtomicAndPhysicalConstants]]
+git-tree-sha1 = "1b65be9932c3c65128919fb2286084701247ccfc"
+uuid = "5c0d271c-5419-4163-b387-496237733d8b"
+version = "0.8.2"
+
+[[deps.Atomix]]
+deps = ["UnsafeAtomics"]
+git-tree-sha1 = "29bb0eb6f578a587a49da16564705968667f5fa8"
+uuid = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
+version = "1.1.2"
+
+    [deps.Atomix.extensions]
+    AtomixCUDAExt = "CUDA"
+    AtomixMetalExt = "Metal"
+    AtomixOpenCLExt = "OpenCL"
+    AtomixoneAPIExt = "oneAPI"
+
+    [deps.Atomix.weakdeps]
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
+    OpenCL = "08131aa3-fb12-5dee-8b74-c09406e224a2"
+    oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b"
+
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+version = "1.11.0"
+
+[[deps.BeamTracking]]
+deps = ["Accessors", "Adapt", "AtomicAndPhysicalConstants", "GTPSA", "HostCPUFeatures", "KernelAbstractions", "MacroTools", "Random", "ReferenceFrameRotations", "SIMD", "SIMDMathFunctions", "SpecialFunctions", "StaticArrays", "Unrolled"]
+path = ".."
+uuid = "8ef5c10a-4ca3-437f-8af5-b84d8af36df0"
+version = "0.4.2"
+
+    [deps.BeamTracking.extensions]
+    BeamTrackingBeamlinesExt = "Beamlines"
+
+    [deps.BeamTracking.weakdeps]
+    Beamlines = "5bb90b03-0719-46b8-8ce4-1ef3afd3cd4b"
+
+[[deps.BenchmarkTools]]
+deps = ["Compat", "JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
+git-tree-sha1 = "7fecfb1123b8d0232218e2da0c213004ff15358d"
+uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+version = "1.6.3"
+
+[[deps.BitTwiddlingConvenienceFunctions]]
+deps = ["Static"]
+git-tree-sha1 = "f21cfd4950cb9f0587d5067e69405ad2acd27b87"
+uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b"
+version = "0.1.6"
+
+[[deps.CPUSummary]]
+deps = ["CpuId", "IfElse", "PrecompileTools", "Preferences", "Static"]
+git-tree-sha1 = "f3a21d7fc84ba618a779d1ed2fcca2e682865bab"
+uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
+version = "0.2.7"
+
+[[deps.CommonWorldInvalidations]]
+git-tree-sha1 = "ae52d1c52048455e85a387fbee9be553ec2b68d0"
+uuid = "f70d9fcc-98c5-4d4a-abd7-e4cdeebd8ca8"
+version = "1.0.0"
+
+[[deps.Compat]]
+deps = ["TOML", "UUIDs"]
+git-tree-sha1 = "9d8a54ce4b17aa5bdce0ea5c34bc5e7c340d16ad"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "4.18.1"
+weakdeps = ["Dates", "LinearAlgebra"]
+
+    [deps.Compat.extensions]
+    CompatLinearAlgebraExt = "LinearAlgebra"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "1.1.1+0"
+
+[[deps.CompositionsBase]]
+git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
+uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
+version = "0.1.2"
+weakdeps = ["InverseFunctions"]
+
+    [deps.CompositionsBase.extensions]
+    CompositionsBaseInverseFunctionsExt = "InverseFunctions"
+
+[[deps.ConstructionBase]]
+git-tree-sha1 = "b4b092499347b18a015186eae3042f72267106cb"
+uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+version = "1.6.0"
+
+    [deps.ConstructionBase.extensions]
+    ConstructionBaseIntervalSetsExt = "IntervalSets"
+    ConstructionBaseLinearAlgebraExt = "LinearAlgebra"
+    ConstructionBaseStaticArraysExt = "StaticArrays"
+
+    [deps.ConstructionBase.weakdeps]
+    IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+    LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[[deps.CpuId]]
+deps = ["Markdown"]
+git-tree-sha1 = "fcbb72b032692610bfbdb15018ac16a36cf2e406"
+uuid = "adafc99b-e345-5852-983c-f28acb93d879"
+version = "0.3.1"
+
+[[deps.Crayons]]
+git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
+uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
+version = "4.1.1"
+
+[[deps.Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+version = "1.11.0"
+
+[[deps.DocStringExtensions]]
+git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.9.5"
+
+[[deps.GTPSA]]
+deps = ["GTPSA_jll", "LinearAlgebra", "MacroTools", "Printf", "SpecialFunctions"]
+git-tree-sha1 = "d9312bc6e02a1fc219af56e78d2be7050519a04d"
+uuid = "b27dd330-f138-47c5-815b-40db9dd9b6e8"
+version = "1.5.1"
+
+[[deps.GTPSA_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenBLAS32_jll"]
+git-tree-sha1 = "421653d4ac27de5b1e34dd674169e8ac29c4adf5"
+uuid = "a4739e29-4b97-5c0b-bbcf-46f08034c990"
+version = "1.6.1+0"
+
+[[deps.HostCPUFeatures]]
+deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Preferences", "Static"]
+git-tree-sha1 = "af9ab7d1f70739a47f03be78771ebda38c3c71bf"
+uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0"
+version = "0.1.18"
+
+[[deps.IfElse]]
+git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
+uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
+version = "0.1.1"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+version = "1.11.0"
+
+[[deps.InverseFunctions]]
+git-tree-sha1 = "a779299d77cd080bf77b97535acecd73e1c5e5cb"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.17"
+
+    [deps.InverseFunctions.extensions]
+    InverseFunctionsDatesExt = "Dates"
+    InverseFunctionsTestExt = "Test"
+
+    [deps.InverseFunctions.weakdeps]
+    Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+    Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[deps.IrrationalConstants]]
+git-tree-sha1 = "b2d91fe939cae05960e760110b328288867b5758"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.2.6"
+
+[[deps.JLLWrappers]]
+deps = ["Artifacts", "Preferences"]
+git-tree-sha1 = "0533e564aae234aff59ab625543145446d8b6ec2"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.7.1"
+
+[[deps.JSON]]
+deps = ["Dates", "Logging", "Parsers", "PrecompileTools", "StructUtils", "UUIDs", "Unicode"]
+git-tree-sha1 = "b3ad4a0255688dcb895a52fafbaae3023b588a90"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "1.4.0"
+
+    [deps.JSON.extensions]
+    JSONArrowExt = ["ArrowTypes"]
+
+    [deps.JSON.weakdeps]
+    ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd"
+
+[[deps.KernelAbstractions]]
+deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs"]
+git-tree-sha1 = "b5a371fcd1d989d844a4354127365611ae1e305f"
+uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
+version = "0.9.39"
+
+    [deps.KernelAbstractions.extensions]
+    EnzymeExt = "EnzymeCore"
+    LinearAlgebraExt = "LinearAlgebra"
+    SparseArraysExt = "SparseArrays"
+
+    [deps.KernelAbstractions.weakdeps]
+    EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
+    LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.LayoutPointers]]
+deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"]
+git-tree-sha1 = "a9eaadb366f5493a5654e843864c13d8b107548c"
+uuid = "10f19ff3-798f-405d-979b-55457f8fc047"
+version = "0.1.17"
+
+[[deps.Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+version = "1.11.0"
+
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+version = "1.11.0"
+
+[[deps.LogExpFunctions]]
+deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "13ca9e2586b89836fd20cccf56e57e2b9ae7f38f"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.29"
+
+    [deps.LogExpFunctions.extensions]
+    LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
+    LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables"
+    LogExpFunctionsInverseFunctionsExt = "InverseFunctions"
+
+    [deps.LogExpFunctions.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+    InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
+
+[[deps.Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+version = "1.11.0"
+
+[[deps.MacroTools]]
+git-tree-sha1 = "1e0228a030642014fe5cfe68c2c0a818f9e3f522"
+uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+version = "0.5.16"
+
+[[deps.ManualMemory]]
+git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd"
+uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667"
+version = "0.1.8"
+
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+version = "1.11.0"
+
+[[deps.OpenBLAS32_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "46cce8b42186882811da4ce1f4c7208b02deb716"
+uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2"
+version = "0.3.30+0"
+
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.27+1"
+
+[[deps.OpenLibm_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+version = "0.8.5+0"
+
+[[deps.OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "1346c9208249809840c91b26703912dff463d335"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.6+0"
+
+[[deps.Parsers]]
+deps = ["Dates", "PrecompileTools", "UUIDs"]
+git-tree-sha1 = "7d2f8f21da5db6a806faf7b9b292296da42b2810"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "2.8.3"
+
+[[deps.PrecompileTools]]
+deps = ["Preferences"]
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
+uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+version = "1.2.1"
+
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "522f093a29b31a93e34eaea17ba055d850edea28"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.5.1"
+
+[[deps.Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+version = "1.11.0"
+
+[[deps.Profile]]
+uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
+version = "1.11.0"
+
+[[deps.Random]]
+deps = ["SHA"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+version = "1.11.0"
+
+[[deps.ReferenceFrameRotations]]
+deps = ["Crayons", "LinearAlgebra", "Printf", "Random", "StaticArrays"]
+git-tree-sha1 = "bbeaa5c8b63c066a134f2cad8cca1b8527cce9b5"
+uuid = "74f56ac7-18b3-5285-802d-d4bd4f104033"
+version = "3.1.2"
+
+    [deps.ReferenceFrameRotations.extensions]
+    ReferenceFrameRotationsZygoteExt = ["ForwardDiff", "Zygote"]
+
+    [deps.ReferenceFrameRotations.weakdeps]
+    ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+    Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+
+[[deps.Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "62389eeff14780bfe55195b7204c0d8738436d64"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.3.1"
+
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
+
+[[deps.SIMD]]
+deps = ["PrecompileTools"]
+git-tree-sha1 = "e24dc23107d426a096d3eae6c165b921e74c18e4"
+uuid = "fdea26ae-647d-5447-a871-4b548cad5224"
+version = "3.7.2"
+
+[[deps.SIMDMathFunctions]]
+deps = ["SIMD", "SLEEFPirates", "VectorizationBase"]
+git-tree-sha1 = "437d24266a7e029fa456036be76654051559fcc2"
+uuid = "d22a7203-ad50-4fbc-abc4-d6ac724cca58"
+version = "0.1.3"
+
+[[deps.SIMDTypes]]
+git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c"
+uuid = "94e857df-77ce-4151-89e5-788b33177be4"
+version = "0.1.0"
+
+[[deps.SLEEFPirates]]
+deps = ["IfElse", "Static", "VectorizationBase"]
+git-tree-sha1 = "456f610ca2fbd1c14f5fcf31c6bfadc55e7d66e0"
+uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
+version = "0.6.43"
+
+[[deps.SciMLPublic]]
+git-tree-sha1 = "0ba076dbdce87ba230fff48ca9bca62e1f345c9b"
+uuid = "431bcebd-1456-4ced-9d72-93c2757fff0b"
+version = "1.0.1"
+
+[[deps.SpecialFunctions]]
+deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+git-tree-sha1 = "f2685b435df2613e25fc10ad8c26dddb8640f547"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "2.6.1"
+
+    [deps.SpecialFunctions.extensions]
+    SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
+
+    [deps.SpecialFunctions.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+
+[[deps.Static]]
+deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools", "SciMLPublic"]
+git-tree-sha1 = "49440414711eddc7227724ae6e570c7d5559a086"
+uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
+version = "1.3.1"
+
+[[deps.StaticArrayInterface]]
+deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Static"]
+git-tree-sha1 = "96381d50f1ce85f2663584c8e886a6ca97e60554"
+uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718"
+version = "1.8.0"
+
+    [deps.StaticArrayInterface.extensions]
+    StaticArrayInterfaceOffsetArraysExt = "OffsetArrays"
+    StaticArrayInterfaceStaticArraysExt = "StaticArrays"
+
+    [deps.StaticArrayInterface.weakdeps]
+    OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[[deps.StaticArrays]]
+deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
+git-tree-sha1 = "eee1b9ad8b29ef0d936e3ec9838c7ec089620308"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.9.16"
+
+    [deps.StaticArrays.extensions]
+    StaticArraysChainRulesCoreExt = "ChainRulesCore"
+    StaticArraysStatisticsExt = "Statistics"
+
+    [deps.StaticArrays.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[deps.StaticArraysCore]]
+git-tree-sha1 = "6ab403037779dae8c514bad259f32a447262455a"
+uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+version = "1.4.4"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "ae3bb1eb3bba077cd276bc5cfc337cc65c3075c0"
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+version = "1.11.1"
+
+    [deps.Statistics.extensions]
+    SparseArraysExt = ["SparseArrays"]
+
+    [deps.Statistics.weakdeps]
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.StructUtils]]
+deps = ["Dates", "UUIDs"]
+git-tree-sha1 = "9297459be9e338e546f5c4bedb59b3b5674da7f1"
+uuid = "ec057cc2-7a8d-4b58-b3b3-92acb9f63b42"
+version = "2.6.2"
+
+    [deps.StructUtils.extensions]
+    StructUtilsMeasurementsExt = ["Measurements"]
+    StructUtilsTablesExt = ["Tables"]
+
+    [deps.StructUtils.weakdeps]
+    Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
+    Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+
+[[deps.TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.3"
+
+[[deps.UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+version = "1.11.0"
+
+[[deps.Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+version = "1.11.0"
+
+[[deps.Unrolled]]
+deps = ["MacroTools"]
+git-tree-sha1 = "6cc9d682755680e0f0be87c56392b7651efc2c7b"
+uuid = "9602ed7d-8fef-5bc8-8597-8f21381861e8"
+version = "0.1.5"
+
+[[deps.UnsafeAtomics]]
+git-tree-sha1 = "b13c4edda90890e5b04ba24e20a310fbe6f249ff"
+uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f"
+version = "0.3.0"
+
+    [deps.UnsafeAtomics.extensions]
+    UnsafeAtomicsLLVM = ["LLVM"]
+
+    [deps.UnsafeAtomics.weakdeps]
+    LLVM = "929cbde3-209d-540e-8aea-75f648917ca0"
+
+[[deps.VectorizationBase]]
+deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static", "StaticArrayInterface"]
+git-tree-sha1 = "d1d9a935a26c475ebffd54e9c7ad11627c43ea85"
+uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
+version = "0.21.72"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+version = "5.11.0+0"
diff --git a/benchmark/Project.toml b/benchmark/Project.toml
new file mode 100644
index 00000000..eb04cdfb
--- /dev/null
+++ b/benchmark/Project.toml
@@ -0,0 +1,4 @@
+[deps]
+BeamTracking = "8ef5c10a-4ca3-437f-8af5-b84d8af36df0"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
diff --git a/benchmark/rk4_kernel_benchmark.jl b/benchmark/rk4_kernel_benchmark.jl
new file mode 100644
index 00000000..35f2be21
--- /dev/null
+++ b/benchmark/rk4_kernel_benchmark.jl
@@ -0,0 +1,142 @@
+using BeamTracking
+using BeamTracking: Species, massof, chargeof, R_to_beta_gamma, R_to_pc, pc_to_R,
+                    RungeKuttaTracking, Bunch, STATE_ALIVE
+using StaticArrays
+using BenchmarkTools
+
+function setup_particle(pc=1e9)
+    species = Species("electron")
+    mc2 = massof(species)
+    p_over_q_ref = pc_to_R(species, pc)
+
+    beta_gamma_0 = R_to_beta_gamma(species, p_over_q_ref)
+    tilde_m = 1 / beta_gamma_0
+    gamsqr_0 = 1 + beta_gamma_0^2
+    beta_0 = beta_gamma_0 / sqrt(gamsqr_0)
+    charge = chargeof(species)
+    p0c = R_to_pc(species, p_over_q_ref)
+
+    return species, p_over_q_ref, beta_0, gamsqr_0, tilde_m, charge, p0c, mc2
+end
+
+function setup_solenoid_benchmark()
+    species, p_over_q_ref, beta_0, gamsqr_0, tilde_m, charge, p0c, mc2 = setup_particle(1e9)
+
+    bunch = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch.coords.v[1, BeamTracking.PXI] = 0.01
+
+    s_span = (0.0, 1.0)
+    ds_step = 0.01
+    g_bend = 0.0
+
+    # Solenoid field
+    Bz_physical = 0.01  # Tesla
+    Bz_normalized = Bz_physical / p_over_q_ref
+    mm = SVector(0)
+    kn = SVector(Bz_normalized)
+    ks = SVector(0.0)
+
+    return bunch, beta_0, tilde_m, charge, p0c, mc2, s_span, ds_step, g_bend, mm, kn, ks, p_over_q_ref
+end
+
+function reset_bunch!(bunch)
+    bunch.coords.v .= 0.0
+    bunch.coords.v[1, BeamTracking.PXI] = 0.01
+    bunch.coords.state[1] = STATE_ALIVE
+end
+
+# Setup
+bunch, beta_0, tilde_m, charge, p0c, mc2, s_span, ds_step, g_bend, mm, kn, ks, p_over_q_ref = setup_solenoid_benchmark()
+
+println("rk4_kernel! benchmark (1 particle)")
+println("=========================================")
+println("s_span: $s_span, ds_step: $ds_step")
+println("n_steps: $(Int(ceil((s_span[2] - s_span[1]) / ds_step)))")
+println()
+
+# Warmup
+reset_bunch!(bunch)
+RungeKuttaTracking.rk4_kernel!(1, bunch.coords, beta_0, tilde_m,
+                               charge, p0c, mc2, s_span, ds_step, g_bend,
+                               mm, kn, ks, p_over_q_ref)
+
+# Benchmark
+reset_bunch!(bunch)
+b = @benchmark begin
+    RungeKuttaTracking.rk4_kernel!(1, $bunch.coords, $beta_0, $tilde_m,
+                                   $charge, $p0c, $mc2, $s_span, $ds_step, $g_bend,
+                                   $mm, $kn, $ks, $p_over_q_ref)
+end setup=(reset_bunch!($bunch))
+
+display(b)
+println()
+
+# Multi-particle benchmark
+println("rk4_kernel! benchmark (1000 particles)")
+println("=========================================")
+
+function setup_multi_particle(n_particles)
+    species, p_over_q_ref, beta_0, gamsqr_0, tilde_m, charge, p0c, mc2 = setup_particle(1e9)
+
+    bunch = Bunch(randn(n_particles, 6) * 0.001, p_over_q_ref=p_over_q_ref, species=species)
+
+    s_span = (0.0, 1.0)
+    ds_step = 0.01
+    g_bend = 0.0
+
+    Bz_physical = 0.01
+    Bz_normalized = Bz_physical / p_over_q_ref
+    mm = SVector(0)
+    kn = SVector(Bz_normalized)
+    ks = SVector(0.0)
+
+    return bunch, beta_0, tilde_m, charge, p0c, mc2, s_span, ds_step, g_bend, mm, kn, ks, p_over_q_ref
+end
+
+function track_all_particles!(bunch, beta_0, tilde_m, charge, p0c, mc2,
+                              s_span, ds_step, g_bend, mm, kn, ks, p_over_q_ref)
+    n = size(bunch.coords.v, 1)
+    for i in 1:n
+        RungeKuttaTracking.rk4_kernel!(i, bunch.coords, beta_0, tilde_m,
+                                       charge, p0c, mc2, s_span, ds_step, g_bend,
+                                       mm, kn, ks, p_over_q_ref)
+    end
+end
+
+n_particles = 1000
+bunch_mp, beta_0_mp, tilde_m_mp, charge_mp, p0c_mp, mc2_mp,
+    s_span_mp, ds_step_mp, g_bend_mp, mm_mp, kn_mp, ks_mp, p_over_q_ref_mp = setup_multi_particle(n_particles)
+
+# Store initial state for reset
+v_init = copy(bunch_mp.coords.v)
+state_init = copy(bunch_mp.coords.state)
+
+function reset_multi!(bunch, v_init, state_init)
+    bunch.coords.v .= v_init
+    bunch.coords.state .= state_init
+end
+
+# Warmup
+track_all_particles!(bunch_mp, beta_0_mp, tilde_m_mp, charge_mp, p0c_mp, mc2_mp,
+                     s_span_mp, ds_step_mp, g_bend_mp, mm_mp, kn_mp, ks_mp, p_over_q_ref_mp)
+
+# Benchmark
+b_mp = @benchmark begin
+    track_all_particles!($bunch_mp, $beta_0_mp, $tilde_m_mp, $charge_mp,
+                         $p0c_mp, $mc2_mp, $s_span_mp, $ds_step_mp, $g_bend_mp,
+                         $mm_mp, $kn_mp, $ks_mp, $p_over_q_ref_mp)
+end setup=(reset_multi!($bunch_mp, $v_init, $state_init))
+
+display(b_mp)
+println()
+
+# Per-particle timing
+median_time_ns = median(b_mp).time
+println("\nPer-particle median time: $(median_time_ns / n_particles) ns")
+
+reset_multi!(bunch_mp, v_init, state_init)
+num_allocs_mp = @allocated track_all_particles!(bunch_mp, beta_0_mp, gamsqr_0_mp, tilde_m_mp, charge_mp,
+                                                 p0c_mp, mc2_mp, s_span_mp, ds_step_mp, g_bend_mp,
+                                                 mm_mp, kn_mp, ks_mp, p_over_q_ref_mp)
+println("Total allocations for $n_particles particles: $num_allocs_mp bytes")
+println("Per-particle allocations: $(num_allocs_mp / n_particles) bytes")
diff --git a/docs/src/runge_kutta.md b/docs/src/runge_kutta.md
new file mode 100644
index 00000000..6b4b194f
--- /dev/null
+++ b/docs/src/runge_kutta.md
@@ -0,0 +1,156 @@
+# Runge-Kutta Tracking
+
+The `RungeKutta` tracking method provides a 4th-order Runge-Kutta (RK4) integration scheme for tracking particles through arbitrary electromagnetic fields. This method is particularly useful for elements with complex fields that cannot be handled analytically.
+
+## Overview
+
+The Runge-Kutta tracking implementation uses a classical 4th-order Runge-Kutta method to numerically integrate the relativistic equations of motion for charged particles in electromagnetic fields. The implementation is optimized for GPU/SIMD compatibility using branchless operations and stack-allocated StaticArrays.
+
+## Configuration
+
+The `RungeKutta` tracking method can be configured with the following parameters:
+
+### Parameters
+
+- **`ds_step`**: The step size for integration (in meters). If not specified, defaults to `0.2` meters when neither parameter is set.
+
+- **`n_steps`** : The number of integration steps. If specified and positive, the step size is calculated as `h = L / n_steps` where `L` is the element length.
+
+**Important**: Only one of `ds_step` or `n_steps` should be specified. If both are set to positive values, an error will be raised. If neither is specified, `ds_step=0.2` is used by default.
+
+### Examples
+
+```julia
+# Use default step size (0.2 meters)
+ele.tracking_method = RungeKutta()
+
+# Specify step size explicitly
+ele.tracking_method = RungeKutta(ds_step=0.1)
+
+# Specify number of steps
+ele.tracking_method = RungeKutta(n_steps=50)
+```
+
+## Usage
+
+### Basic Usage with Beamlines.jl
+
+```julia
+using BeamTracking, Beamlines
+
+# Create an element with Runge-Kutta tracking
+ele = Quadrupole(L=1.0, k1=0.1)
+ele.tracking_method = RungeKutta(ds_step=0.1)
+bl = Beamline([ele], R_ref=1e6)
+
+# Create a bunch and track
+bunch = Bunch(zeros(100, 6), R_ref=1e6, species=Species("electron"))
+track!(bunch, ele)
+```
+The Runge-Kutta method works with all thick elements that have BMultipoleParams. Other thick elements are treated as drifts.
+
+### Low-Level Kernel Interface
+
+For advanced use cases, the `rk4_kernel!` function provides direct access to the integration routine:
+
+```julia
+rk4_kernel!(i, coords, beta_0, gamsqr_0, tilde_m, charge, p0c, mc2,
+            s_span, n_steps, g_bend, mm, kn, ks)
+```
+
+**Parameters:**
+- `i`: Particle index
+- `coords`: Particle coordinate array
+- `beta_0`, `gamsqr_0`, `tilde_m`: Reference particle parameters
+- `charge`: Particle charge (in units of elementary charge)
+- `p0c`, `mc2`: Reference momentum and rest mass energy
+- `s_span`: Integration range `(s_start, s_end)`
+- `n_steps`: Number of integration steps
+- `g_bend`: Curvature parameter (1/ρ for bends, 0 otherwise)
+- `mm`: StaticArray of multipole orders
+- `kn`: Normal multipole strengths
+- `ks`: Skew multipole strengths
+
+The multipole arrays define the magnetic field configuration. For example, a quadrupole with `k1=0.5` would use `mm=SA[2]`, `kn=SA[0.5]`, `ks=SA[0.0]`.
+
+## Physics
+
+The Runge-Kutta method integrates the relativistic equations of motion expressed in terms of arc length `s` as the independent variable. The state vector is:
+
+```math
+\mathbf{u} = [x, p_x, y, p_y, z, p_z].
+```
+
+The derivatives `du/ds` are calculated from the relativistic electromagnetism:
+
+### Position Derivatives
+
+```math
+\frac{dx}{ds} = v_x \frac{dt}{ds}, \quad \frac{dy}{ds} = v_y \frac{dt}{ds}
+```
+
+where `v_x, v_y` are the transverse velocity components and `dt/ds` accounts for the relationship between arc length and time.
+
+### Momentum Derivatives
+
+```math
+\frac{dp_x}{ds} = \frac{F_x}{p₀c} \frac{dt}{ds} + g_{bend} \frac{p_z}{p₀}, \quad \frac{dp_y}{ds} = \frac{F_y}{p₀c} \frac{dt}{ds}
+```
+
+where `F_x, F_y` are the Lorentz force components and `g_bend` is the curvature parameter (nonzero only in bend elements).
+
+### Longitudinal Coordinate
+
+```math
+\frac{dz}{ds} = \text{rel\_dir} \left(\frac{\beta}{\beta₀} - 1\right) + \text{rel\_dir} \frac{\sqrt{1-v_t²} - 1 - dh_{bend}}{\sqrt{1-v_t²}} + \frac{d\beta}{ds} \frac{z}{\beta}
+```
+At the moment `rel_dir` is hardcoded to be 1. In the future this could be extended to track particles going backwards.
+
+### Energy Deviation
+
+```math
+\frac{dp_z}{ds} = \frac{dp}{ds} / p₀c
+```
+
+where `dp/ds` is calculated from the work done by the electromagnetic fields.
+
+## Implementation Details
+
+### RK4 Algorithm
+
+The 4th-order Runge-Kutta method uses the standard four-stage scheme:
+
+1. **k₁**: Evaluate derivatives at current state `u(s)`
+2. **k₂**: Evaluate derivatives at `u(s) + (h/2) k₁`
+3. **k₃**: Evaluate derivatives at `u(s) + (h/2) k₂`
+4. **k₄**: Evaluate derivatives at `u(s) + h k₃`
+
+The final update is:
+
+```math
+u(s+h) = u(s) + \frac{h}{6}(k₁ + 2k₂ + 2k₃ + k₄)
+```
+
+
+### Particle Loss Condition
+
+The implementation includes detection of unphysical particle states:
+
+- **Unphysical momenta**: If the transverse velocity squared `v_t² ≥ 1`, the particle is marked as lost (`STATE_LOST_PZ`) and won't be tracked.
+
+### Multipole Field Calculation
+
+Electromagnetic fields are computed using the `multipole_em_field` function, which calculates field components from magnetic multipole parameters:
+
+```julia
+multipole_em_field(x, y, z, s, mm, kn, ks) -> (Ex, Ey, Ez, Bx, By, Bz)
+```
+
+**Parameters:**
+- **`mm`**: StaticArray of magnetic multipole orders (0=solenoid, 1=dipole, 2=quadrupole, 3=sextupole, etc.)
+- **`kn`**: Normal multipole strengths (normalized and not integrated)
+- **`ks`**: Skew multipole strengths (normalized and not integrated)
+
+**Field computation:**
+- For `m=0` (solenoid): Returns longitudinal field `Bz`
+- For `m≥1` (dipole, quadrupole, etc.): Computes transverse fields `Bx`, `By` using a Horner-like scheme for efficient polynomial evaluation
diff --git a/ext/BeamTrackingBeamlinesExt/BeamTrackingBeamlinesExt.jl b/ext/BeamTrackingBeamlinesExt/BeamTrackingBeamlinesExt.jl
index 5c674839..ffbd1857 100644
--- a/ext/BeamTrackingBeamlinesExt/BeamTrackingBeamlinesExt.jl
+++ b/ext/BeamTrackingBeamlinesExt/BeamTrackingBeamlinesExt.jl
@@ -43,5 +43,6 @@ include("unpack.jl")
 include("scibmadstandard.jl")
 include("exact.jl")
 include("yoshida.jl")
+include("rungekutta.jl")
 
 end
\ No newline at end of file
diff --git a/ext/BeamTrackingBeamlinesExt/rungekutta.jl b/ext/BeamTrackingBeamlinesExt/rungekutta.jl
new file mode 100644
index 00000000..a7336f3d
--- /dev/null
+++ b/ext/BeamTrackingBeamlinesExt/rungekutta.jl
@@ -0,0 +1,164 @@
+# =========== BYPASS UNPACKING FOR RUNGEKUTTA ============= #
+
+"""
+Specialized _track! for RungeKutta that bypasses the unpacking system.
+Gets field function directly from Beamlines.field_calc and passes the full element.
+"""
+function _track!(
+  coords::Coords,
+  bunch::Bunch,
+  ele::LineElement,
+  tm::RungeKutta,
+  scalar_params::Bool,
+  ramp_without_rf;
+  kwargs...
+)
+  # Get basic element properties (type-unstable unpacking)
+  L = float(ele.L)
+  ap = deval(ele.AlignmentParams)
+  bp = deval(ele.BendParams)
+  dp = deval(ele.ApertureParams)
+  patch = deval(ele.PatchParams)
+  bm = deval(ele.BMultipoleParams)
+  p_over_q_ref = bunch.p_over_q_ref
+
+  if scalar_params
+    L = scalarize(L)
+    ap = scalarize(ap)
+    bp = scalarize(bp)
+    bm = scalarize(bm)
+    pp = scalarize(pp)
+    dp = scalarize(dp)
+    mp = scalarize(mp)
+    rp = scalarize(rp)
+    lp = scalarize(lp)
+    p_over_q_ref = scalarize(p_over_q_ref)
+  end
+  
+  # Function barrier
+  runge_kutta_universal!(coords, tm, ramp_without_rf, bunch, L, p_over_q_ref, ap, bp, dp, patch, bm; kwargs...)
+end
+
+# Step 2: Type-stable computation -----------------------------------------
+function runge_kutta_universal!(
+  coords,
+  tm,
+  ramp_without_rf,
+  bunch,
+  L,
+  p_over_q_ref,
+  alignmentparams,
+  bendparams,
+  apertureparams,
+  patchparams,
+  bmultipoleparams;
+  kwargs...
+)
+  species = bunch.species
+  # Setup reference state
+  beta_gamma_ref = R_to_beta_gamma(species, p_over_q_ref)
+  kc = KernelChain(Val{6}(), RefState(bunch.t_ref, beta_gamma_ref))
+
+  # Evolve time through whole element
+  bunch.t_ref += L/beta_gamma_to_v(beta_gamma_ref)
+
+  # Handle reference momentum ramping
+  if p_over_q_ref isa TimeDependentParam
+    p_over_q_ref_initial = p_over_q_ref
+    p_over_q_ref_final = p_over_q_ref(bunch.t_ref)
+    if !(p_over_q_ref_initial ≈ p_over_q_ref_final)
+      kc = push(kc, KernelCall(BeamTracking.update_P0!, (p_over_q_ref_initial, p_over_q_ref_final, ramp_without_rf)))
+      setfield!(bunch, :p_over_q_ref, p_over_q_ref_final)
+    end
+  end
+
+  # Error conditions
+  if L <= 0.0
+    error("RungeKutta tracking does not support zero-length elements") 
+  end
+  if isactive(patchparams)
+    error("RungeKutta tracking does not support patch elements")
+  end
+
+  # Entrance aperture and alignment
+  if isactive(alignmentparams)
+    if isactive(apertureparams)
+      if apertureparams.aperture_shifts_with_body
+        kc = push(kc, @inline(alignment(tm, bunch, alignmentparams, bendparams, L, true)))
+        kc = push(kc, @inline(aperture(tm, bunch, apertureparams, true)))
+      else
+        kc = push(kc, @inline(aperture(tm, bunch, apertureparams, true)))
+        kc = push(kc, @inline(alignment(tm, bunch, alignmentparams, bendparams, L, true)))
+      end
+    else
+      kc = push(kc, @inline(alignment(tm, bunch, alignmentparams, bendparams, L, true)))
+    end
+  elseif isactive(apertureparams)
+    kc = push(kc, @inline(aperture(tm, bunch, apertureparams, true)))
+  end
+
+  # Setup physics parameters
+  p_over_q_ref = bunch.p_over_q_ref
+  tilde_m, _, beta_0 = BeamTracking.drift_params(species, p_over_q_ref)
+  charge = chargeof(species)
+  p0c = BeamTracking.R_to_pc(species, p_over_q_ref)
+  mc2 = massof(species)
+
+  # Determine step size to use
+  if tm.ds_step > 0
+    ds_step = tm.ds_step
+  elseif tm.n_steps > 0
+    ds_step = L / tm.n_steps
+  else
+    ds_step = BeamTracking.DEFAULT_RK4_DS_STEP
+  end
+
+  s_span = (0.0, L)
+
+  # Get curvature from BendParams if present
+  g_bend = isactive(bendparams) ? bendparams.g_ref : 0.0
+
+  # Extract multipole parameters
+  if isactive(bmultipoleparams)
+    mm = bmultipoleparams.order
+    kn, ks = get_strengths(bmultipoleparams, L, p_over_q_ref)
+  else
+    # Default to drift
+    mm = SVector{0, Int}()
+    kn = SVector{0, typeof(L)}()
+    ks = SVector{0, typeof(L)}()
+  end
+
+  # Build RK4 kernel call
+  params = (beta_0, tilde_m, charge, p0c, mc2, s_span, ds_step, g_bend, mm, kn, ks, p_over_q_ref)
+  kc = push(kc, KernelCall(BeamTracking.RungeKuttaTracking.rk4_kernel!, params))
+
+  # Exit aperture and alignment
+  if isactive(alignmentparams)
+    if isactive(apertureparams)
+      if apertureparams.aperture_shifts_with_body
+        kc = push(kc, @inline(aperture(tm, bunch, apertureparams, false)))
+        kc = push(kc, @inline(alignment(tm, bunch, alignmentparams, bendparams, L, false)))
+      else
+        kc = push(kc, @inline(alignment(tm, bunch, alignmentparams, bendparams, L, false)))
+        kc = push(kc, @inline(aperture(tm, bunch, apertureparams, false)))
+      end
+    else
+      kc = push(kc, @inline(alignment(tm, bunch, alignmentparams, bendparams, L, false)))
+    end
+  elseif isactive(apertureparams)
+    kc = push(kc, @inline(aperture(tm, bunch, apertureparams, false)))
+  end
+
+  # Launch kernels
+  @noinline launch!(coords, kc; kwargs...)
+  return nothing
+end
+
+# =========== ALIGNMENT AND APERTURE ============= #
+
+@inline alignment(tm::RungeKutta, bunch, alignmentparams, bendparams, L, entering) =
+  alignment(Exact(), bunch, alignmentparams, bendparams, L, entering)
+
+@inline aperture(tm::RungeKutta, bunch, apertureparams, entering) =
+  aperture(Exact(), bunch, apertureparams, entering)
diff --git a/src/BeamTracking.jl b/src/BeamTracking.jl
index cb00b020..da4a456a 100644
--- a/src/BeamTracking.jl
+++ b/src/BeamTracking.jl
@@ -18,7 +18,7 @@ using KernelAbstractions
 import GTPSA: sincu, sinhcu
 
 export Bunch, State, ParticleView, Time, TimeDependentParam
-export Yoshida, Yoshida, MatrixKick, BendKick, SolenoidKick, DriftKick, Exact
+export Yoshida, Yoshida, MatrixKick, BendKick, SolenoidKick, DriftKick, Exact, RungeKutta
 export track!
 
 
@@ -49,6 +49,7 @@ include("kernels/spin.jl")
 include("kernels/transforms.jl")
 include("kernels/yoshida.jl")
 
+include("modules/RungeKuttaTracking.jl") #; TRACKING_METHOD(::RungeKuttaTracking) = RungeKutta
 
 # Empty tracking method to be imported+implemented by package extensions
 function track! end
diff --git a/src/modules/RungeKuttaTracking.jl b/src/modules/RungeKuttaTracking.jl
new file mode 100644
index 00000000..c2a29be1
--- /dev/null
+++ b/src/modules/RungeKuttaTracking.jl
@@ -0,0 +1,269 @@
+"""
+  RungeKuttaTracking
+
+Module implementing particle tracking through arbitrary electromagnetic fields using a 4th order Runge-Kutta method.
+"""
+module RungeKuttaTracking
+using ..BeamTracking, ..StaticArrays
+using ..BeamTracking: @makekernel, Coords
+using ..BeamTracking: XI, PXI, YI, PYI, ZI, PZI, STATE_ALIVE, STATE_LOST_PZ
+using ..BeamTracking: C_LIGHT, E_CHARGE, vifelse, normalized_field
+
+
+"""
+  multipole_em_field(x, y, z, s, mm, kn, ks, p_over_q_ref)
+
+Compute EM field from multipole moments for RK4 tracking.
+Handles ALL multipole orders:
+- m=0: solenoid (longitudinal Bz)
+- m=1: dipole (transverse By, Bx)
+- m≥2: higher-order multipoles (quadrupole, sextupole, etc.)
+
+Returns (Ex, Ey, Ez, Bx, By, Bz) in physical units (Tesla for B, V/m for E) where:
+- Bx, By: transverse field from all orders except m=0
+- Bz: longitudinal field from m=0 term if present
+- Ex, Ey, Ez: zero (static magnetic elements only)
+"""
+@inline function multipole_em_field(x, y, z, s, mm::SVector{0}, kn, ks, p_over_q_ref)
+  return (zero(x), zero(x), zero(x), zero(x), zero(x), zero(x))
+end
+
+@inline function multipole_em_field(x, y, z, s, mm::SVector{N}, kn, ks, p_over_q_ref) where N
+  bx, by = normalized_field(mm, kn, ks, x, y, 0)
+  is_solenoid = (mm[1] == 0)
+  bz = vifelse(is_solenoid, kn[1], zero(x))
+
+  # Convert from normalized (field/Bρ) to physical units (Tesla)
+  return (zero(x), zero(x), zero(x), bx * p_over_q_ref, by * p_over_q_ref, bz * p_over_q_ref)
+end
+
+"""
+  kick_vector(x, px, y, py, z, pz, s, Ex, Ey, Ez, Bx, By, Bz,
+        charge, tilde_m, beta_0, g_bend, p0c, mc2)
+
+Calculate the derivative vector du/ds for relativistic particle tracking.
+Returns an SVector{6} containing [dx/ds, dpx/ds, dy/ds, dpy/ds, dz/ds, dpz/ds].
+
+Uses branchless operations for GPU/SIMD compatibility. For unphysical momenta,
+returns zero derivatives (caller should mark particle as lost).
+
+# Arguments
+- `x, px, y, py, z, pz`: State vector components
+- `s`: Arc length position
+- `Ex, Ey, Ez`: Electric field components (V/m)
+- `Bx, By, Bz`: Magnetic field components (T)
+- `charge`: Particle charge in units of e
+- `tilde_m`: Normalized mass mc²/(p₀c)
+- `beta_0`: Reference velocity β₀ = v₀/c
+- `g_bend`: Curvature (0 for drift, 1/ρ for bends)
+- `p0c`: Reference momentum × c (eV)
+- `mc2`: Rest mass energy (eV)
+"""
+@inline function kick_vector(x, px, y, py, z, pz, s, Ex, Ey, Ez, Bx, By, Bz,
+                charge, tilde_m, beta_0, g_bend, p0c, mc2)
+  # Relative momentum
+  rel_p = 1 + pz
+
+  # Transverse velocity components (normalized)
+  vt_x = px / rel_p
+  vt_y = py / rel_p
+  vt2 = vt_x^2 + vt_y^2
+
+  # Check for unphysical momenta (branchless)
+  vt2_1 = one(vt2)
+  good_momenta = (vt2 < vt2_1)
+  vt2_safe = vifelse(good_momenta, vt2, zero(vt2))
+
+  # Particle beta and velocity
+  rel_p2 = rel_p^2
+  inv_gamma_v = sqrt(rel_p2 + tilde_m^2)
+  beta = rel_p / inv_gamma_v
+  
+  inv_beta_c = 1 / (beta * C_LIGHT)
+
+  # Longitudinal velocity component
+  rel_dir = 1  # +1 for forward tracking
+  vz_norm = sqrt(1 - vt2_safe) * rel_dir
+  vx = beta * C_LIGHT * vt_x
+  vy = beta * C_LIGHT * vt_y
+  vz = beta * C_LIGHT * vz_norm
+
+  # Lorentz force: F = q*(E + v×B)
+  E_force_x = charge * Ex
+  E_force_y = charge * Ey
+  E_force_z = charge * Ez
+  B_force_x = charge * (vy*Bz - vz*By)
+  B_force_y = charge * (vz*Bx - vx*Bz)
+  B_force_z = charge * (vx*By - vy*Bx)
+
+  # Time derivative w.r.t. arc length
+  dh_bend = x * g_bend  # Longitudinal distance deviation
+  abs_vz = abs(vz)
+  abs_vz_safe = vifelse(good_momenta, abs_vz, one(abs_vz))  # Avoid division by zero
+  dt_ds = rel_dir * (1 + dh_bend) / abs_vz_safe
+
+  # Longitudinal momentum (normalized)
+  pz_p0 = rel_p * rel_dir * abs_vz * inv_beta_c
+
+  # Energy derivative: dp/ds = (F · v) * dt/ds * inv_beta_c
+  F_dot_v = E_force_x*vx + E_force_y*vy + E_force_z*vz
+  dp_ds = F_dot_v * dt_ds * inv_beta_c
+
+  # Total energy for dbeta_ds calculation
+  e_tot = p0c * rel_p / beta
+  dbeta_ds = mc2^2 * dp_ds * C_LIGHT / e_tot^3
+
+  # Position derivatives: dr/ds = v * dt/ds
+  dx_ds = vx * dt_ds
+  dy_ds = vy * dt_ds
+
+  # Momentum derivatives: dp_i/ds = F_i * dt/ds / p0c + corrections
+  p0 = p0c / C_LIGHT
+  dpx_ds = (E_force_x + B_force_x) * dt_ds / p0 + g_bend * pz_p0
+  dpy_ds = (E_force_y + B_force_y) * dt_ds / p0
+
+  # Longitudinal coordinate z derivative
+  sqrt_1mvt2 = sqrt(1 - vt2_safe)
+  dz_ds = rel_dir * (beta / beta_0 - 1) + rel_dir * (sqrt_1mvt2 - 1 - dh_bend) / sqrt_1mvt2 + dbeta_ds * z / beta
+
+  # Energy deviation derivative
+  dpz_ds = dp_ds / p0
+
+  # Return zero derivatives if momenta are unphysical (branchless)
+  zero_deriv = zero(dx_ds)
+  return SVector(
+    vifelse(good_momenta, dx_ds, zero_deriv),
+    vifelse(good_momenta, dpx_ds, zero_deriv),
+    vifelse(good_momenta, dy_ds, zero_deriv),
+    vifelse(good_momenta, dpy_ds, zero_deriv),
+    vifelse(good_momenta, dz_ds, zero_deriv),
+    vifelse(good_momenta, dpz_ds, zero_deriv)
+  )
+end
+
+"""
+  rk4_step!(coords, i, s, h, mm, kn, ks, tracking_params)
+
+Perform a single RK4 step for particle i, updating coordinates in-place.
+Only updates state if particle is alive.
+
+# Arguments
+- `coords`: Coordinates structure
+- `i`: Particle index
+- `s`: Current arc length
+- `h`: Step size
+- `mm`: Multipole orders (StaticArray)
+- `kn`: Normal multipole strengths (StaticArray)
+- `ks`: Skew multipole strengths (StaticArray)
+- `charge`: Particle charge in units of e
+- `tilde_m`: Normalized mass mc²/(p₀c)
+- `beta_0`: Reference velocity β₀ = v₀/c
+- `g_bend`: Curvature (0 for drift, 1/ρ for bends)
+- `p0c`: Reference momentum × c (eV)
+- `mc2`: Rest mass energy (eV)
+- `p_over_q_ref`: Reference magnetic rigidity Bρ = p₀c/(c·charge)
+"""
+@inline function rk4_step!(coords, i, s, h, mm, kn, ks, charge, tilde_m, beta_0, g_bend, p0c, mc2, p_over_q_ref)
+  # Check if particle is alive
+  alive = (coords.state[i] == STATE_ALIVE)
+  
+  # Extract current particle
+  v = coords.v
+  x = v[i, XI]
+  px = v[i, PXI]
+  y = v[i, YI]
+  py = v[i, PYI]
+  z = v[i, ZI]
+  pz = v[i, PZI]
+
+  # k1 = f(u, s)
+  Ex, Ey, Ez, Bx, By, Bz = multipole_em_field(x, y, z, s, mm, kn, ks, p_over_q_ref)
+  k1 = kick_vector(x, px, y, py, z, pz, s, Ex, Ey, Ez, Bx, By, Bz,
+                charge, tilde_m, beta_0, g_bend, p0c, mc2)
+
+  # k2 = f(u + h/2 * k1, s + h/2)
+  h2 = h / 2
+  x2 = x + h2 * k1[1]
+  px2 = px + h2 * k1[2]
+  y2 = y + h2 * k1[3]
+  py2 = py + h2 * k1[4]
+  z2 = z + h2 * k1[5]
+  pz2 = pz + h2 * k1[6]
+  Ex, Ey, Ez, Bx, By, Bz = multipole_em_field(x2, y2, z2, s + h2, mm, kn, ks, p_over_q_ref)
+  k2 = kick_vector(x2, px2, y2, py2, z2, pz2, s + h2, Ex, Ey, Ez, Bx, By, Bz,
+                charge, tilde_m, beta_0, g_bend, p0c, mc2)
+
+  # k3 = f(u + h/2 * k2, s + h/2)
+  x3 = x + h2 * k2[1]
+  px3 = px + h2 * k2[2]
+  y3 = y + h2 * k2[3]
+  py3 = py + h2 * k2[4]
+  z3 = z + h2 * k2[5]
+  pz3 = pz + h2 * k2[6]
+  Ex, Ey, Ez, Bx, By, Bz = multipole_em_field(x3, y3, z3, s + h2, mm, kn, ks, p_over_q_ref)
+  k3 = kick_vector(x3, px3, y3, py3, z3, pz3, s + h2, Ex, Ey, Ez, Bx, By, Bz,
+                charge, tilde_m, beta_0, g_bend, p0c, mc2)
+
+  # k4 = f(u + h * k3, s + h)
+  x4 = x + h * k3[1]
+  px4 = px + h * k3[2]
+  y4 = y + h * k3[3]
+  py4 = py + h * k3[4]
+  z4 = z + h * k3[5]
+  pz4 = pz + h * k3[6]
+  Ex, Ey, Ez, Bx, By, Bz = multipole_em_field(x4, y4, z4, s + h, mm, kn, ks, p_over_q_ref)
+  k4 = kick_vector(x4, px4, y4, py4, z4, pz4, s + h, Ex, Ey, Ez, Bx, By, Bz,
+                charge, tilde_m, beta_0, g_bend, p0c, mc2)
+
+  # Update state: u += h/6 * (k1 + 2*k2 + 2*k3 + k4)
+  # Only update if particle is alive
+  h6 = h / 6
+  v[i, XI] = vifelse(alive, x + h6 * (k1[1] + 2*k2[1] + 2*k3[1] + k4[1]), v[i, XI])
+  v[i, PXI] = vifelse(alive, px + h6 * (k1[2] + 2*k2[2] + 2*k3[2] + k4[2]), v[i, PXI])
+  v[i, YI] = vifelse(alive, y + h6 * (k1[3] + 2*k2[3] + 2*k3[3] + k4[3]), v[i, YI])
+  v[i, PYI] = vifelse(alive, py + h6 * (k1[4] + 2*k2[4] + 2*k3[4] + k4[4]), v[i, PYI])
+  v[i, ZI] = vifelse(alive, z + h6 * (k1[5] + 2*k2[5] + 2*k3[5] + k4[5]), v[i, ZI])
+  v[i, PZI] = vifelse(alive, pz + h6 * (k1[6] + 2*k2[6] + 2*k3[6] + k4[6]), v[i, PZI])
+end
+
+"""
+  rk4_kernel!(i, coords, beta_0, tilde_m, charge, p0c, mc2,
+              s_span, ds_step, g_bend, mm, kn, ks, p_over_q_ref)
+
+Kernelized RK4 tracking through multipole fields.
+Compatible with @makekernel and the package's kernel architecture.
+
+The electromagnetic field is computed from multipole moments (mm, kn, ks) using
+the multipole_em_field function.
+"""
+@makekernel function rk4_kernel!(i, coords::Coords, beta_0, tilde_m,
+                                charge, p0c, mc2, s_span, ds_step, g_bend, mm, kn, ks, p_over_q_ref)
+  s_start = s_span[1]
+  s_end = s_span[2]
+  s = s_start
+
+  v = coords.v
+  
+  # Calculate number of steps for deterministic iteration
+  total_distance = s_end - s_start
+  n_steps = ceil(Int, total_distance / ds_step)
+  
+  for step in 1:n_steps
+    remaining = s_end - s
+    h = min(ds_step, remaining)
+
+    # Chck if particle is lost
+    rel_p = 1 + v[i, PZI]
+    inv_rel_p = 1 / rel_p
+    vt2 = (v[i, PXI] * inv_rel_p)^2 + (v[i, PYI] * inv_rel_p)^2
+    alive = (coords.state[i] == STATE_ALIVE)
+    # Mark particle as lost
+    coords.state[i] = vifelse((vt2 >= 1) & alive, STATE_LOST_PZ, coords.state[i])
+
+    # Perform RK4 step (check for alive status is now inside rk4_step!)
+    rk4_step!(coords, i, s, h, mm, kn, ks, charge, tilde_m, beta_0, g_bend, p0c, mc2, p_over_q_ref)
+    s += h
+  end
+end
+
+end
diff --git a/src/tracking_methods.jl b/src/tracking_methods.jl
index f2d75f60..52830ccb 100644
--- a/src/tracking_methods.jl
+++ b/src/tracking_methods.jl
@@ -40,4 +40,35 @@ end
 
 
 # ========== Exact ===========================
-struct Exact end
\ No newline at end of file
+struct Exact end
+
+# ========== Explicit RK4 Tracking ==========
+struct RungeKutta
+  ds_step::Float64
+  n_steps::Int
+end
+
+DEFAULT_RK4_DS_STEP = 0.2
+function RungeKutta(; ds_step::Union{Float64, Nothing}=nothing, n_steps::Union{Int, Nothing}=nothing)
+  # Get actual values (use provided or sentinel)
+  _ds_step = ds_step === nothing ? -1.0 : ds_step
+  _n_steps = n_steps === nothing ? -1 : n_steps
+  
+  # Error if both are explicitly set to positive values
+  if _ds_step > 0 && _n_steps > 0
+    error("Only one of ds_step or n_steps should be specified")
+  end
+  
+  # If user sets n_steps (and it's positive), set ds_step to negative
+  if _n_steps > 0
+    return RungeKutta(-1.0, _n_steps)
+  end
+  
+  # If user sets ds_step (and it's positive), set n_steps to -1
+  if _ds_step > 0
+    return RungeKutta(_ds_step, -1)
+  end
+  
+  # Fallback: use defaults if both are negative/not set
+  return RungeKutta(DEFAULT_RK4_DS_STEP, -1)
+end
\ No newline at end of file
diff --git a/test/RungeKuttaTracking_test.jl b/test/RungeKuttaTracking_test.jl
new file mode 100644
index 00000000..476c3def
--- /dev/null
+++ b/test/RungeKuttaTracking_test.jl
@@ -0,0 +1,265 @@
+@testset "RungeKuttaTracking" begin
+  using BeamTracking
+  using BeamTracking: Species, massof, chargeof, R_to_beta_gamma, R_to_pc, pc_to_R,
+                      RungeKuttaTracking, Bunch, STATE_ALIVE, STATE_LOST_PZ, E_CHARGE, C_LIGHT
+  using StaticArrays
+
+  # Helper function to setup tracking parameters
+  function setup_particle(pc=1e9)  # pc in eV, default corresponds to 1 GeV
+    species = Species("electron")
+    mc2 = massof(species)  # eV
+    p_over_q_ref = pc_to_R(species, pc)
+
+    # Calculate tracking parameters
+    beta_gamma_0 = R_to_beta_gamma(species, p_over_q_ref)
+    tilde_m = 1 / beta_gamma_0
+    gamsqr_0 = 1 + beta_gamma_0^2
+    beta_0 = beta_gamma_0 / sqrt(gamsqr_0)
+    charge = chargeof(species)
+    p0c = R_to_pc(species, p_over_q_ref)
+
+    return species, p_over_q_ref, beta_0, gamsqr_0, tilde_m, charge, p0c, mc2
+  end
+
+  @testset "RungeKutta constructor" begin
+    using BeamTracking: RungeKutta
+
+    # Test default constructor (no arguments)
+    rk_default = RungeKutta()
+    @test rk_default.ds_step == 0.2
+    @test rk_default.n_steps == -1
+
+    # Test constructor with ds_step only
+    rk_ds = RungeKutta(ds_step=0.1)
+    @test rk_ds.ds_step == 0.1
+    @test rk_ds.n_steps == -1
+
+    # Test constructor with n_steps only
+    rk_ns = RungeKutta(n_steps=50)
+    @test rk_ns.ds_step == -1.0
+    @test rk_ns.n_steps == 50
+
+    # Test constructor with both ds_step and n_steps (should error)
+    @test_throws ErrorException RungeKutta(ds_step=0.1, n_steps=50)
+
+    # Test constructor with explicit nothing values (should use defaults)
+    rk_nothing = RungeKutta(ds_step=nothing, n_steps=nothing)
+    @test rk_nothing.ds_step == 0.2
+    @test rk_nothing.n_steps == -1
+  end
+
+  @testset "Pure drift" begin
+    species, p_over_q_ref, beta_0, gamsqr_0, tilde_m, charge, p0c, mc2 = setup_particle()
+
+    # Create bunch with small transverse momentum
+    bunch = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch.coords.v[1, BeamTracking.PXI] = 0.01
+
+    s_span = (0.0, 1.0)
+    ds_step = 0.01
+    g_bend = 0.0
+    
+    # Empty multipole vectors for drift
+    mm = SVector{0, Int}()
+    kn = SVector{0, Float64}()
+    ks = SVector{0, Float64}()
+
+    RungeKuttaTracking.rk4_kernel!(1, bunch.coords, beta_0, tilde_m,
+                                   charge, p0c, mc2, s_span, ds_step, g_bend,
+                                   mm, kn, ks, p_over_q_ref)
+
+    # Regression test
+    solution = [0.0100005  0.01  0.0  0.0  -5.00038e-5  0.0]
+    @test isapprox(bunch.coords.v, solution, rtol=1e-6)
+    @test bunch.coords.state[1] == STATE_ALIVE
+  end
+
+  @testset "Solenoid" begin
+    species, p_over_q_ref, beta_0, gamsqr_0, tilde_m, charge, p0c, mc2 = setup_particle(1e9)
+
+    bunch = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch.coords.v[1, BeamTracking.PXI] = 0.01
+
+    s_span = (0.0, 1.0)
+    ds_step = 0.01
+    g_bend = 0.0
+    
+    # Solenoid field
+    Bz_physical = 0.01  # Tesla
+    Bz_normalized = Bz_physical / p_over_q_ref
+    mm = SVector(0)  # Solenoid (m=0)
+    kn = SVector(Bz_normalized)
+    ks = SVector(0.0)
+
+    RungeKuttaTracking.rk4_kernel!(1, bunch.coords, beta_0, tilde_m,
+                                   charge, p0c, mc2, s_span, ds_step, g_bend,
+                                   mm, kn, ks, p_over_q_ref)
+
+    # In uniform B-field, particle should follow circular path
+    # Total transverse momentum should be conserved
+    pt2 = bunch.coords.v[1, 2]^2 + bunch.coords.v[1, 4]^2
+    @test isapprox(pt2, 0.01^2, rtol=1e-4)
+    # Regression test
+    solution = [0.010000485056009705 0.009999955057780502 1.4991110783291216e-5 2.9980699961334158e-5 -5.000375031233078e-5 0.0]
+    @test isapprox(bunch.coords.v, solution, rtol=1e-6)
+    @test bunch.coords.state[1] == STATE_ALIVE
+  end
+
+  @testset "Dipole" begin
+    species, p_over_q_ref, beta_0, gamsqr_0, tilde_m, charge, p0c, mc2 = setup_particle(1e9)
+
+    bunch = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch.coords.v[1, BeamTracking.PXI] = 0.01
+
+    s_span = (0.0, 1.0)
+    ds_step = 0.01
+    g_bend = 0.0
+    
+    # Dipole field
+    By_physical = 0.01  # Tesla
+    By_normalized = By_physical / p_over_q_ref
+    mm = SVector(1)  # Dipole (m=1)
+    kn = SVector(By_normalized)
+    ks = SVector(0.0)
+
+    RungeKuttaTracking.rk4_kernel!(1, bunch.coords, beta_0, tilde_m,
+                                   charge, p0c, mc2, s_span, ds_step, g_bend,
+                                   mm, kn, ks, p_over_q_ref)
+
+    # Regression test
+    solution = [0.011499735519796054 0.012997924579999955 0.0 0.0 -6.649432859025015e-5 0.0]
+    @test isapprox(bunch.coords.v, solution, rtol=1e-6)
+    @test bunch.coords.state[1] == STATE_ALIVE
+  end
+
+  @testset "Particle loss detection" begin
+    species, p_over_q_ref, beta_0, gamsqr_0, tilde_m, charge, p0c, mc2 = setup_particle(1e9)
+
+    bunch = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch.coords.v[1, BeamTracking.PXI] = 1.5 # Unphysical initial momentum
+
+    s_span = (0.0, 1.0)
+    ds_step = 0.1  # 10 cm step size
+    g_bend = 0.0
+    
+    # Empty multipole vectors for drift
+    mm = SVector{0, Int}()
+    kn = SVector{0, Float64}()
+    ks = SVector{0, Float64}()
+
+    RungeKuttaTracking.rk4_kernel!(1, bunch.coords, beta_0, tilde_m,
+                                   charge, p0c, mc2, s_span, ds_step, g_bend,
+                                   mm, kn, ks, p_over_q_ref)
+
+    # Particle should not track
+    solution = [0.0  1.5  0.0  0.0  0.0  0.0]
+    @test isapprox(bunch.coords.v, solution, rtol=1e-6)
+    @test bunch.coords.state[1] == STATE_LOST_PZ
+  end
+
+  @testset "Convergence test" begin
+    species, p_over_q_ref, beta_0, gamsqr_0, tilde_m, charge, p0c, mc2 = setup_particle(1e9)
+
+    bunch1 = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch2 = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch1.coords.v[1, BeamTracking.PXI] = 0.01
+    bunch2.coords.v[1, BeamTracking.PXI] = 0.01
+
+    s_span = (0.0, 1.0)
+    g_bend = 0.0
+    
+    # Empty multipole vectors for drift
+    mm = SVector{0, Int}()
+    kn = SVector{0, Float64}()
+    ks = SVector{0, Float64}()
+
+    # Track with different step sizes
+    RungeKuttaTracking.rk4_kernel!(1, bunch1.coords, beta_0, tilde_m,
+                                   charge, p0c, mc2, s_span, 0.1, g_bend,
+                                   mm, kn, ks, p_over_q_ref)
+    RungeKuttaTracking.rk4_kernel!(1, bunch2.coords, beta_0, tilde_m,
+                                   charge, p0c, mc2, s_span, 0.05, g_bend,
+                                   mm, kn, ks, p_over_q_ref)
+
+    # Results should be identical
+    @test isapprox(bunch1.coords.v, bunch2.coords.v, rtol=1e-2)
+  end
+
+  @testset "Beamlines integration - Drift" begin
+    using Beamlines
+
+    species, p_over_q_ref, _, _, _, _, _, _ = setup_particle()
+    bunch = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch.coords.v[1, BeamTracking.PXI] = 0.01
+
+    drift_ele = Drift(L=1.0)
+    drift_ele.tracking_method = RungeKutta()
+
+    track!(bunch, drift_ele)
+
+    # Regression test
+    solution = [0.0100005  0.01  0.0  0.0  -5.00038e-5  0.0]
+    @test isapprox(bunch.coords.v, solution, rtol=1e-6)
+  end
+
+  @testset "Beamlines integration - SBend" begin
+    using Beamlines
+
+    species, p_over_q_ref, _, _, _, _, _, _ = setup_particle()
+    bunch = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch.coords.v[1, BeamTracking.PXI] = 0.01
+
+    sbend_ele = SBend(L=1.0, angle=pi/132)
+    sbend_ele.tracking_method = RungeKutta()
+
+    track!(bunch, sbend_ele)
+
+    # Regression test
+    solution = [0.010000150630002367 0.009995978032305387 0.0 0.0 -0.00016899908120890584 0.0]
+    @test isapprox(bunch.coords.v, solution, rtol=1e-6)
+  end
+
+  @testset "RungeKutta with different step configurations" begin
+    using Beamlines
+
+    species, p_over_q_ref, _, _, _, _, _, _ = setup_particle()
+
+    # Test with ds_step
+    drift_ds = Drift(L=1.0)
+    drift_ds.tracking_method = RungeKutta(ds_step=0.1)
+    bunch_ds = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch_ds.coords.v[1, BeamTracking.PXI] = 0.01
+    track!(bunch_ds, drift_ds)
+
+    # Test with n_steps
+    drift_ns = Drift(L=1.0)
+    drift_ns.tracking_method = RungeKutta(n_steps=10)
+    bunch_ns = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    bunch_ns.coords.v[1, BeamTracking.PXI] = 0.01
+    track!(bunch_ns, drift_ns)
+
+    # Both should give the same results
+    @test isapprox(bunch_ds.coords.v, bunch_ns.coords.v, rtol=1e-2)
+  end
+
+  @testset "Zero-length elements" begin
+    using Beamlines
+
+    species, p_over_q_ref, _, _, _, _, _, _ = setup_particle()
+
+    # Test zero-length drift should throw an error
+    drift_zero = Drift(L=0.0)
+    drift_zero.tracking_method = RungeKutta()
+    bunch_drift = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    
+    @test_throws ErrorException track!(bunch_drift, drift_zero)
+
+    # Test negative length should also throw an error
+    drift_negative = Drift(L=-0.1)
+    drift_negative.tracking_method = RungeKutta()
+    bunch_negative = Bunch(zeros(1, 6), p_over_q_ref=p_over_q_ref, species=species)
+    
+    @test_throws ErrorException track!(bunch_negative, drift_negative)
+  end
+
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 7fefaea6..5b2ad11f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -200,4 +200,5 @@ include("alignment_tracking_test.jl")
 include("aperture_tracking_test.jl")
 include("ExactTracking_test.jl")
 include("IntegrationTracking_test.jl")
-include("time_test.jl")
\ No newline at end of file
+include("time_test.jl")
+include("RungeKuttaTracking_test.jl")
\ No newline at end of file