Skip to content

Commit

Permalink
More ring tutorial devel.
Browse files Browse the repository at this point in the history
  • Loading branch information
DavidSagan committed Aug 2, 2024
1 parent 482af37 commit 6550b21
Show file tree
Hide file tree
Showing 8 changed files with 387 additions and 24 deletions.
18 changes: 18 additions & 0 deletions bmad-doc/tutorial_ring_design/doc/spin.tex
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
%------------------------------------------------------------------------------
%------------------------------------------------------------------------------

%------------------------------------------------------------------------------
%\subsection{Tracking with Ramping}
%The ESR will not be used to accelerate, however to show what ramping in Bmad might look like we include this section.

%------------------------------------------------------------------------------
%\subsection{Adding Siberian Snakes}
%Taylor element overview



%------------------------------------------------------------------------------
%\subsection{Invariant Spin Field Calculations}
%\sodom $\rightarrow$ \ltt

\newpage
115 changes: 110 additions & 5 deletions bmad-doc/tutorial_ring_design/doc/tutorial_ring_design.tex
Original file line number Diff line number Diff line change
Expand Up @@ -2299,7 +2299,7 @@ \subsection{Example Lattice}
is used, only one control variable (in this case \vn{hh}) may be used by the overlay.
\item
The initial values for control variables may be set when defining the control element. For example,
\vn{hh} of \vn{ov2} is set to 0.2. Control variables default to a value of zero.
\vn{hh} of \vn{ov2} is set to 0.01. Control variables default to a value of zero.
\end{itemize}

%----------------------------------------------------------
Expand Down Expand Up @@ -3037,9 +3037,7 @@ \subsection{Example: Adding Corrector Coils and BPMs}
CH: hkicker, L = 0.2
CV: vkicker, L = 0.2
\end{code}
\item Add corrector coils near quadrupoles. Then, add horizontal kickers near focusing quadrupoles and vertical kickers near defocusing quadrupoles.

After adding the kicker, correctors in the D1 drift have to be added. This means placing the correctors after quadrupoles in the forward geometry and before quadrupoles in the reverse geometry.
\item Then, add horizontal kickers near focusing quadrupoles and vertical kickers near defocusing quadrupoles. Correctors are added in the D1 drifts. This means placing the correctors after quadrupoles in the forward geometry and before quadrupoles in the reverse geometry.

Because we have used D1 drifts repeatedly throughout the ring, adding correctors manually is a tedious task. The easiest approach is to use \vn{regular expressions} in Find and Replace. The recipe is the following:
\begin{code}
Expand Down Expand Up @@ -3187,7 +3185,7 @@ \subsection{Example: Orbit Correction with Quadrupole Misalignments}
%
\item Correct the orbit using the same optimization process. You can reuse the \vn{tao.init} we wrote for example \sref{s:sawtooth}.
%
\item If the misalignment amplitude is large, the orbit could become unstable and the optimizer will have a hard time finding a solution. Try larger misalignments of $\sigma=100\unit{\mu m}$ with "\vn{set ele quadrupole::* x_offset = 1e-5*ran_gauss(3)} and correct the orbit. This will likely fail and you'll see many error messages.
\item If the misalignment amplitude is large, the orbit could become unstable and the optimizer will have a hard time finding a solution. Try larger misalignments of $\sigma=100\unit{\mu m}$ with "\vn{set ele quadrupole::* x_offset = 1e-4*ran_gauss(3)} and correct the orbit. This will likely fail and you'll see many error messages.
%
\item One approach to avoid unstable orbits is correcting the orbit with an open lattice first. Reinitialize \tao to start fresh from the beginning. Use command \vn{cut} to change the lattice geometry from closed to open. \tao now uses the single-turn orbit instead, so it won't complain even if the closed orbit is unstable.
%
Expand Down Expand Up @@ -3218,6 +3216,113 @@ \subsection{Exercises}
%
\end{enumerate}

\newpage

