NonlinearNormalForm

Generic symplectic skew symmetric S matrix (size inferred from other matrix in operations) using SkewLinearAlgebra's JMatrix

Generic symplectic skew symmetric S matrix (size inferred from other matrix in operations) using SkewLinearAlgebra's JMatrix

TaylorMap{V0,V,Q,S}

Abstract type for TPSAMap and DAMap used for normal form analysis. DAMaps have a coordinate system chosen so that the expansion point is always around zero, e.g. Δv = v, while TPSAMaps can have any coordinate system/expansion point, but therefore will accrue truncation error when trying to compose two TPSAMaps with differing expansion points. For normal form analysis of periodic maps, using a DAMap ensures no truncation error up to the chosen truncation order.

Any truncated power series (TPS) type supported by TPSAInterface.jl is allowed for use in a TaylorMap. Henceforth we will generically refer to this type as a TPS

Fields

• v0::V0 – Reference orbit. The entrance coordinates of the map as scalars, or equivalently the Taylor map expansion point.
• v::V – Orbital ray as truncated power series, expansion around v0, with scalar part equal to EXIT coordinates of map
• q::Q – Quaternion as truncated power series if spin is included, else nothing
• s::S – Matrix of the envelope for stochastic kicks as scalars if included, else nothing

Type Requirements

• V0 <: AbstractVector{<:Number} where ismutabletype(V0) == true
• V <: AbstractVector{<:TPS} where TPSAInterface.numtype(V) == eltype(V0)
• Q <: Union{Quaternion{<:TPS},Nothing} and if Q != Nothing then eltype(Q) == eltype(V)
• S <: Union{AbstractMatrix{<:Number},Nothing} where S != Nothing then eltype(S) == TPSAInterface.numtype(eltype(V)) AND ismutabletype(S) == true

Because the TPS type is mutable and TPSAInterface.jl provides in-place functions for modifying TPSs, at the lowest level, all operations on TaylorMaps are in-place for performance. Therefore, the v and q arrays which contain TPSs may be immutable, e.g. the orbital ray v may be an SVector from the StaticArrays.jl package, and the Quaternion type which is taken from ReferenceFrameRotations.jl is already immutable. The default for the orbital ray is SVector. The v0 and s arrays contain immutable number types, and so these arrays MUST be mutable.

VectorField{V,Q}

Lie operator which acts on DAMaps e.g. dM/dt = FM where F is the VectorField and M is the DAMap. F can be promoted to a DAMap using exp(F).

Fields

• v::V – Orbital ray as truncated power series, expansion with scalar part equal to EXIT coordinates of map
• q::Q – Quaternion as truncated power series if spin is included, else nothing

Type Requirements

• V <: AbstractVector{<:TPS}
• Q <: Union{Quaternion{<:TPS},Nothing} and if Q != Nothing then eltype(Q) == eltype(V)

one(m::TaylorMap)

Creates an identity m with the same properties as m, including GTPSA Descriptor, spin, and stochasticity.

zero(m::TaylorMap)

Creates a zero m with the same properties as m including TPSA definiton, spin, and stochasticity.

checksymp(M)

Returns tranpose(M)*S*M - S, where S is the skew-symmetric matrix S = blkdiag([0 1; -1 0], ...). If M is symplectic, then the result should be a matrix containing all zeros. The non-symplectic parts of the matrix can be identified by those nonzero elements in the result.

exp!(
  m::DAMap, 
  F::VectorField, 
  m1::DAMap; 
  work_maps::NTuple{2, DAMap}=ntuple(t->zero(m1), Val{2}()), 
  work_q::Union{Quaternion,Nothing}=isnothing(m.q) ? nothing : zero(m.q)
)

Computes exp(F)*m1, and stores the result in m. Explicity, this is exp(F)*m1 = m1 + F*m1 + 1/2*F*(F*m1) + 1/6*F*(F*(F*m1)) + ..., where * is the operation of a VectorField on the map. See the documentation for mul! for more details of this operation. If the linear part of F is zero, then the number of iterations is equal to the number of higher orders left.

is_orbital_resonance(varidx, ords, nhv, res)

Checks if the monomial corresponds to a particular resonance (and resonance in the same family)

Note that for 3Q_x = integer, we also have 6Q_x though where

lambda(3Q_x) = integer where lambda*3 <= MAXORD + 1 according to notes are in the same family and so these must be left in too. These are the same family

each column of resonance corresponds to one res in the family

Each row is a multiple of previous

So For example for the 1Q_1 + 2Q2 + 3*Q3 = integer resonance: m = [ 1 2 3 ...; 2 4 6 ...; 3 6 9 ...];

nu_s + j dot Q = n

0Qx + 1Qy + 2*Qs = n

IMPORTANT:::: For spin resonances, there is only one resonance in each resonance family:

m = [0 ; 1 ]

ms = [1];

In the code, check - of m + ms

is_tune_shift(varidx, ords, nhv)

Checks if the monomial corresponds to a tune shift.

Input

• varidx – Current variable index (e.g. 1 is v, 2 is px, etc)
• ords – Array of monomial index as orders
• nhv – Number harmonic variables
• hamiltonian – Default is false, if the monomial is in a vector field and not a hamitlonian then this should be false.

lb!(G::VectorField, F::VectorField, H::VectorField; work_q::Union{Nothing,Quaternion}=nothing) -> G

