final beambeam. (#1132)

* Fix beambeam tracking case with sig_z = 0.

bmad/doc/elements.tex

@@ -314,7 +314,8 @@ \section{BeamBeam}
   ks             = <Real>    ! Solenoid strength.
   bs_field       = <Real>    ! Solenoid field strength.
   field_master   = <T/F>     ! Is ks or bs_field value the master (\sref{s:field.master})?
-  z_crossing     = <Real>    ! Weak particle phase space z when strong beam center reaches element.
+  s_beta_ref     = <Real>    ! Reference position of strong beam Twiss.
+  z_crossing     = <Real>    ! Weak particle phase space z when strong beam center reaches IP.
   repetition_frequency = <Real>  ! Strong beam repetition rate.
 \end{example}
 
@@ -347,13 +348,14 @@ \section{BeamBeam}
 
 \index{sig_x}\index{sig_y}
 \index{beta_a_strong}\index{beta_b_strong}\index{alpha_a_strong}\index{alpha_b_strong}
-\vn{sig_x}, \vn{sig_y} are the transverse sigmas of the strong beam at $s_0$ which is the
-$s$-position where the \vn{beambeam} element is located.  For calculating the sigmas of any given
-slice, \vn{sig_x} and \vn{sig_y} are extrapolated using the Twiss parameters at $s_0$. The Twiss
-parameters at $s_0$ are set by \vn{beta_a_strong}, \vn{beta_b_strong}, \vn{alpha_a_strong}, and
-\vn{alpha_b_strong}.  If \vn{beta_a_strong} is zero (the default), the $a$-mode Twiss parameters as
-calculated from the lattice is used. Similarly, if \vn{beta_b_strong} is zero (the default), the
-$b$-mode Twiss parameters as calculated from the lattice is used.
+The strong beam Twiss parameters \vn{beta_a_strong}, \vn{beta_b_strong}, \vn{alpha_a_strong}, and
+\vn{alpha_b_strong} are the Twiss parameters of the strong beam at the $s$-position given by
+\vn{s_twiss_ref}. Additionally, \vn{sig_x}, \vn{sig_y} are the transverse sigmas of the strong beam
+at this point. \vn{S_beta_ref} is measured relative to the position of the\vn{beambeam} element.  If
+\vn{beta_a_strong} is zero (the default), the $a$-mode Twiss parameters as calculated from the
+lattice is used. Similarly, if \vn{beta_b_strong} is zero (the default), the $b$-mode Twiss
+parameters as calculated from the lattice is used.  To calculate the sigmas of any given slice,
+\vn{sig_x} and \vn{sig_y} are extrapolated using the Twiss parameters at \vn{s_twiss_ref}.
 
 \index{x_offset}\index{y_offset}\index{z_offset}
 \index{x_pitch}\index{y_pitch}\index{tilt}
@@ -405,11 +407,12 @@ \section{BeamBeam}
 strong beam. This is only relevant if the velocity of the strong beam is not equal to the velocity
 of the weak beam.
 
-The \vn{z_crossing} parameter sets when in time the center of the strong beam is at the plane of the
-beambeam element. For example, if tracking is done with radiation damping on, the closed orbit will
-have a finite phase space $z$ value at the \vn{beambeam} element. To have the weak beam and strong beam
-centers cross the plane of the \vn{beambeam} element at the same time, the value of \vn{z_crossing} should
-be set to the value of the weak beam $z$.
+The \vn{z_crossing} parameter sets where the center of the strong beam is relative to the plane of the
+beambeam element (the IP) at the time when a weak particle with $z = 0$ is at the IP. For example,
+if tracking is done with radiation damping on, the (weak beam) closed orbit will have a finite phase
+space $z$ value at the \vn{beambeam} element. To have the weak beam and strong beam centers cross
+the plane of the \vn{beambeam} element at the same time, the value of \vn{z_crossing} should be set
+to the value of the weak beam closed orbit $z$ at the IP.
 
 When with absolute time tracking (\sref{s:rf.time}) is in use, the \vn{repetition_frequency}
 parameter (along with the \vn{z_crossing} parameter) is used to calculate the time that the strong

bmad/low_level/bbi_slice_calc.f90

@@ -1,7 +1,8 @@
 !+
 ! Subroutine bbi_slice_calc (ele, n_slice, z_slice)
 !
-! Routine to compute the longitudinal positions of the slices of a beambeam element.
+! Routine to compute the longitudinal positions of the slices of a beambeam element in the
+! beambeam element frame of referece. 
 !
 ! Input:
 !   ele        -- ele_struct: beambeam element
@@ -24,17 +25,17 @@ subroutine bbi_slice_calc (ele, n_slice, z_slice)
 real(rp) z_slice(:), y
 real(rp) :: z_norm
 
-!
+! Since calc is in beambeam element frame of reference, z_slice is independent of ele%value(z_offset$).
 
 z_slice = 0  
 
 if (n_slice == 1) then
-  z_slice(1) = ele%value(z_offset$)
+  z_slice(1) = 0
 elseif (n_slice > 1) then
   do i = 1, n_slice
     y = (i - 0.5) / n_slice - 0.5
     z_norm = inverse(probability_funct, y, -5.0_rp, 5.0_rp, 1.0e-5_rp)
-    z_slice(i) = ele%value(sig_z$) * z_norm + ele%value(z_offset$)
+    z_slice(i) = ele%value(sig_z$) * z_norm
   enddo
 else
   print *, 'ERROR IN BBI_SLICE_CALC: N_SLICE IS NEGATIVE:', n_slice

bmad/low_level/strong_beam_sigma_calc.f90

@@ -23,13 +23,16 @@ subroutine strong_beam_sigma_calc (ele, s_pos, sigma, bbi_const, dsigma_ds)
 type (ele_struct) ele
 type (ele_struct), pointer :: ele0
 
-real(rp) s_pos, bbi_const, x_center, y_center, sigma(2), dsigma_ds(2)
+real(rp) s_pos, bbi_const, x_center, y_center, sigma(2), dsigma_ds(2), ds
 real(rp) beta, sig_x0, sig_y0, beta_a0, beta_b0, alpha_a0, alpha_b0, gamma0
 
 !
 
 sig_x0 = ele%value(sig_x$)
 sig_y0 = ele%value(sig_y$)
+ds = s_pos - ele%value(s_twiss_ref$)
+
+!
 
 if (ele%value(beta_a_strong$) == 0) then
   ele0 => pointer_to_next_ele(ele, -1)
@@ -40,6 +43,18 @@ subroutine strong_beam_sigma_calc (ele, s_pos, sigma, bbi_const, dsigma_ds)
   alpha_a0 = ele%value(alpha_a_strong$)
 endif
 
+if (beta_a0 == 0) then   ! Can happen when testing.
+  sigma(1) = sig_x0
+  dsigma_ds(1) = 0
+else
+  gamma0 = (1 + alpha_a0**2) / beta_a0
+  beta = beta_a0 - 2 * alpha_a0 * ds + gamma0 * ds**2
+  sigma(1) = sig_x0 * sqrt(beta / beta_a0)
+  dsigma_ds(1) = -(alpha_a0 - ds * gamma0) * sigma(1) / beta
+endif
+
+!
+
 if (ele%value(beta_b_strong$) == 0) then
   ele0 => pointer_to_next_ele(ele, -1)
   beta_b0 = ele%b%beta
@@ -49,21 +64,18 @@ subroutine strong_beam_sigma_calc (ele, s_pos, sigma, bbi_const, dsigma_ds)
   alpha_b0 = ele%value(alpha_b_strong$)
 endif
 
-if (beta_a0 == 0) then
-  sigma = [sig_x0, sig_y0]
-  dsigma_ds = 0
+if (beta_b0 == 0) then   ! Can happen when testing.
+  sigma(2) = sig_y0
+  dsigma_ds(2) = 0
 else
-  gamma0 = (1 + alpha_a0**2) / beta_a0
-  beta = beta_a0 - 2 * alpha_a0 * s_pos + gamma0 * s_pos**2
-  sigma(1) = sig_x0 * sqrt(beta / beta_a0)
-  dsigma_ds(1) = -(alpha_a0 - s_pos * gamma0) * sigma(1) / beta
-
   gamma0 = (1 + alpha_b0**2) / beta_b0
-  beta = beta_b0 - 2 * alpha_b0 * s_pos + gamma0 * s_pos**2
+  beta = beta_b0 - 2 * alpha_b0 * ds + gamma0 * ds**2
   sigma(2) = sig_y0 * sqrt(beta / beta_b0)
-  dsigma_ds(2) = -(alpha_b0 - s_pos * gamma0) * sigma(2) / beta
+  dsigma_ds(2) = -(alpha_b0 - ds * gamma0) * sigma(2) / beta
 endif
 
+!
+
 bbi_const = -strong_beam_strength(ele) * classical_radius_factor /  &
                                       (2 * pi * ele%value(p0c$) * (sigma(1) + sigma(2)))
 

bmad/low_level/track_a_beambeam.f90

@@ -39,9 +39,9 @@ subroutine track_a_beambeam (orbit, ele, param, track, mat6, make_matrix)
 
 real(rp), optional :: mat6(6,6)
 real(rp) sigma(2), dsigma_ds(2), ff, ds_slice
-real(rp) xmat(6,6), del_s, bbi_const, dx, dy, dcoef, z, d, coef, s0
+real(rp) xmat(6,6), del_s, bbi_const, dx, dy, dcoef, z, d, coef, ds
 real(rp) om(3), quat(0:3), beta_strong, s_body_save, p_rel, nk(2), dnk(2,2)
-real(rp) f_factor, new_beta, px_old, py_old, e_factor, ds_dz
+real(rp) f_factor, new_beta, px_old, py_old, e_factor, ds_dz, s0
 real(rp), allocatable :: z_slice(:)
 real(rp), target :: slice_center(3), s_body, s_lab
 real(rp), pointer :: center_ptr(:), s_body_ptr, s_lab_ptr ! To get around ifort bug.
@@ -78,8 +78,6 @@ subroutine track_a_beambeam (orbit, ele, param, track, mat6, make_matrix)
 else
   beta_strong = ele%value(p0c$) / ele%value(E_tot$)
 endif
-s0 = 0.1_rp * ele%value(sig_z$) + &
-        abs(particle_rf_time(orbit, ele, rf_freq = ele%value(repetition_frequency$)) * c_light * beta_strong)
 
 call offset_particle (ele, set$, orbit, s_pos = s_lab, s_out = s_body, set_spin = .true., mat6 = mat6, make_matrix = make_matrix)
 
@@ -95,7 +93,9 @@ subroutine track_a_beambeam (orbit, ele, param, track, mat6, make_matrix)
   final_calc = .false.; make_mat = .false.
   s_body_save = s_body
   orb_save = orbit
-  s_lab = super_zbrent(at_slice_func, -abs(z)-s0, abs(z)+s0, 1e-12_rp, 1e-12_rp, status)
+  s0 = 0.5_rp * (orbit%vec(5) + ele%value(z_offset_tot$) - ele%value(z_crossing$) + slice_center(3))
+  ds = abs(s_body - 0.5_rp * slice_center(3)) + 0.5_rp * abs(orbit%vec(5) - ele%value(z_crossing$) - ele%value(z_offset$))
+  s_lab = super_zbrent(at_slice_func, s0-ds, s0+ds, 1e-12_rp, 1e-12_rp, status)
 
   final_calc = .true.; make_mat = logic_option(.false., make_matrix)
   s_body = s_body_save
@@ -181,14 +181,22 @@ subroutine track_a_beambeam (orbit, ele, param, track, mat6, make_matrix)
 !-------------------------------------------------------
 contains
 
+!+
+! Start in body frame
+! Convert to lab frame
+! Move to s_lab_target
+! Convert back to body frame
+! Return ds_slice which is the distance between the particle and the strong slice.
+!-
+
 function at_slice_func(s_lab_target, status) result (ds_slice)
 
 real(rp), intent(in) :: s_lab_target
 real(rp) ds_slice
 real(rp) s_body_target, s_lab_slice, del_s, s_lab, s_beam_center_strong, s_weak
 integer status
 
-!
+! 
 
 call offset_particle(ele, unset$, orbit, s_pos = s_body, s_out = s_lab, set_spin = final_calc, mat6 = mat6, make_matrix = make_mat)
 del_s = s_lab_target - s_lab

bmad/modules/attribute_mod.f90

@@ -879,6 +879,7 @@ subroutine init_attribute_name_array ()
 call init_attribute_name1 (beambeam$, crab_x5$,                     'CRAB_X5')
 call init_attribute_name1 (beambeam$, crab_tilt$,                   'CRAB_TILT')
 call init_attribute_name1 (beambeam$, z_crossing$,                  'Z_CROSSING')
+call init_attribute_name1 (beambeam$, s_twiss_ref$,                 'S_TWISS_REF')
 call init_attribute_name1 (beambeam$, repetition_frequency$,        'REPETITION_FREQUENCY')
 call init_attribute_name1 (beambeam$, rf_clock_harmonic$,           'rf_clock_harminic', private$)
 call init_attribute_name1 (beambeam$, species_strong$,              'SPECIES_STRONG')
@@ -2041,7 +2042,7 @@ function attribute_units (attrib_name, unrecognized_units) result (attrib_units)
       'X_REF', 'Y_REF', 'Z', 'Z0', 'Z1', 'Z_OFFSET', 'Z_OFFSET_TOT', 'Z_POSITION', 'Z_REF', 'MOSAIC_THICKNESS', &
       'C12_MAT0', 'C12_MAT1', 'X_GAIN_CALIB', 'Y_GAIN_CALIB', 'X_GAIN_ERR', 'Y_GAIN_ERR', 'RADIUS', &
       'Z_APERTURE_WIDTH2', 'Z_APERTURE_CENTER', 'RF_WAVELENGTH', 'Z_CROSSING', &
-      'X1_EDGE', 'X2_EDGE', 'Y1_EDGE', 'Y2_EDGE', 'L_RECTANGLE', &
+      'X1_EDGE', 'X2_EDGE', 'Y1_EDGE', 'Y2_EDGE', 'L_RECTANGLE', 'S_TWISS_REF', &
       'X_DISPERSION_ERR', 'Y_DISPERSION_ERR', 'X_DISPERSION_CALIB', 'Y_DISPERSION_CALIB')
   attrib_units = 'm'
 

