working on spin normal form

fpp-ptc-sandbox/code/CMakeLists.txt

@@ -32,6 +32,7 @@ set (EXE_SPECS
   cmake_files/cmake.z_radiation_matt_fake_maps
   cmake_files/cmake.z_canonize_matt_fake_maps
   cmake_files/cmake.z_resonance
+  cmake_files/cmake.z_spin1
 )
 
 

fpp-ptc-sandbox/code/cmake_files/cmake.z_spin1

@@ -0,0 +1,14 @@
+set (EXENAME z_spin1)
+
+FILE (GLOB SRC_FILES z_spin1.f90)
+
+set (INC_DIRS
+)
+
+set (LINK_LIBS
+  forest
+  bsim
+  bmad
+  sim_utils
+  ${ACC_BMAD_LINK_LIBS}
+)

fpp-ptc-sandbox/code/z_spin1.f90

@@ -0,0 +1,160 @@
+
+Conversation opened. 1 unread message.
+
+Skip to content
+Using Cornell Mail with screen readers
+1 of 11,666
+spin
+External
+Inbox
+Patrice Nishikawa
+
+Attachments11:11 AM (0 minutes ago)
+
+to me
+
+ One attachment  •  Scanned by Gmail
+
+Compose:
+Quaternion
+	MinimizePop-outClose
+
+           4  Dimensional TPSA/DA map 
+
+          1, NO =    3, NV =    4, INA =  206
+ *********************************************
+
+    I  COEFFICIENT          ORDER   EXPONENTS
+      NO =     3      NV =     4
+   1  -1.076999868651706       0.000000000000000       1  0  0  0
+   1  0.5007864172415770       0.000000000000000       0  1  0  0
+    -2   0.000000000000000       0.000000000000000       0  0  0  0
+
+          1, NO =    3, NV =    4, INA =  205
+ *********************************************
+
+    I  COEFFICIENT          ORDER   EXPONENTS
+      NO =     3      NV =     4
+   1 -0.1099413460130740       0.000000000000000       1  0  0  0
+   1 -0.8773843848341096       0.000000000000000       0  1  0  0
+    -2   0.000000000000000       0.000000000000000       0  0  0  0
+
+          1, NO =    3, NV =    4, INA =  204
+ *********************************************
+
+    I  COEFFICIENT          ORDER   EXPONENTS
+      NO =     3      NV =     4
+   1  0.8286937949539068       0.000000000000000       0  0  1  0
+   1  0.9383342212004021       0.000000000000000       0  0  0  1
+    -2   0.000000000000000       0.000000000000000       0  0  0  0
+
+          1, NO =    3, NV =    4, INA =  203
+ *********************************************
+
+    I  COEFFICIENT          ORDER   EXPONENTS
+      NO =     3      NV =     4
+   1 -0.6705862447736176       0.000000000000000       0  0  1  0
+   1  0.4474101055423855       0.000000000000000       0  0  0  1
+    -2   0.000000000000000       0.000000000000000       0  0  0  0
+  Spin Matrix is identity 
+  Quaternion  
+  c_quaternion 
+
+          1, NO =    3, NV =    4, INA =  186
+ *********************************************
+
+    I  COEFFICIENT          ORDER   EXPONENTS
+      NO =     3      NV =     4
+   0  0.9050015613244589       0.000000000000000       0  0  0  0
+   1 -0.9503967950294451E-01   0.000000000000000       1  0  0  0
+   1  -1.112362312442405       0.000000000000000       0  1  0  0
+   2  -1.543019756088252       0.000000000000000       2  0  0  0
+   2  0.6629121268802496       0.000000000000000       1  1  0  0
+   2  -4.761727132588617       0.000000000000000       0  2  0  0
+   2  -10.84057903687762       0.000000000000000       0  0  2  0
+   2  -7.625306718426536       0.000000000000000       0  0  1  1
+   2  -10.12597288417461       0.000000000000000       0  0  0  2
+   3 -0.7035748778319140       0.000000000000000       3  0  0  0
+   3  -8.022246878574116       0.000000000000000       2  1  0  0
+   3   6.226975993860259       0.000000000000000       1  2  0  0
+   3  -8.261773184839884       0.000000000000000       0  3  0  0
+   3   4.891464788238377       0.000000000000000       1  0  2  0
+   3  -73.30062999530527       0.000000000000000       0  1  2  0
+   3  -31.68508851608486       0.000000000000000       1  0  1  1
+   3  -47.51630813007174       0.000000000000000       0  1  1  1
+   3  -16.25997300096116       0.000000000000000       1  0  0  2
+   3  -53.61950386253088       0.000000000000000       0  1  0  2
+   -19   0.000000000000000       0.000000000000000       0  0  0  0
+
+          1, NO =    3, NV =    4, INA =  185
+ *********************************************
+
+    I  COEFFICIENT          ORDER   EXPONENTS
+      NO =     3      NV =     4
+   1  -1.640737988472106       0.000000000000000       0  0  1  0
+   1 -0.5580901068366162       0.000000000000000       0  0  0  1
+   2   4.495757658080316       0.000000000000000       1  0  1  0
+   2  -10.21919374003079       0.000000000000000       0  1  1  0
+   2  -15.73046762538523       0.000000000000000       1  0  0  1
+   2  -1.409846817974418       0.000000000000000       0  1  0  1
+   3   6.527589613109422       0.000000000000000       2  0  1  0
+   3   11.47144226575744       0.000000000000000       1  1  1  0
+   3  -24.60780806613684       0.000000000000000       0  2  1  0
+   3   8.459672725402381       0.000000000000000       2  0  0  1
+   3  -30.76909316739637       0.000000000000000       1  1  0  1
+   3  -7.445167915902111       0.000000000000000       0  2  0  1
+   3   18.88043769734887       0.000000000000000       0  0  3  0
+   3   17.12656111216577       0.000000000000000       0  0  2  1
+   3   28.81440459189866       0.000000000000000       0  0  1  2
+   3   10.15488346067025       0.000000000000000       0  0  0  3
+   -16   0.000000000000000       0.000000000000000       0  0  0  0
+
+          1, NO =    3, NV =    4, INA =  184
+ *********************************************
+
+    I  COEFFICIENT          ORDER   EXPONENTS
+      NO =     3      NV =     4
+   0 -0.4254082439261045       0.000000000000000       0  0  0  0
+   1 -0.2021847473949788       0.000000000000000       1  0  0  0
+   1  -2.366408370999338       0.000000000000000       0  1  0  0
+   2  -3.223914208566652       0.000000000000000       2  0  0  0
+   2   2.783459801598894       0.000000000000000       1  1  0  0
+   2  -2.093873957716453       0.000000000000000       0  2  0  0
+   2  -19.86617375271195       0.000000000000000       0  0  2  0
+   2  -13.85555150150932       0.000000000000000       0  0  1  1
+   2  -20.81529816009441       0.000000000000000       0  0  0  2
+   3  0.3801948418494159       0.000000000000000       3  0  0  0
+   3   3.430992680594370       0.000000000000000       2  1  0  0
+   3  -1.910813103176785       0.000000000000000       1  2  0  0
+   3   6.522696900714276       0.000000000000000       0  3  0  0
+   3  0.3629422487711790       0.000000000000000       1  0  2  0
+   3   22.88711306753415       0.000000000000000       0  1  2  0
+   3  -23.02332398486435       0.000000000000000       1  0  1  1
+   3   21.98983984271880       0.000000000000000       0  1  1  1
+   3  -12.87324854557045       0.000000000000000       1  0  0  2
+   3   48.06116540649838       0.000000000000000       0  1  0  2
+   -19   0.000000000000000       0.000000000000000       0  0  0  0
+
+          1, NO =    3, NV =    4, INA =  183
+ *********************************************
+
+    I  COEFFICIENT          ORDER   EXPONENTS
+      NO =     3      NV =     4
+   1 -0.1642941931757383       0.000000000000000       0  0  1  0
+   1 -0.5536958539643637       0.000000000000000       0  0  0  1
+   2   11.82602203326790       0.000000000000000       1  0  1  0
+   2  -1.438884425557026       0.000000000000000       0  1  1  0
+   2   8.508296470635484       0.000000000000000       1  0  0  1
+   2  -13.84023855236655       0.000000000000000       0  1  0  1
+   3  -16.07892195940298       0.000000000000000       2  0  1  0
+   3   82.94352763492097       0.000000000000000       1  1  1  0
+   3  -10.27052716192868       0.000000000000000       0  2  1  0
+   3   1.559247096195197       0.000000000000000       2  0  0  1
+   3   56.55124485074391       0.000000000000000       1  1  0  1
+   3  -58.36553225410486       0.000000000000000       0  2  0  1
+   3   1.423571285444478       0.000000000000000       0  0  3  0
+   3  -72.24414550183792       0.000000000000000       0  0  2  1
+   3  -62.81269926655657       0.000000000000000       0  0  1  2
+   3  -60.19180222963241       0.000000000000000       0  0  0  3
+   -16   0.000000000000000       0.000000000000000       0  0  0  0
+ No Stochastic Radiation 