\section{Spin Tracking with Ramping}
\subsection{Introduction}
\bmad is one of the most powerful modern accelerator physics libraries for polarization analyses and spin-tracking simulations. The most important tools of polarization analysis, such as the invariant spin field, the amplitude-dependent spin tune, and spin-orbit resonance strengths, can be computed using \bmad. The invariant spin field (ISF) and the amplitude-dependent spin tune (ADST) can be computed in \bmad using normal form via \tao, K. Yokoya's SODOM-2 algorithm via \sodom, or stroboscopic averaging via \spinstrob, but these applications will not be covered in this brief introduction. In this chapter, we will ramp our storage ring energy through a spin-orbit resonance, calculate the polarization transmission, and use it to verify the resonance strength computed by \tao using the Froissart-Stora formula.

In a perfectly-flat ring (a ring in which the only magnetic fields encountered on the closed orbit are vertical), the closed-orbit spin tune is $\nu_0=G\gamma$, where $G=(g-2)/2$ is the anomalous gyromagnetic $g$-factor and $\gamma$ is the reference gamma. We will ramp our ring through the $\nu_0=60-Q_y \approx 51.81$ resonance.

\subsection*{Example: Ramping Through a Spin-Orbit Resonance}
\begin{enumerate}[leftmargin=*]
\item Start with the \vn{spin_lat.bmad} file, located in \vn{lattices/15_SpinTracking}. It is a FODO ring based on the layout of the AGS.

\item Let's define a \vn{ramper} element which ramps $G\gamma$ from $51.71$ to $51.82$ in $5000$ turns. To ramp the \vn{e_tot} attribute for all elements in such a time, we need to obtain the revolution time, which is done by calling \vn{show uni} and looking for the
\vn{Lattice branch transit time}. For our ring, this is $2.67024$ microseconds. Thus, our \vn{ramper} element is

\begin{code}
rampe: ramper = {*[e_tot] : 51.71 / anomalous_moment_of(proton) * m_proton +
(.2/5000/anomalous_moment_of(proton) * m_proton) / (2.67024e-6)*time}, var = {time}
\end{code}

We can define this element anywhere in the lattice file to initiate the ramping.

\item Create a \vn{long_term_tracking.init} file, containing the following:
\begin{code}
&params
ltt%lat_file = 'spin_lat.bmad'
ltt%ele_start = '0'
ltt%ele_stop = ''
ltt%phase_space_output_file = 'tracking.dat'
ltt%averages_output_file = 'ramp.dat'
ltt%averages_output_every_n_turns = 1
ltt%particle_output_every_n_turns = 1

ltt%ramping_on = T
ltt%ramp_update_each_particle = T
ltt%simulation_mode = 'SINGLE'
ltt%tracking_method = 'BMAD'
ltt%n_turns = 5000
ltt%rfcavity_on = T
ltt%split_bends_for_stochastic_rad = F
ltt%random_seed = 0
ltt%add_closed_orbit_to_init_position = T

bmad_com%spin_tracking_on = T
bmad_com%radiation_damping_on = F
bmad_com%radiation_fluctuations_on = F

beam_init%n_particle = 1
beam_init%spin = 0, 1, 0
beam_init%center = 0, 0, 1e-5, 0, 2.66, 1.15e-5
/
/
\end{code}

\item Plot $S_y$ versus turn number with the following \textit{gnuplot} command:
\begin{code}
plot 'tracking.dat' u 1:13
\end{code}
Observe the asymptotic value of $S_y$ as the number of turns increases.

\item Load Tao, and set the energy such that $G\gamma=51.81$ with the command
\begin{code}
set ele * e_tot = 51.81 / anomalous_moment_of(proton) * mass_of(proton)
\end{code}
Then use the command \vn{show spin -ele 0}
to observe the resonance strength. We are working with a B-mode (vertical) resonance, and whether $\text{frac}(Q_s+Q)$ or $\text{frac}(Q_s-Q)$ is smaller tells you whether to use \vn{Xi_sum} or \vn{Xi_diff}. The correct resonance strength corresponds to whichever $\text{frac}(\dots)$ expression is closer to zero. The resonance strength given by \tao is normalized, so we have to multiply by $\sqrt{J_y}$ to find the true resonance strength. We can obtain $J_y$ using the equation
\begin{equation}
y(s)=\sqrt{2J_y\beta_y(s)}\cos[\Phi_y(s)],
\end{equation}
along with the command \vn{show twiss}.