bmad/modules/bmad_struct.f90

@@ -1711,7 +1711,7 @@ module bmad_struct
 integer, parameter :: tilt$ = 2, roll$ = 2, n_part$ = 2, inherit_from_fork$ = 2 ! Important: tilt$ = roll$
 integer, parameter :: ref_tilt$ = 3, direction$ = 3, repetition_frequency$ = 3, deta_ds_master$ = 3, &
                       kick$ = 3, x_gain_err$ = 3, taylor_order$ = 3, r_solenoid$ = 3, final_charge$ = 3
-integer, parameter :: k1$ = 4, kx$ = 4, harmon$ = 4, h_displace$ = 4, y_gain_err$ = 4, &
+integer, parameter :: k1$ = 4, kx$ = 4, harmon$ = 4, h_displace$ = 4, y_gain_err$ = 4, s_twiss_ref$ = 4, &
                       critical_angle_factor$ = 4, tilt_corr$ = 4, ref_coords$ = 4, dt_max$ = 4
 integer, parameter :: graze_angle$ = 5, k2$ = 5, b_max$ = 5, v_displace$ = 5, gradient_tot$ = 5, harmon_master$ = 5, &
                       ks$ = 5, flexible$ = 5, crunch$ = 5, ref_orbit_follows$ = 5, pc_out_min$ = 5