src/map/compose.jl

@@ -3,12 +3,12 @@
 # --- compose! ---
 
 """
-compose!(m::DAMap, m2::DAMap, m1::DAMap; keep_scalar::Bool=true, work_ref::Union{Nothing,Vector{<:Union{Float64,ComplexF64}}}=nothing, dospin::Bool=true, work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS}}}}=prep_comp_work_prom(m,m2,m1))
+compose!(m::DAMap, m2::DAMap, m1::DAMap; keep_scalar::Bool=true, work_ref::Union{Nothing,Vector{<:Union{Float64,ComplexF64}}}=nothing, dospin::Bool=true, work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS64}}}}=prep_comp_work_prom(m,m2,m1))
 
 In-place `DAMap` composition which calculates `m = m2 ∘ m1`, ignoring the scalar part of `m1`.
 
 Aliasing `m === m2` is allowed, but NOT `m === m1`. `m2 === m1` is fine. The destination map `m` should be 
-properly set up (with correct types promoted if necessary), and have `m.x[1:nv]` (and `m.Q.q` if spin) 
+properly set up (with correct types promoted if necessary), and have `m.x[1:nv]` (and `m.Q` if spin) 
 containing allocated TPSs.
 
 If `keep_scalar` is set to `true`, then the scalar part of `m1` is retained. In this case, a temporary 
@@ -24,9 +24,9 @@ See the documentation for `compose_it!` for information on `work_prom`.
 - `keep_scalar` -- Specify whether to keep the scalar part in `m1` or throw it away. If `true`, a temporary vector storing the scalar part must be used. Default is true
 - `work_ref` -- If `keep_scalar` is true, the temporary vector can be provided in this keyword argument. If `nothing` is provided, the temporary will be created internally. Default is `nothing`
 - `dospin` -- Specify whether or not to include the quaternions in the concatenation. Default is `true`
-- `work_prom` -- Temporary vector of allocated `ComplexTPS`s when there is implicit promotion. See the `compose_it!` documentation for more details. Default is output from `prep_comp_work_prom(m, m2, m1)`
+- `work_prom` -- Temporary vector of allocated `ComplexTPS64`s when there is implicit promotion. See the `compose_it!` documentation for more details. Default is output from `prep_comp_work_prom(m, m2, m1)`
 """