\item We can find what the asymptotic value of $S_y$ should have been using the Froissart-Stora formula:
\begin{equation}
S_y(\theta\to\infty) = 2\exp\left(-\frac{\pi\epsilon^2}{2\tilde{\alpha}}\right)-1.
\end{equation}
Here, $\theta$ is the machine azimuth, $\epsilon$ is the resonance strength, and
\begin{equation}
\tilde{\alpha} = \frac{\mathrm{d}}{\mathrm{d}\theta}G\gamma = \frac{1}{2\pi} \times [\text{Change in }G\gamma\text{ per turn}].
\end{equation}
Compare the result of the Froissart-Stora formula to the result obtained from tracking.

\end{enumerate}

\subsection{Exercises}
\begin{enumerate}[leftmargin=*]
\item The AGS has two special spin rotators known as partial Siberian snakes. Add the partial Siberian snakes to \vn{spin_lat.bmad}. They are represented by elements called CSNK (cold snake) and WSNK (warm snake), which are markers in the lattice we have provided to you. We will represent the ideal (point-like) Siberian snakes using \vn{Taylor} elements. Here is the example syntax:
\begin{code}
CSNK: Taylor, {S1: q_0|}, {Sx: q_1|}, {Sy: q_2|}, {Sz: q_3|}
\end{code}
In this case, the \vn{Taylor} element produces an instantaneous rotation represented by the unit quaternion $(q_0,q_1,q_2,q_3)$. Recall that the unit quaternion $(\cos(\vartheta/2),\sin(\vartheta/2)\vec{e})$ produces a rotation by angle $\vartheta$ around the unit vector $\vec{e}$. The cold snake has a rotation angle of $18.3^\circ$, and the warm snake has a rotation angle of $10.57^\circ$. Both snakes rotate the spin around the longitudinal axis.

Play around with the \vn{show spin} command at different energies, and verify that the spin tune is \emph{not} $G\gamma$. Why is this the case?

\item If you have done the previous exercise, comment out the \vn{Taylor} elements and return the snakes to markers for this exercise. Use the script \vn{DepolTao.py} to plot the resonance spectrum in the range $G\gamma = 30$ to $G\gamma =60$. It is run using the following syntax:
\begin{code}
python DepolTao.py <lat_file> <output_file> <plot_file> <particle> <SPRINT>
<min Ggamma> <max Ggamma> <J_y>
\end{code}
For this case, set <SPRINT> to T. This boolean flag can take the values T or F, and if T, the spin tracking mode is set to SPRINT. SPRINT tracking is fast, but can be inaccurate when magnets are misaligned or there are strong nonlinearities.

What do you observe about the spectrum? Can you explain any of its features? \emph{Hint: look at the structure of the lattice.}
\end{enumerate}





\newpage

%------------------------------------------------------------------------------
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -272,26 +272,27 @@ RF0[voltage] = 3557943.69553301670



QFSS_2[K1] = 3.32203875341996524E-001
QDSS_2[K1] = -3.82037558156482993E-001
QFF2_2[K1] = 3.59013350200977244E-001
QDF2_2[K1] = -3.59058256797260722E-001
QFF3_2[K1] = 3.70341883503764902E-001
QDF3_2[K1] = -3.82980378703924262E-001
QFF3_6[K1] = 3.78364742603524318E-001
QDF3_6[K1] = -3.60005652920728170E-001
QEF1[K1] = 4.70887920003225435E-001
QEF2[K1] = -1.06759448298209547E-001
QEF3[K1] = 1.09247965928191637E-001
QEF4[K1] = -2.24282292780295017E-001
QDR2_6[K1] = -3.74852087251277888E-001
QFR3_6[K1] = 3.72072224553677777E-001
QDR3_6[K1] = -3.78318987293625675E-001
QER1[K1] = -3.81490113592264146E-001
QER2[K1] = -6.98044300333951120E-003
QER3[K1] = 2.34941240777968713E-001
QER4[K1] = -2.37161404652608571E-001