bmad/output/type_ele.f90

@@ -1622,7 +1622,7 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c
 character(40) a_name, a2_name
 logical is_2nd_col_attrib
 
-character(41), parameter :: att_name(96) = [character(40):: 'X_PITCH', 'Y_PITCH', 'X_OFFSET', &
+character(41), parameter :: att_name(97) = [character(40):: 'X_PITCH', 'Y_PITCH', 'X_OFFSET', &
                 'Y_OFFSET', 'Z_OFFSET', 'REF_TILT', 'TILT', 'ROLL', 'X1_LIMIT', 'Y1_LIMIT', &
                 'FB1', 'FQ1', 'LORD_PAD1', 'HKICK', 'VKICK', 'KICK', 'FRINGE_TYPE', 'DS_STEP', 'R0_MAG', &
                 'KS', 'K1', 'K2', 'G', 'DG', 'G_TOT', 'H1', 'E1', 'FINT', 'HGAP', &
@@ -1636,9 +1636,9 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c
                 'C11_MAT0', 'C12_MAT0', 'C21_MAT0', 'C22_MAT0', 'HARMON', 'FINAL_CHARGE', &
                 'MODE_FLIP0', 'BETA_A_STRONG', 'BETA_B_STRONG', 'REF_TIME_START', 'THICKNESS', &
                 'PX_KICK', 'PY_KICK', 'PZ_KICK', 'E_TOT_OFFSET', 'FLEXIBLE', 'CRUNCH', 'NOISE', &
-                'F_FACTOR', 'EXACT_MULTIPOLES']
+                'F_FACTOR', 'EXACT_MULTIPOLES', 'Z_CROSSING']
 
