Skip to content

Commit

Permalink
Fix spin tune degeneracy problem. (#1059)
Browse files Browse the repository at this point in the history 
* Fix spin tune degeneracy problem.
  • Loading branch information
@DavidSagan
DavidSagan authored Jul 14, 2024
1 parent c8d13b3 commit dc5d0a5
Show file tree
Hide file tree
Showing 13 changed files with 186 additions and 69 deletions.
2 changes: 1 addition & 1 deletion bmad/doc/cover-page.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@

\begin{flushright}
\large
Revision: July 10, 2024 \\
Revision: July 14, 2024 \\
\end{flushright}

\pdfbookmark[0]{Preamble}{Preamble}
Expand Down
18 changes: 9 additions & 9 deletions bmad/doc/hdf5-programming.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,15 +52,15 @@ \section{HDF5 Particle Beam Data Storage}
\begin{center}
\tt
\begin{tabular}{lll} \toprule
{\em Beam Physics Parameter} & {\em Bmad Equivalent} & {\em Notes} \\ \midrule
time & -\%vec(5) / (c \%beta) & time - ref_time. See \Eq{zbctt} \\
timeOffset & \%t - time & reference time \\
totalMomentumOffset & \%p0c & \\
sPosition & \%s & See Fig.~\ref{f:local.coords} \\
weight & \%charge & Macro bunch charge \\
branchIndex & \%ix_branch & \\
elementIndex & \%ix_ele & \\
locationInElement & \%location & See below \\
{\em Beam Physics Parameter} & {\em Bmad Equivalent} & {\em Notes} \\ \midrule
time & -\%vec(5) / (c \%beta) & time - ref_time. See \Eq{zbctt} \\
timeOffset & \%t - time (beam physics) & reference time \\
totalMomentumOffset & \%p0c & \\
sPosition & \%s & See Fig.~\ref{f:local.coords} \\
weight & \%charge & Macro bunch charge \\
branchIndex & \%ix_branch & \\
elementIndex & \%ix_ele & \\
locationInElement & \%location & See below \\
particleStatus & \%state & See the \%state table in \sref{s:coord.struct} \\ \bottomrule
\end{tabular}
\end{center}
Expand Down
26 changes: 26 additions & 0 deletions bmad/doc/spin.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -317,6 +317,7 @@ \section{Polarization Limits and Polarization/Depolarization Rates}
P_{bks} = \pm \frac{8}{5 \, \sqrt{3}}
\frac{\displaystyle \oint ds \, g^3 \, \what\bfb \dotproduct \bfn_0}
{\displaystyle \oint ds \, g^3 \left( 1 - \frac{2}{9} (\bfn_0 \dotproduct \what\bfs)^2 \right)}
\label{pbks}
\end{equation}
where $g = 1/\rho$ is the bending strength ($\rho$ is the bending radius), $\what\bfs$ is the unit
vector in the direction of motion, and $\what\bfb$ is defined to be
Expand Down Expand Up @@ -414,6 +415,31 @@ \section{Polarization Limits and Polarization/Depolarization Rates}
where $\tau_{dep}^{-1}$ is calculated via particle tacking and the other quantities are calculated by
integrals over the closed orbit ignoring beam size effects.

%-----------------------------------------------------------------
\section{Bmad Tune and $n_0$ convention}
\label{s:spin.tune}

As mentioned above, at any point in a ring the direction of the closed orbit invariant spin $\bfn_0$
is ambiguous since, if $\bfn_0$ is a valid invariant spin direction, then so is $-\bfn_0$. The same
is true of the spin tune $\nu_0$ and if $\nu_0$ is the spin tune associated with $\bfn_0$, $-\nu_0$
is the spin tune associated with $-\bfn_0$. This ambiguity complicates various calculations. For
example, to do the integral in \Eq{pbks}, it is necessary to make sure that the direction of
$\bfn_0(s_1)$, the invariant spin at $s_1$, must be consistent with the direction of the invariant
spin at $s_2$, $\bfn_0(s_2)$. That is, the two must be related via
\begin{equation}
\bfn_0(s_2) = \bfq_{21} \, \bfn_0(s_1) \, \BAR\bfq_{21}
\label{n0qn0q}
\end{equation}
where $\bfq_{21}$ is the closed orbit spin transport quaternion from $s_1$ to $s_2$. The tune must
also be computed consistent with the choice of invariant spin direction. This is important when
calculating sum and difference resonance strengths (\Eqs{x12pq1} and \eq{x12pq2}).

To ensure a consistent invariant spin direction, \Eq{n0qn0q} is used when calculating integrals
involving $\bfn_0$. To ensure a consistent tune, \bmad uses the following convention that the spin
tune will always be in the range $[0, \pi]$, and the direction of $\bfn_0$ will be chosen to be
consistent with this choice in tune (it is left as an exercise for the reader to prove that there
will always be exactly one spin tune in the range $[0, \pi]$.

%-----------------------------------------------------------------
\section{Linear \texorpdfstring{$\partial\bfn/\partial\delta$}{dn/dpz} Calculation}
\label{s:dn.calc}
Expand Down
2 changes: 1 addition & 1 deletion bmad/output/write_lattice_in_julia.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,7 @@ subroutine write_lattice_in_julia(bmad_file, lat, julia_file)

if (ie == 0) then
line = trim(line) // ', pc_ref = ' // re_str(ele%value(p0c$))
line = trim(line) // ', species_ref = ' // trim(species_name(ele%ref_species))
line = trim(line) // ', species_ref = ' // trim(openpmd_species_name(ele%ref_species))
if (ele%a%beta /= 0) line = trim(line) // ', twiss.a.beta = ' // re_str(ele%a%beta)
if (ele%b%beta /= 0) line = trim(line) // ', twiss.b.beta = ' // re_str(ele%b%beta)
if (ele%a%alpha /= 0) line = trim(line) // ', twiss.a.alpha = ' // re_str(ele%a%alpha)
Expand Down
13 changes: 10 additions & 3 deletions bmad/spin/spin_quat_resonance_strengths.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,15 +22,22 @@ subroutine spin_quat_resonance_strengths (orb_evec, spin_q, xi_sum, xi_diff)

implicit none

real(rp) spin_q(0:3,0:6), xi_sum, xi_diff, nn0(3)
real(rp) spin_q(0:3,0:6), xi_sum, xi_diff, nn0(3), anorm
complex(rp) orb_evec(6), qv(0:3), qv2(0:3), np(0:3), nm(0:3)

integer k

!

nn0 = spin_q(1:3,0)
nn0 = nn0 / norm2(nn0)
nn0 = spin_q(1:3,0)
anorm = norm2(nn0)
if (anorm == 0) then
xi_sum = 0
xi_diff = 0
return
endif

nn0 = sign_of(spin_q(0,0)) * nn0 / norm2(nn0)

np(0) = 0.5_rp
np(1:3) = i_imag * nn0 / 2
Expand Down
1 change: 1 addition & 0 deletions regression_tests/cmake_files/cmake.spin_general_test
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ set (INC_DIRS
)

set (LINK_LIBS
tao
bmad
sim_utils
${ACC_BMAD_LINK_LIBS}
Expand Down
35 changes: 34 additions & 1 deletion regression_tests/spin_general_test/output.correct
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,40 @@
"Spin-Res-2" ABS 1E-9 1.943274E-01 4.581804E-01 1.462823E-01
"Spin-Res-3" ABS 1E-9 -1.742571E-02 1.404553E-02 1.579951E-02
"dPTC-Quad" ABS 0 0.000049048 0.000000000 -0.000029432
"dBmad-Quad" ABS 0 0.000049048 0.000000000 -0.000029432
"dBmad-Quad" ABS 0 0.000049048 -0.000000000 -0.000029432
"dRot" ABS 1e-10 2.93E-01 4.13E-01 4.94E-01
"Taylor-Taylor" ABS 1e-10 -0.50000000 0.40000000 0.30000000
"PTC-Taylor" ABS 1e-10 -0.50000000 0.40000000 0.30000000
"Polarization Limit ST" ABS 1e-8 0.00000000
"Polarization Limit DK" ABS 1e-8 0.00000000
"Polarization Limits DK (a,b,c-modes)" ABS 1e-8 0.00000000 0.00000000 0.00000000
"Polarization Limits DK (bc,ac,ab-modes)" ABS 1e-8 0.00000000 0.00000000 0.00000000
"Polarization Rate BKS" REL 1e-8 1.44802509E-05
"Depolarization Rate" REL 1e-8 5.69044334E-04
"Depolarization Rate Partial" REL 1e-8 5.10277794E-04 4.10236457E-06 1.04048639E-06
"Depolarization Rate Partial2" REL 1e-8 1.81658388E-06 5.06163317E-04 5.76485078E-04
"Integral g^3 * b_hat * n_0" REL 1e-8 0.00000000E+00
"Integral g^3 * b_hat * dn/ddelta" REL 1e-8 2.71191375E-18
"Integral g^3 (1 - 2(n * s_hat)/9)" REL 1e-8 1.34575854E-02
"Integral g^3 * 11 (dn/ddelta)^2 / 9" REL 1e-8 5.28855663E-01
"orb-evec-1" ABS 1e-8 2.406797 0.399089 0.703623 0.002991
"spin-re-evec-1" REL 1e-8 -4.68985050E+00 4.22193710E-01
"spin-im-evec-1" REL 1e-8 4.03529761E-01 4.90355786E+00
"orb-evec-2" ABS 1e-8 2.406797 0.399089 0.703623 0.002991
"spin-re-evec-2" REL 1e-8 -4.68985050E+00 4.22193710E-01
"spin-im-evec-2" REL 1e-8 -4.03529761E-01 -4.90355786E+00
"orb-evec-3" ABS 1e-8 0.487345 0.059887 0.251203 0.001242
"spin-re-evec-3" REL 1e-8 1.21103352E+00 -7.14952982E-01
"spin-im-evec-3" REL 1e-8 1.75457148E-02 -2.99685682E-02
"orb-evec-4" ABS 1e-8 0.487345 0.059887 0.251203 0.001242
"spin-re-evec-4" REL 1e-8 1.21103352E+00 -7.14952982E-01
"spin-im-evec-4" REL 1e-8 -1.75457148E-02 2.99685682E-02
"orb-evec-5" ABS 1e-8 0.559385 0.081892 3.237229 0.155199
"spin-re-evec-5" REL 1e-8 2.97259086E-04 -2.79231130E-03
"spin-im-evec-5" REL 1e-8 -2.16290507E-02 5.17073042E-02
"orb-evec-6" ABS 1e-8 0.559385 0.081892 3.237229 0.155199
"spin-re-evec-6" REL 1e-8 2.97259086E-04 -2.79231130E-03
"spin-im-evec-6" REL 1e-8 2.16290507E-02 -5.17073042E-02
"Res-Strength-1" REL 1e-8 0.0665726 0.0392010
"Res-Strength-2" REL 1e-8 0.4413666 0.1248851
"Res-Strength-3" REL 1e-8 0.0130688 0.0146930
20 changes: 20 additions & 0 deletions regression_tests/spin_general_test/small_ring.bmad
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
parameter[p0c] = 1e9
bmad_com[radiation_damping_on] = f
bmad_com[radiation_fluctuations_on] = f
bmad_com[spin_tracking_on] = T

q1: quad, l = 0.5, k1 = 0.6, tilt = 0.05
q2: quad, l = 0.5, k1 = -0.5

d: drift, l = 0.3
b: sbend, l = 3, angle = pi/8, e1 = 0.07, e2 = 0.07

s1: sextupole, l = 0.01
s2: sextupole, l = 0.01

rf: rfcavity, l = 1.2, rf_frequency = 250e6, voltage = 1e6

sector: line = (s1, d, q1, d, b, d, q2, d, s2, rf, b)
ring: line = (1*sector)

use, ring
60 changes: 57 additions & 3 deletions regression_tests/spin_general_test/spin_general_test.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,23 +2,29 @@ program spin_general_test

use bmad
use pointer_lattice, only: c_linear_map, operator(*), assignment(=)
use tao_interface

implicit none

type (lat_struct), target :: lat
type (branch_struct), pointer :: branch
type (ele_struct), pointer :: ele
type (ele_struct) t_ele
type (coord_struct) orb0, orb_start, orb_end, orb1, orb2
type (spin_orbit_map1_struct) map1, inv_map1
type (tao_lattice_branch_struct), target :: tao_branch
type (tao_spin_map_struct) :: sm

real(rp) spin_a(3), spin_b(3), spin0(3), dr(6), a_quat(0:3), n_vec(3)
real(rp) mat6(6,6), smap(0:3,0:6), n0(3), q0(0:3), q1(0:3), t, q, xi_sum, xi_diff
real(rp) mat6(6,6), smap(0:3,0:6), q0(0:3), q1(0:3), t, q, vec(6)
real(rp) xi_sum, xi_diff, n0(3)

complex(rp) orb_eval(6), orb_evec(6,6), spin_evec(6,3)
integer i, j, nargs
integer i, j, nargs, n_eigen
logical print_extra, err, err_flag

character(40) :: lat_file = 'spin_general_test.bmad'
logical print_extra, err, err_flag
character(100) excite_zero(3)

namelist / param / dr

Expand Down Expand Up @@ -186,6 +192,54 @@ program spin_general_test
call track1 (orb_start, lat%ele(2), lat%param, orb_end)
write (1, '(a, 3f12.8)') '"PTC-Taylor" ABS 1e-10 ', orb_end%spin

!---------------------------------

call bmad_parser('small_ring.bmad', lat)
branch => lat%branch(0)
call closed_orbit_calc(lat, tao_branch%orbit, 6)

excite_zero = ''
call tao_spin_polarization_calc (branch, tao_branch, excite_zero, '')
call spin_concat_linear_maps(err, sm%map1, branch, 0, 0, orbit = tao_branch%orbit, excite_zero = excite_zero)

call spin_mat_to_eigen (sm%map1%orb_mat, sm%map1%spin_q, orb_eval, orb_evec, n0, spin_evec, err)
call spin_quat_resonance_strengths(orb_evec(j,:), sm%map1%spin_q, xi_sum, xi_diff)

write(1, '(a, f12.8)') '"Polarization Limit ST" ABS 1e-8 ', tao_branch%spin%pol_limit_st
write(1, '(a, f12.8)') '"Polarization Limit DK" ABS 1e-8 ', tao_branch%spin%pol_limit_dk
write(1, '(a, 3f12.8)') '"Polarization Limits DK (a,b,c-modes)" ABS 1e-8 ', tao_branch%spin%pol_limit_dk_partial
write(1, '(a, 3f12.8)') '"Polarization Limits DK (bc,ac,ab-modes)" ABS 1e-8 ', tao_branch%spin%pol_limit_dk_partial2

write(1, '(a, es16.8)') '"Polarization Rate BKS" REL 1e-8 ', tao_branch%spin%pol_rate_bks
write(1, '(a, es16.8)') '"Depolarization Rate" REL 1e-8 ', tao_branch%spin%depol_rate
write(1, '(a, 3es16.8)') '"Depolarization Rate Partial" REL 1e-8 ', tao_branch%spin%depol_rate_partial
write(1, '(a, 3es16.8)') '"Depolarization Rate Partial2" REL 1e-8', tao_branch%spin%depol_rate_partial2

write(1, '(a, es16.8)') '"Integral g^3 * b_hat * n_0" REL 1e-8 ', tao_branch%spin%integral_bn
write(1, '(a, es16.8)') '"Integral g^3 * b_hat * dn/ddelta" REL 1e-8 ', tao_branch%spin%integral_bdn
write(1, '(a, es16.8)') '"Integral g^3 (1 - 2(n * s_hat)/9)" REL 1e-8 ', tao_branch%spin%integral_1ns
write(1, '(a, es16.8)') '"Integral g^3 * 11 (dn/ddelta)^2 / 9" REL 1e-8', tao_branch%spin%integral_dn2

do i = 1, 6
vec = abs(orb_evec(i, :))
select case (i)
case (1, 2); vec(1:4) = vec([1,2,5,6])
case (3, 4); vec(1:4) = vec([1,2,5,6])
case (5, 6); vec(1:4) = vec([1,2,5,6])
end select
write (1, '(a, 6f12.6)') '"orb-evec-' // int_str(i) // '" ABS 1e-8', vec(1:4)
write (1, '(a, 3es16.8)') '"spin-re-evec-' // int_str(i) // '" REL 1e-8', real(spin_evec(i, 1:3:2), 8)
write (1, '(a, 3es16.8)') '"spin-im-evec-' // int_str(i) // '" REL 1e-8', aimag(spin_evec(i, 1:3:2))
enddo

do i = 1, 3
j = 2 * i - 1
call spin_quat_resonance_strengths(orb_evec(j,:), sm%map1%spin_q, xi_sum, xi_diff)
write (1, '(a, f13.7, 2(f17.7, es13.5))') '"Res-Strength-' // int_str(i) // '" REL 1e-8', xi_sum, xi_diff
enddo

!----------------------------

close (1)

end program
47 changes: 12 additions & 35 deletions sim_utils/math/rotation_3d_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,57 +27,30 @@ module rotation_3d_mod
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
!+
! Subroutine w_mat_to_axis_angle (w_mat, axis, angle, ref_axis)
! Subroutine w_mat_to_axis_angle (w_mat, axis, angle)
!
! Routine to find the rotation axis and rotation angle corresponding to a given
! 3D rotation matrix.
!
! The rotation axis is chosen to have a non-negative dot product with the
! reference axis ref_axis.
!
! The rotation angle is chosen in the range [-pi, pi].
! The rotation angle is chosen in the range [0, pi].
!
! Input:
! w_mat(3,3) -- real(rp): Rotation matrix.
! ref_axis(3) -- real(rp), optional: Reference axis. Default is [1, 1, 1].
!
! Output:
! axis(3) -- real(rp): Rotation axis. Normalized to 1.
! angle -- real(rp): Rotation angle in the range [-pi, pi].
! angle -- real(rp): Rotation angle in the range [0, pi].
!-

subroutine w_mat_to_axis_angle (w_mat, axis, angle, ref_axis)
subroutine w_mat_to_axis_angle (w_mat, axis, angle)

real(rp) w_mat(3,3), axis(3), angle
real(rp), optional :: ref_axis(3)
real(rp) sin2_ang, quat(0:3)

!

quat = w_mat_to_quat(w_mat)
if (all(quat(1:3) == 0)) then
axis = [1, 0, 0]
angle = 0
return
endif

!

sin2_ang = norm2(quat(1:3))
axis = quat(1:3) / sin2_ang
angle = modulo2(2 * atan2(sin2_ang, quat(0)), pi)

! Align to axis to point in the general direction of (1,1,1)

if (present(ref_axis)) then
if (dot_product(ref_axis, axis) < 0) then
axis = -axis
angle = -angle
endif
elseif (sum(axis) < 0) then
axis = -axis
angle = -angle
endif
call quat_to_axis_angle(quat, axis, angle)

end subroutine w_mat_to_axis_angle

Expand Down Expand Up @@ -303,14 +276,15 @@ end function omega_to_quat
!+
! Subroutine quat_to_axis_angle (quat, axis, angle)
!
! Routine to convert from axis + angle representation to a quaternion.
! Routine to convert from quaternion to axis + angle representation.
! The angle will be in the range 0 <= angle <= pi.
!
! Input:
! quat(0:3) -- real(rp): Rotation quaternion. Assumed normalized.
!
! Output:
! axis(3) -- real(rp): Axis of rotation.
! angle -- real(rp): angle of rotation.
! angle -- real(rp): angle of rotation in range [0, pi].
!-

subroutine quat_to_axis_angle (quat, axis, angle)
Expand All @@ -325,9 +299,12 @@ subroutine quat_to_axis_angle (quat, axis, angle)
if (anorm == 0) then
angle = 0
axis = [0, 0, 1] ! Arbitrary.
else
elseif (quat(0) > 0) then
angle = 2 * atan2(anorm, quat(0))
axis = quat(1:3) / anorm
else
angle = 2 * atan2(anorm, -quat(0))
axis = -quat(1:3) / anorm
endif

end subroutine quat_to_axis_angle
Expand Down
Loading

0 comments on commit dc5d0a5

Please sign in to comment.