-function compose!(m::DAMap, m2::DAMap, m1::DAMap; keep_scalar::Bool=true, work_ref::Union{Nothing,Vector{<:Union{Float64,ComplexF64}}}=nothing, dospin::Bool=true, work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS}}}}=prep_comp_work_prom(m,m2,m1))
+function compose!(m::DAMap, m2::DAMap, m1::DAMap; keep_scalar::Bool=true, work_ref::Union{Nothing,Vector{<:Union{Float64,ComplexF64}}}=nothing, dospin::Bool=true, work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS64}}}}=prep_comp_work_prom(m,m2,m1))
   checkinplace(m, m2, m1)
 
   # DAMap setup:
@@ -72,12 +72,12 @@ function compose!(m::DAMap, m2::DAMap, m1::DAMap; keep_scalar::Bool=true, work_r
 end
 
 """
-    compose!(m::TPSAMap, m2::TPSAMap, m1::TPSAMap; dospin::Bool=true, work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS}}}}=prep_comp_work_prom(m,m2,m1))
+    compose!(m::TPSAMap, m2::TPSAMap, m1::TPSAMap; dospin::Bool=true, work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS64}}}}=prep_comp_work_prom(m,m2,m1))
 
 In-place `TPSAMap` composition which calculates `m = m2 ∘ m1`, including the scalar part of `m1`.
 
 Aliasing `m === m2` is allowed, but NOT `m === m1`. `m2 === m1` is fine. The destination map `m` should be 
-properly set up (with correct types promoted if necessary), and have `m.x[1:nv]` (and `m.Q.q` if spin) 
+properly set up (with correct types promoted if necessary), and have `m.x[1:nv]` (and `m.Q` if spin) 
 containing allocated TPSs.
 
 If `dospin` is `true`, then the quaternion part of the maps will be composed as well. Default is `true`.
@@ -86,9 +86,9 @@ See the documentation for `compose_it!` for information on `work_prom`.
 
 ### Keyword Arguments
 - `dospin` -- Specify whether or not to include the quaternions in the concatenation. Default is `true`
-- `work_prom` -- Temporary vector of allocated `ComplexTPS`s when there is implicit promotion. See the `compose_it!` documentation for more details. Default is output from `prep_comp_work_prom(m, m2, m1)`
+- `work_prom` -- Temporary vector of allocated `ComplexTPS64`s when there is implicit promotion. See the `compose_it!` documentation for more details. Default is output from `prep_comp_work_prom(m, m2, m1)`
 """
-function compose!(m::TPSAMap, m2::TPSAMap, m1::TPSAMap; dospin::Bool=true, work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS}}}}=prep_comp_work_prom(m,m2,m1))
+function compose!(m::TPSAMap, m2::TPSAMap, m1::TPSAMap; dospin::Bool=true, work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS64}}}}=prep_comp_work_prom(m,m2,m1))
   checkinplace(m, m2, m1)
 
   # TPSAMap setup:

