More ring tutorial devel. (#1114)

* More ring tutorial devel.

bmad-doc/tao_examples/custom_tao_with_measured_data/RRNOVAMU2E11172016.bmad

@@ -5500,7 +5500,7 @@ mps638au: multipole, type = "RGF036_U_SHIM", k1l = k1lendsfu+k1lfromk2lfendu, k2
  ffac*larcmag*dir*arck2f-0.1021426e-2*dir, k3l = 0.1297873
 mps638ad: multipole, type = "RGF036_D_SHIM", k1l = -0.221980598e-3, k2l = &
  0.01609515*dir-0.1021426e-2*dir, k3l = 0.1297873
- 
+
 mp638bn: multipole, type = "SGF073_N_BODY", k3l = -0.017814, k4l = 4.4084*dir
 mp638bs: multipole, type = "SGF073_S_BODY", k1l = 0.24e-4, t1 = 0.785398163397, k2l = &
  -0.366e-3*dir, t2 = 0.523598775598, k3l = 0.022904, t3 = 0.392699081699, k4l = &

bmad-doc/tutorial_ring_design/doc/spin.tex

@@ -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

bmad-doc/tutorial_ring_design/doc/tutorial_ring_design.tex

@@ -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}
 
 %----------------------------------------------------------
@@ -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}
@@ -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. 
 %
@@ -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
 
 %------------------------------------------------------------------------------

bmad-doc/tutorial_ring_design/lattices/12_DynamicAperture/solution/ring_w_correction.bmad

@@ -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
 
 
 

bmad-doc/tutorial_ring_design/lattices/15_SpinTracking/DepolTao.py

@@ -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()

bmad-doc/tutorial_ring_design/lattices/15_SpinTracking/solution/long_term_tracking.init

@@ -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 
+/
+/