QFSS_2[K1] = 3.37974220242816981E-001
QDSS_2[K1] = -3.86068025281473437E-001
QFF2_2[K1] = 3.58163349981193480E-001
QDF2_2[K1] = -3.59534202666213554E-001
QFF3_2[K1] = 3.73186818651336460E-001
QDF3_2[K1] = -3.85507583316477442E-001
QDF2_6[K1] = -3.43367011901662356E-001
QFF3_6[K1] = 3.82685494718714492E-001
QDF3_6[K1] = -4.70125485018425049E-001
QEF1[K1] = 4.98569440085783133E-001
QEF2[K1] = -7.49160193280077602E-002
QEF3[K1] = 1.08628573978235388E-001
QEF4[K1] = -2.22875512209623011E-001
QDR2_6[K1] = -3.75207659158594498E-001
QFR3_6[K1] = 3.60552988686989928E-001
QDR3_6[K1] = -3.75816347627969616E-001
QER1[K1] = -3.78378195523198291E-001
QER2[K1] = -1.45851666549070350E-002
QER3[K1] = 2.34797351415774680E-001
QER4[K1] = -2.37167456221783396E-001



Expand Down
90 changes: 90 additions & 0 deletions bmad-doc/tutorial_ring_design/lattices/15_SpinTracking/DepolTao.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
from pytao import Tao
import numpy as np
import matplotlib.pyplot as plt
import os
import time
import math
import sys
import csv



def generate_spin_resonant_tunes(NU, Ki, Kf):
# Initialize a dictionary to hold valid K values and their annotations
valid_K_values_with_annotations = {}

# Estimate a reasonable range for N
N_min = int(Ki - abs(NU)) - 1
N_max = int(Kf + abs(NU)) + 1

for N in range(N_min, N_max):
# Calculate potential K values and check if they meet the criteria
K_values = [
(NU + N, 2 * math.pi),
(abs(NU - N), -2 * math.pi if N > NU else 2 * math.pi),
(N + NU, 2 * math.pi),
]

for K, annotation in K_values:
if Ki <= K <= Kf and K > 0 and K not in valid_K_values_with_annotations:
valid_K_values_with_annotations[K] = annotation

# Sort the dictionary by key (K value) to return a sorted list of tuples
sorted_K_values_and_annotations = sorted(valid_K_values_with_annotations.items())

# Return the sorted list of unique K values and their annotations
return sorted_K_values_and_annotations

def generate_strengths(inf,out,plot_file,particle,SPRINT,lo,hi,J):
tao=Tao(f'-lat {inf} -noplot')
#tao.cmd('set element overlay::* is_on = F')
#tao.cmd('set element * field_master = F')
if SPRINT:
tao.cmd('set ele * spin_tracking_method = sprint')
tao.cmd('set bmad_com spin_tracking_on = T')
#tao.cmd('show -write bmad.twiss lat beginning:end -att l -att ang -att beta_a -att alpha_a -att beta_b -att alpha_b -att phi_b/(2*pi) -att orbit_y -att k1*l -att k2*l -att h1 -att h2 -att orbit_x -att orbit_px -att orbit_py')
data = tao.cmd('show val lat::tune.b[0]/twopi')
Qy =float(data[0])
Kn = generate_spin_resonant_tunes(Qy, lo, hi)

strengths = []

for i in range(0,len(Kn)):
tao.cmd(f'set element 0 e_tot = {Kn[i][0]} / anomalous_moment_of({particle}) * mass_of({particle})')
s = tao.cmd('python spin_resonance')
d1 = float(s[6][s[6].find('dq_b_sum;REAL;F;')+len('dq_b_sum;REAL;F;'):])
d2 = float(s[7][s[7].find('dq_b_diff;REAL;F;')+len('dq_b_diff;REAL;F;'):])
e1 = float(s[8][s[8].find('xi_res_b_sum;REAL;F;')+len('xi_res_b_sum;REAL;F;'):])
e2 = float(s[9][s[9].find('xi_res_b_diff;REAL;F;')+len('xi_res_b_diff;REAL;F;'):])
if abs(d1) > abs(d2):
strengths.append(e2*np.sqrt(J))
else:
strengths.append(e1*np.sqrt(J))

lst = zip([Kn[i][0] for i in range(len(Kn))],strengths)
open(out, 'w').close()
with open(out,"w") as f:
writer = csv.writer(f,delimiter='\t')
writer.writerows(lst)