src/map/compose_it.jl

@@ -9,23 +9,23 @@ Assumes the destination map is properly set up (with correct types promoted if n
 that `m.x[1:nv]` (and `m.Q` if spin) contain allocated TPSs. The parameters part of `m.x` (`m.x[nv+1:nn]`) 
 does not need to contain allocated TPSs.
 
-If promotion is occuring, then one of the maps is `ComplexTPS` and the other `TPS`, with output map `ComplexTPS`. 
+If promotion is occuring, then one of the maps is `ComplexTPS64` and the other `TPS`, with output map `ComplexTPS64`. 
 Note that the spin part is required to agree with the orbital part in terms of type by definition of the `TaylorMap` 
-struct. `work_prom` can optionally be passed as a tuple containing the temporary `ComplexTPS`s if promotion is occuring:
+struct. `work_prom` can optionally be passed as a tuple containing the temporary `ComplexTPS64`s if promotion is occuring:
 
 If `eltype(m.x) != eltype(m1.x)` (then `m1` must be promoted):
-`work_prom[1] = m1x_prom  # Length >= nv+np, Vector{ComplexTPS}`
+`work_prom[1] = m1x_prom  # Length >= nv+np, Vector{ComplexTPS64}`
 
 else if `eltype(m.x) != eltype(m2.x)` (then `m2` must be promoted):
 ```
-work_prom[1] = m2x_prom  # Length >= nv, Vector{ComplexTPS}
-work_prom[2] = m2Q_prom  # Length >= 4, Vector{ComplexTPS}
+work_prom[1] = m2x_prom  # Length >= nv, Vector{ComplexTPS64}
+work_prom[2] = m2Q_prom  # Length >= 4, Vector{ComplexTPS64}
 ```
-Note that the `ComplexTPS`s in the vector(s) must be allocated and have the same `Descriptor`.
+Note that the `ComplexTPS64`s in the vector(s) must be allocated and have the same `Descriptor`.
 
 If spin is included, not that the final quaternion concatenation step mul! will create allocations
 """ 
-function compose_it!(m::$t, m2::$t, m1::$t; dospin::Bool=true, dostochastic::Bool=true, work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS}}}}=prep_comp_work_prom(m,m2,m1))
+function compose_it!(m::$t, m2::$t, m1::$t; dospin::Bool=true, dostochastic::Bool=true, work_prom::Union{Nothing,Tuple{Vararg{Vector{<:ComplexTPS64}}}}=prep_comp_work_prom(m,m2,m1))
   checkinplace(m, m2, m1)
   @assert !(m === m1) "Cannot compose_it!(m, m2, m1) with m === m1"
 
@@ -69,9 +69,9 @@ function compose_it!(m::$t, m2::$t, m1::$t; dospin::Bool=true, dostochastic::Boo
 
   # Spin:
   if !isnothing(m.Q) && dospin
-    compose!(m.Q.q, m2.Q.q, m1.x, work=Q_prom)
-    mul!(m.Q, m1.Q, m.Q)
-  end
+    compose!(m.Q, m2.Q, m1.x, work=Q_prom)
+    mul!(m.Q, m.Q, m1.Q) # This will be made faster eventually
+  end 
 
   # Stochastic
   # MAKE THIS FASTER!

src/map/ctors.jl

@@ -26,10 +26,10 @@ function $t(m::Union{TaylorMap{S,T,U,V,W},Probe{S,T,U,V,W}}; use::UseType=m, idp
 
   # set quaternion
   if !isnothing(outm.Q)
-    setTPS!(outm.Q.q[1], m.Q.q[1], change=true)
-    setTPS!(outm.Q.q[2], m.Q.q[2], change=true)
-    setTPS!(outm.Q.q[3], m.Q.q[3], change=true)
-    setTPS!(outm.Q.q[4], m.Q.q[4], change=true)
+    setTPS!(outm.Q.q0, m.Q.q0, change=true)
+    setTPS!(outm.Q.q1, m.Q.q1, change=true)
+    setTPS!(outm.Q.q2, m.Q.q2, change=true)
+    setTPS!(outm.Q.q3, m.Q.q3, change=true)
   end
 
   # set the reference orbit properly
@@ -57,7 +57,7 @@ end
 """
     $($t)(;use::UseType=GTPSA.desc_current, x::Vector=vars(getdesc(use)), x0::Vector=zeros(eltype(eltype(x)), numvars(use)), Q::Union{Quaternion,Nothing}=nothing, E::Union{Matrix,Nothing}=nothing, idpt::Union{Bool,Nothing}=nothing, spin::Union{Bool,Nothing}=nothing, FD::Union{Bool,Nothing}=nothing) 
 
-Constructs a $($t) with the passed vector of `TPS`/`ComplexTPS` as the orbital ray, and optionally the entrance 
+Constructs a $($t) with the passed vector of `TPS`/`ComplexTPS64` as the orbital ray, and optionally the entrance 
 coordinates `x0`, `Quaternion` for spin `Q`, and FD matrix `E` as keyword arguments. The helper keyword 
 arguments `spin` and `FD` may be set to `true` to construct a $($t) with an identity quaternion/FD 
 matrix, or `false` for no spin/FD. Note that setting `spin`/`FD` to any `Bool` value without `Q` or `E` 
@@ -129,12 +129,12 @@ function $t(;use::UseType=GTPSA.desc_current, x::Vector=vars(getdesc(use)), x0::
   # set quaternion
   if !isnothing(outm.Q)
     if !isnothing(Q)
-      setTPS!(outm.Q.q[1], Q.q[1], change=true)
-      setTPS!(outm.Q.q[2], Q.q[2], change=true)
-      setTPS!(outm.Q.q[3], Q.q[3], change=true)
-      setTPS!(outm.Q.q[4], Q.q[4], change=true)
+      setTPS!(outm.Q.q0, Q.q0, change=true)
+      setTPS!(outm.Q.q1, Q.q1, change=true)
+      setTPS!(outm.Q.q2, Q.q2, change=true)
+      setTPS!(outm.Q.q3, Q.q3, change=true)
     else
-      outm.Q.q[1][0] = 1
+      outm.Q.q0[0] = 1
     end
   end
 
@@ -168,7 +168,7 @@ function $t(M::AbstractMatrix; use::UseType=GTPSA.desc_current, x0::Vector=zeros
   nv >= size(M,1) || error("Number of rows in matrix > number of variables in GTPSA!")
 
   if eltype(M) <: Complex
-    T = ComplexTPS
+    T = ComplexTPS64
   else
     T = TPS{Float64}
   end
@@ -245,10 +245,10 @@ end
 function zero_op(m2::Union{$t,Number}, m1::Union{$t,Number})
   outtype = promote_type(typeof(m1),typeof(m2))
 
-  # If either inputs are ComplexTPS, use those parameters
+  # If either inputs are ComplexTPS64, use those parameters
   if m2 isa $t
     if m1 isa $t
-      if eltype(m2.x) == ComplexTPS
+      if eltype(m2.x) == ComplexTPS64
         return zero(outtype,use=m2,idpt=m2.idpt)
       else
         return zero(outtype,use=m1,idpt=m1.idpt)
@@ -282,7 +282,7 @@ function one(t::Type{$t{S,T,U,V,W}}; use::UseType=GTPSA.desc_current, idpt::W=no
   end
 
   if !isnothing(m.Q)
-    @inbounds m.Q.q[1][0] = 1
+    @inbounds m.Q.q0[0] = 1
   end
 
   return m

src/map/inv.jl

@@ -76,7 +76,7 @@ function inv!(m::TaylorMap, m1::TaylorMap; dospin::Bool=true, work_ref::Union{No
   # Now do quaternion: inverse of q(z0) is q^-1(M^-1(zf))
   if !isnothing(m.Q) && dospin
     inv!(m.Q, m1.Q)
-    compose!(m.Q.q, m.Q.q, m.x)
+    compose!(m.Q, m.Q, m.x)
   end
 
   if m1 === m