Sets G equal to the Lie bracket of the vector fields F and H. Explicitly, including spin (with the lower case letter for the quaternion of the vector field), this is:

(G,g) = ⟨(F,f) , (H,h)⟩ = (F⋅∇H-H⋅∇F , [h,f]+F⋅∇h-G⋅∇f)

where [h,f] = h*f - f*h is just a quaternion commutator. See Equation 44.52 in the Bmad manual

locate_modes!(evecs, evals=nothing; sort=true, modes=nothing)

For each mode (every pair of rows in evecs), determines the eigenvector (each column in evecs) with the largest norm in that mode. The eigenvectors must have dimensionality of 2n, n ∈ ℤ and be next to each other pairs.

If sort is true, then the pairs of eigenvectors are sorted in-place. The eigenvalues will also be sorted if evals is provided. If sort is false, then the modes vector is filled with a mapping of the current eigenvector position to its mode (index of modes is the current eigenvector position, and the element at that index in modes is the mode that it corresponds to). For sort=true, modes == 1:n after calling this function.

If the mode locating is successful, true is returned. For unsuccessful mode locating, false is returned, modes will contain junk, and the eigenvectors/values may be partially sorted.

Assumes the eigenvectors (and eigenvalues if provided) are already in pairs (next to each other) starting at the first column.

Input:

• evecs – 2n x 2n matrix of eigenvectors already in pairs
• evals – Vector of length 2n containing the eigenvalues corresponding to the eigenvectors in evecs
• sort – (Optional) kwarg to specify whether or not to sort the eigenvectors, default is true

Output:

• ret – Returns true if the mode locating is successful, false if otherwise
• modes – (Optional) kwarg for eigvec -> mode mapping, length(evecs) == Int(size(evecs,1)/2)

log!(
  F::VectorField, 
  m1::DAMap; 
  work_maps::NTuple{3, DAMap}=ntuple(t->zero(m1), Val{3}()),
  work_vfs::NTuple{2, VectorField}=ntuple(t->zero(F), Val{2}()),
  work_q::Union{Quaternion,Nothing}=isnothing(m1.q) ? nothing : zero(m1.q)
)

Computes the log of the map m1 - that is, calculates the VectorField F that would represent the map m1 as a Lie exponent exp(F) - and stores the result in F. The map m1 should be close to the identity for this to converge quickly.

mat_eigen!(mat; sort=true, phase_modes=true)

Same as mat_eigen, but mutates mat for speed. See the documentation for mat_eigen for more details.

mat_eigen(mat; sort=true, phase_modes=true)

Given a matrix mat with an even number of rows/cols, calculates the eigenvectors and eigenvalues. For stable eigenmodes (complex conjugate eigenvector/eigenvalue pairs), the eigenvectors are normalized so that vⱼ'*S*vⱼ = +im for odd j, and -im for even j.

If sort is true, then each eigenvector/eigenvalue pair will be sorted according to the mode it best identifies with. A warning will be printed if the mode sorting fails. Mode sorting will automatically fail if more than 1 mode is unstable. Default is true.

If phase_modes is true, then each stable mode will be multiplied by a phase factor to make the first component of that mode (e.g. the x, y, or z component, NOT the px, py, or pz component) be fully real. This both makes the eigenvector "pretty", and in this case of weak coupling, ensures that vⱼ'*vⱼ₊₁ for j in (1, 3, 5) has a positive imaginary part so that the eigenvectors/eigenvalues for odd j are associated with the positive tune and even j the negative tune. This phase factor is harmless/useless to include in a highly-coupled matrix.

For complex matrices, Julia's eigen, which is called by mat_eigen, is type-unstable.

moveback_unstable!(F::Eigen) -> Int

Moves back eigenvectors with eigenvalues having a zero imaginary component to the end of the values and vectors arrays in the Eigen struct, and returns the number of unstable eigenvectors.

Note that if more than 1 mode is unstable, the pair of eigenvectors corresponding to a mode will not necessarily be next to each other at the end of the matrix.

mul!(m::DAMap, F::VectorField, m1::DAMap; work_q::Union{Quaternion,Nothing}=isnothing(m.q) ? nothing : zero(m.q)) -> m

Computes the Lie operator F acting on a DAMap m1, and stores the result in m. Explicity, that is F * m = (F.v, F.q) * (m.v, m.q) = (F.v ⋅ ∇ m.v , F.v ⋅ ∇ m.q + m.q*F.q)

normalize_eigenmode!(evec_pair, eval_pair, phase_mode::Integer=-1)

Normalizes the complex-conjugate eigenvector pair so that vⱼ'*S*vⱼ = +im for j=1, and -im for j=2. This may involve swapping the complex-conjugate eigenvectors/eigenvalues.

If phase_mode is set to a particular mode (1, 2, 3), then a phase will be multiplied by the eigenvectors to make them "pretty"; the phase will make the first component of that mode (e.g. the x, y, or z component) be fully real. This gives the eigenvectors a simple form in the weakly-coupled case, with the product v₁'*v₂ having a positive imaginary part, and ensures the odd numbered eigenvector/eigenvalues are associated with the tune and the even numbered ones have the negative of the tune. In cases with a lot of coupling, there is no simple way to define the positive or negative tune, and multiplying by this phase factor is harmless. See the "Tunes From One-Turn Matrix Eigen Analysis" section in the Bmad manual for more details.