-character(41), parameter :: att2_name(96) = [character(40):: 'X_PITCH_TOT', 'Y_PITCH_TOT', 'X_OFFSET_TOT', &
+character(41), parameter :: att2_name(97) = [character(40):: 'X_PITCH_TOT', 'Y_PITCH_TOT', 'X_OFFSET_TOT', &
                 'Y_OFFSET_TOT', 'Z_OFFSET_TOT', 'REF_TILT_TOT', 'TILT_TOT', 'ROLL_TOT', 'X2_LIMIT', 'Y2_LIMIT', &
                 'FB2', 'FQ2', 'LORD_PAD2', 'BL_HKICK', 'BL_VKICK', 'BL_KICK', 'FRINGE_AT', 'NUM_STEPS', 'R0_ELEC', &
                 'BS_FIELD', 'B1_GRADIENT', 'B2_GRADIENT', 'B_FIELD', 'DB_FIELD', 'B_FIELD_TOT', 'H2', 'E2', 'FINTX', 'HGAPX', &
@@ -1652,7 +1652,7 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c
                 'C11_MAT1', 'C12_MAT1', 'C21_MAT1', 'C22_MAT1', 'HARMON_MASTER', 'SCATTER', &
                 'MODE_FLIP1', 'ALPHA_A_STRONG', 'ALPHA_B_STRONG', 'DELTA_REF_TIME', 'DTHICKNESS_DX', &
                 'X_KICK', 'Y_KICK', 'Z_KICK', 'E_TOT_START', 'REF_COORDS', 'CRUNCH_CALIB', 'N_SAMPLE', &
-                'SCATTER_METHOD', 'FIDUCIAL_PT']
+                'SCATTER_METHOD', 'FIDUCIAL_PT', 'S_BETA_MIN']
 
 ! Exceptional cases