Nostatic

Project.toml

@@ -1,7 +1,7 @@
 name = "NonlinearNormalForm"
 uuid = "05e19671-dec8-4f15-984f-54eaa6ca64be"
 authors = ["Matt Signorelli"]
-version = "0.5.1"
+version = "0.5.2"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

src/NonlinearNormalForm.jl

@@ -83,15 +83,18 @@ export        TaylorMap,
 # initially (before nv and np are known) is higher
 
 macro _DEFAULT_X0(NV)
-  return :(MVector{$(esc(NV))})
+  # return :(MVector{$(esc(NV))})
+  return :(Vector)
 end
 
 macro _DEFAULT_X(NN)
-  return :(MVector{$(esc(NN))})
+  # return :(MVector{$(esc(NN))})
+  return :(Vector)
 end
 
 macro _DEFAULT_S(NV)
-  return :(MMatrix{$(esc(NV)),$(esc(NV))})
+  # return :(MMatrix{$(esc(NV)),$(esc(NV))})
+  return :(Matrix)
 end
 
 const DEFAULT_NVARS = 6

src/normal.jl

@@ -593,12 +593,12 @@ end
 # The phase advance can be acquired from this map by atan(r11,r22), etc etc
 # Factors the map to: as ∘ a0 ∘ a1 ∘ a2 ∘ r . Spin canonise not yet implemented
 function factorise(
-  a::DAMap{V0};
+  a::DAMap;
   canonise::Integer=-1, # 0 = fully linear, 1 = linear w parameters, 2 = fully nonlinear
   phase=nothing, 
   damping::Bool=!isnothing(a.s),
   damp=nothing
-) where {V0<:StaticArray}
+)
   nn = ndiffs(a)
   nv = nvars(a)
   nhv = nhvars(a)
@@ -661,20 +661,11 @@ function factorise(
         error("Canonization including damping and coasting beam is not supported")
       end
       #a_matrix = jacobian(a, VARS)*jacobian(canoniser, VARS)
-      tmp = StaticArrays.sacollect(SMatrix{div(nhv, 2),div(nhv, 2)}, begin
-          a_matrix[2*i-1,2*j-1]*a_matrix[2*i,2*j]-a_matrix[2*i-1,2*j]*a_matrix[2*i,2*j-1]
-      end for j in 1:div(nhv, 2) for i in 1:div(nhv, 2))
+      tmp = block_det(a, a_matrix) 
       tmp2 = inv(tmp)
 
       dampt = sqrt.(sum(tmp2, dims=2))
-      Λi = StaticArrays.sacollect(SMatrix{nhv,nhv}, begin
-        if j == i
-          dampt[floor(Int, (i-1)/2)+1]
-        else
-          0
-        end
-      end for j in 1:nhv for i in 1:nhv)
-
+      Λi = double_diag(dampt)
       # now multiply canoniser matrix by this one
       a_rot = jacobian(canoniser, VARS)*Λi
       if !isnothing(damp)

src/utils.jl

@@ -512,4 +512,44 @@ function checksymp(m::TaylorMap)
   end
   return out-S
 end
-# =================================================================================== #
+# =================================================================================== #
+# Block determinants and double diag (for damping)
+
+function block_det(m::DAMap{V0}, a_matrix::SMatrix) where {V0<:StaticArray}
+  nhv = nhvars(m)
+  return StaticArrays.sacollect(SMatrix{div(nhv, 2),div(nhv, 2)}, begin
+    a_matrix[2*i-1,2*j-1]*a_matrix[2*i,2*j]-a_matrix[2*i-1,2*j]*a_matrix[2*i,2*j-1]
+  end for j in 1:div(nhv, 2) for i in 1:div(nhv, 2))
+end
+
+function block_det(m::DAMap, a_matrix)
+  nhv = nhvars(m)
+  nd = div(nhv, 2)
+  out = zeros(eltype(a_matrix), nd, nd)
+  for i in 1:nd
+    for j in 1:nd
+      out[i,j] = det(view(a_matrix, (2*i-1):(2*i), (2*j-1):(2*j)))
+    end
+  end
+  return out
+end
+
+function double_diag(v::V) where {V<:StaticArray}
+  n = length(V)
+  return StaticArrays.sacollect(SMatrix{2*n,2*n}, begin
+    if j == i
+      v[floor(Int, (i-1)/2)+1]
+    else
+      0
+    end
+  end for j in 1:(2*n) for i in 1:(2*n))
+end
+
+function double_diag(v)
+  n = length(v)
+  out = zeros(eltype(v), 2*n)
+  for i in 1:(2*n)
+    out[i] = v[floor(Int, (i-1)/2)+1]
+  end
+  return Diagonal(out)
+end