plt.scatter([k[0] for k in Kn],strengths)
plt.xlabel(r"$G\gamma$")
plt.ylabel("$\epsilon$")
plt.savefig(plot_file, dpi=1200)


def main():
inf = sys.argv[1]
out = sys.argv[2]
plot_file = sys.argv[3]
particle = sys.argv[4]
SPRINT = True if sys.argv[5].upper() == "T" else False
lo = float(sys.argv[6])
hi = float(sys.argv[7])
J = float(sys.argv[8])
generate_strengths(inf,out,plot_file,particle,SPRINT,lo,hi,J)



if __name__ == "__main__":
main()
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
&params
ltt%lat_file = 'spin_lat.bmad' ! Lattice file
ltt%ele_start = '0' ! Where to start in the lattice
ltt%ele_stop = '' ! Where to start in the lattice
ltt%phase_space_output_file = 'tracking.dat'
ltt%averages_output_file = 'ramp.dat'
ltt%averages_output_every_n_turns = 1
ltt%particle_output_every_n_turns = 1

ltt%ramping_on = T
ltt%ramp_update_each_particle = T
ltt%simulation_mode = 'SINGLE'
ltt%tracking_method = 'BMAD' !
ltt%n_turns = 5000 ! Number of turns to track
ltt%rfcavity_on = T
ltt%split_bends_for_stochastic_rad = F
ltt%random_seed = 0 ! Random number seed. 0 => use system clock.
ltt%add_closed_orbit_to_init_position = T

bmad_com%spin_tracking_on = T ! See Bmad manual for bmad_com parameters.
bmad_com%radiation_damping_on = F
bmad_com%radiation_fluctuations_on = F

beam_init%n_particle = 1
beam_init%spin = 0, 1, 0 ! See Bmad manual for beam_init_struct parameters.
beam_init%center = 0, 0, 1e-5, 0, 2.66, 1.15e-5
/
/
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
parameter[geometry] = Closed

parameter[e_tot] = 51.82 / anomalous_moment_of(proton) * mass_of(proton)
parameter[particle] = Proton
bmad_com[spin_tracking_on] = T

beginning[beta_a] = 10.4690053690602
beginning[alpha_a] = 8.47834323010459E-17
beginning[eta_x] = 1.59632886512312
beginning[etap_x] = -4.15786647847316E-16
beginning[beta_b] = 24.1898493928511
beginning[alpha_b] = 9.8485584267628E-16
particle_start[spin_x] = -8.05383665788996E-9
particle_start[spin_y] = 0.999999999999998
particle_start[spin_z] = 6.54209117387931E-8


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

rampe: ramper = {*[e_tot] : 51.71 / anomalous_moment_of(proton) * mass_of(proton) + (.2/5000 / anomalous_moment_of(proton) * mass_of(proton)) / (2.67024e-6)*time}, var = {time}

Q1H: Quadrupole, L = 1.33333333333333, K1 = -0.0541431006647352
DR: Drift, L = 0.666666666666667
B: SBend, L = 2.66697130844605, G = 0.0196327112309049, E1 = 0.0261799387799149, E2 = 0.0261799387799149
Q2: Quadrupole, L = 2.66666666666667, K1 = 0.054724356256391
RF: Rfcavity, L = 0.0, VOLTAGE = 1e6, harmon = 1
CSNK: Taylor, {S1: cos(18.30*pi/360)|}, {Sx: 0|}, {Sy: 0|}, {Sz: sin(18.30*pi/360)|}
WSNK: Taylor, {S1: cos(10.57*pi/360)|}, {Sx: 0|}, {Sy: 0|}, {Sz: sin(10.57*pi/360)|}

!-------------------------------------------------------
! Overlays, groups, rampers, and superimpose


!-------------------------------------------------------
! Lattice lines


RING: line = ( Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR,
Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B,
DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR,
B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H,
DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H,
Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR,
Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B,
DR, Q1H, CSNK, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR,
Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B,
DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR,
B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H,
DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H,
Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR,
Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B,
DR, Q1H, WSNK, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR,
Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B,
DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR,
B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H,
DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H,
Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR,
Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B, DR, Q1H, Q1H, DR, B, DR, Q2, DR, B,
DR, Q1H, RF)

use, RING
Loading

0 comments on commit 6550b21

Please sign in to comment.