Fix beginning element dispersion bookkeeping.

bmad/code/lat_make_mat6.f90

@@ -119,7 +119,6 @@ recursive subroutine lat_make_mat6 (lat, ix_ele, ref_orb, ix_branch, err_flag)
  n_taylor = 0 ! number of taylor map found
 
  do i = 1, branch%n_ele_track
-
  ele => branch%ele(i)
 
  ! Check if transfer matrix needs to be recomputed

bmad/code/make_mat6.f90

@@ -43,6 +43,13 @@ recursive subroutine make_mat6 (ele, param, start_orb, end_orb, err_flag)
 character(*), parameter :: r_name = 'make_mat6'
 
 !--------------------------------------------------------
+! The beginning element is handled specially.
+! Also see twiss_propagate1.
+
+if (ele%key == beginning_ele$) then
+ return
+endif
+
 ! Some init.
 ! If start_orb is in its not_set state (can happen if a particle is lost in 
 ! tracking and ele is downstream from the loss point), init the orbit to zero.

bmad/code/twiss_propagate1.f90

@@ -30,7 +30,7 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag)
 integer key2, geometry
 
 real(rp) :: mat6(6,6)
-real(rp) v_mat(4,4), v_inv_mat(4,4), det, mat2_a(2,2), mat2_b(2,2)
+real(rp) det, mat2_a(2,2), mat2_b(2,2)
 real(rp) big_M(2,2), small_m(2,2), big_N(2,2), small_n(2,2)
 real(rp) c_conj_mat(2,2), E_inv_mat(2,2), F_inv_mat(2,2)
 real(rp) mat2(2,2), eta1_vec(6), eta_vec(6), vec(6), dpz2_dpz1, rel_p1, rel_p21, rel_p2
@@ -50,10 +50,12 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag)
  if (is_true(ele1%value(deta_ds_master$))) then
  ele1%x%etap = ele1%x%deta_ds * rel_p1 + ele1%map_ref_orb_out%vec(2) / rel_p1
  ele1%y%etap = ele1%y%deta_ds * rel_p1 + ele1%map_ref_orb_out%vec(4) / rel_p1
- elseif (ele1%x%deta_ds == real_garbage$) then
+ else
  ele1%x%deta_ds = ele1%x%etap / rel_p1 - ele1%map_ref_orb_out%vec(2) / rel_p1**2
  ele1%y%deta_ds = ele1%y%etap / rel_p1 - ele1%map_ref_orb_out%vec(4) / rel_p1**2
  endif
+
+ call normal_mode_dispersion(ele1)
 endif
 
 !
@@ -244,17 +246,7 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag)
 ele2%z%etap = 1
 ele2%z%deta_ds = 1
 
-call make_v_mats (ele2, v_mat, v_inv_mat)
-eta_vec(1:4) = matmul (v_inv_mat, eta_vec(1:4))
-vec(1:4) = matmul(v_inv_mat, orb_out%vec(1:4))
-
-ele2%a%eta = eta_vec(1)
-ele2%a%etap = eta_vec(2)
-ele2%a%deta_ds = eta_vec(2) / rel_p2 - vec(2) / rel_p2**2
-
-ele2%b%eta = eta_vec(3)
-ele2%b%etap = eta_vec(4)
-ele2%b%deta_ds = eta_vec(4) / rel_p2 - vec(4) / rel_p2**2
+call normal_mode_dispersion(ele2)
 
 !
 

bmad/doc/cover-page.tex

@@ -3,7 +3,7 @@
 
 \begin{flushright}
 \large
- Revision: August 18, 2024 \\
+ Revision: August 31, 2024 \\
 \end{flushright}
 
 \pdfbookmark[0]{Preamble}{Preamble} 

bmad/doc/linear-optics.tex

@@ -424,7 +424,9 @@ \section{Dispersion Calculation}
 \end{equation}
 and this vector is propagated via
 \begin{equation}
+ \bfeta(s_2) = \bfM_{21} \, \bfeta(s_1)
 \end{equation}
+where $\bfM_{21} is the transfer matrix between points $s_1$ and $s_2$.
 
 For an open geometry lattice branch, there are two ways one can imagine
 defining the dispersion: Either with respect to changes in energy at the beginning of the machine or

bmad/low_level/normal_mode_dispersion.f90

@@ -0,0 +1,63 @@
+!+
+! Subroutine normal_mode_dispersion(ele, reverse) 
+!
+! Routine to calculate the normal mode dispersion from the x,y dispersions.
+! Or vice versa if reverse = True.
+!
+! Input:
+! ele -- ele_struct: Element whose dispersions are to be adjusted.
+! reverse -- logical, optional: Default is False. If True, calculate the x,y dispersions 
+! from the normal mode ones.
+!
+! Output:
+! ele -- ele_struct: Element with adjusted dispersions.
+!-
+
+subroutine normal_mode_dispersion(ele, reverse) 
+
+use bmad_interface
+
+implicit none
+
+type (ele_struct) ele
+
+real(rp) v_mat(4,4), v_inv_mat(4,4), eta_vec(4), orb_vec(4), rel_p
+logical, optional :: reverse
+
+! Normal mode to x,y
+
+if (logic_option(.false., reverse)) then
+
+ call make_v_mats (ele, v_mat = v_mat)
+ eta_vec = [ele%a%eta, ele%a%etap, ele%b%eta, ele%b%etap]
+ eta_vec = matmul (v_mat, eta_vec)
+ orb_vec = ele%map_ref_orb_out%vec(1:4)
+ rel_p = 1 + ele%map_ref_orb_out%vec(6)
+
+ ele%x%eta = eta_vec(1)
+ ele%x%etap = eta_vec(2)
+ ele%x%deta_ds = eta_vec(2) / rel_p - orb_vec(2) / rel_p**2
+
+ ele%y%eta = eta_vec(3)
+ ele%y%etap = eta_vec(4)
+ ele%y%deta_ds = eta_vec(4) / rel_p - orb_vec(4) / rel_p**2
+
+! x,y to normal mode
+
+else
+ call make_v_mats (ele, v_inv_mat = v_inv_mat)
+ eta_vec = [ele%x%eta, ele%x%etap, ele%y%eta, ele%y%etap]
+ eta_vec = matmul (v_inv_mat, eta_vec)
+ orb_vec = matmul(v_inv_mat, ele%map_ref_orb_out%vec(1:4))
+ rel_p = 1 + ele%map_ref_orb_out%vec(6)
+
+ ele%a%eta = eta_vec(1)
+ ele%a%etap = eta_vec(2)
+ ele%a%deta_ds = eta_vec(2) / rel_p - orb_vec(2) / rel_p**2
+
+ ele%b%eta = eta_vec(3)
+ ele%b%etap = eta_vec(4)
+ ele%b%deta_ds = eta_vec(4) / rel_p - orb_vec(4) / rel_p**2
+endif
+
+end subroutine

bmad/modules/bmad_routine_interface.f90

@@ -1580,6 +1580,12 @@ subroutine new_control (lat, ix_ele, ele_name)
  character(*), optional :: ele_name
 end subroutine
 
+subroutine normal_mode_dispersion(ele, reverse) 
+ import
+ type (ele_struct) ele
+ logical, optional :: reverse
+end subroutine
+
 function num_field_eles (ele) result (n_field_ele)
  import
  implicit none

bmad/modules/bmad_struct.f90

@@ -19,7 +19,7 @@ module bmad_struct
 ! IF YOU CHANGE THE LAT_STRUCT OR ANY ASSOCIATED STRUCTURES YOU MUST INCREASE THE VERSION NUMBER !!!
 ! THIS IS USED BY BMAD_PARSER TO MAKE SURE DIGESTED FILES ARE OK.
 
-integer, parameter :: bmad_inc_version$ = 322
+integer, parameter :: bmad_inc_version$ = 323
 
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

bmad/modules/changed_attribute_bookkeeper.f90

@@ -281,7 +281,7 @@ subroutine set_flags_for_changed_real_attribute (ele, attrib, set_dependent)
 
 real(rp), optional, target :: attrib
 real(rp), pointer :: a_ptr
-real(rp) v_mat(4,4), v_inv_mat(4,4), eta_vec(4), eta_xy_vec(4), p0c_factor, ff, rel_p1
+real(rp) p0c_factor, ff, rel_p1
 real(rp), target :: unknown_attrib, dangle
 
 integer i
@@ -526,27 +526,15 @@ subroutine set_flags_for_changed_real_attribute (ele, attrib, set_dependent)
 
  if (associated(a_ptr, ele%x%eta) .or. associated(a_ptr, ele%y%eta) .or. coupling_change) then
  if (dep_set) then
- call make_v_mats (ele, v_mat, v_inv_mat)
- eta_xy_vec = [ele%x%eta, ele%x%etap, ele%y%eta, ele%y%etap]
- eta_vec = matmul (v_inv_mat, eta_xy_vec)
- ele%a%eta = eta_vec(1)
- ele%a%etap = eta_vec(2)
- ele%b%eta = eta_vec(3)
- ele%b%etap = eta_vec(4)
+ call normal_mode_dispersion(ele)
  endif
  return
  endif
 
  if (associated(a_ptr, ele%a%eta) .or. associated(a_ptr, ele%a%etap) .or. &
  associated(a_ptr, ele%b%eta) .or. associated(a_ptr, ele%b%etap)) then 
  if (dep_set) then
- call make_v_mats (ele, v_mat, v_inv_mat)
- eta_vec = [ele%a%eta, ele%a%etap, ele%b%eta, ele%b%etap]
- eta_xy_vec = matmul (v_mat, eta_vec)
- ele%x%eta = eta_xy_vec(1)
- ele%x%etap = eta_xy_vec(2)
- ele%y%eta = eta_xy_vec(3)
- ele%y%etap = eta_xy_vec(4)
+ call normal_mode_dispersion(ele, .true.)
  endif
  return
  endif

bmad/output/type_ele.f90

@@ -117,7 +117,7 @@ subroutine type_ele (ele, type_zero_attrib, type_mat6, type_taylor, twiss_out, t
 logical, optional, intent(in) :: type_taylor, type_wake
 logical, optional, intent(in) :: type_zero_attrib
 logical, optional, intent(in) :: type_floor_coords, type_wall, type_rad_kick
-logical type_zero, err_flag, print_it, is_default, has_it, has_been_added, z1, z2
+logical type_zero, err_flag, print_it, is_default, has_it, has_been_added, is_zero1, is_zero2
 
 ! init
 
@@ -234,9 +234,10 @@ subroutine type_ele (ele, type_zero_attrib, type_mat6, type_taylor, twiss_out, t
  endif
  end select
 
- z1 = ((attrib%kind == is_real$ .or. attrib%kind == is_integer$) .and. attrib%value == 0)
- z2 = ((attrib2%kind == is_real$ .or. attrib2%kind == is_integer$) .and. attrib2%value == 0)
- if (z1 .and. z2 .and. .not. type_zero) cycle
+ is_zero1 = ((attrib%kind == is_real$ .or. attrib%kind == is_integer$) .and. attrib%value == 0)
+ is_zero2 = (((attrib2%kind == is_real$ .or. attrib2%kind == is_integer$) .and. attrib2%value == 0) .or. &
+ attrib2%name == null_name$)
+ if (is_zero1 .and. is_zero2 .and. .not. type_zero) cycle
 
  line = ''
  call write_this_attribute (attrib, ia, n_att, line(3:))
@@ -1622,7 +1623,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(98) = [character(40):: 'X_PITCH', 'Y_PITCH', 'X_OFFSET', &
+character(41), parameter :: att_name(100) = [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 +1637,10 @@ 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', 'Z_CROSSING', 'SPIN_TRACKING_MODEL']
+ 'F_FACTOR', 'EXACT_MULTIPOLES', 'Z_CROSSING', 'SPIN_TRACKING_MODEL', &
+ 'SPIN_DN_DPZ_X', 'INHERIT_FROM_FORK']
 
-character(41), parameter :: att2_name(98) = [character(40):: 'X_PITCH_TOT', 'Y_PITCH_TOT', 'X_OFFSET_TOT', &
+character(41), parameter :: att2_name(100) = [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 +1654,8 @@ 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', 'S_BETA_MIN', 'RECALC']
+ 'SCATTER_METHOD', 'FIDUCIAL_PT', 'S_BETA_MIN', 'RECALC', &
+ 'SPIN_DN_DPZ_Y', 'MODE_FLIP']
 
 ! Exceptional cases
 

bmad/parsing/bmad_parser_mod.f90

@@ -7320,21 +7320,13 @@ subroutine settable_dep_var_bookkeeping (ele)
  endif
 
 !------------------
+! Note: Dispersion will be handled by twiss_propagate1.
+
 case (beginning_ele$)
 
  if (ele%a%beta /= 0) ele%a%gamma = (1 + ele%a%alpha**2) / ele%a%beta
  if (ele%b%beta /= 0) ele%b%gamma = (1 + ele%b%alpha**2) / ele%b%beta
 
- ele%gamma_c = sqrt(1 - ele%c_mat(1,1)*ele%c_mat(2,2) + ele%c_mat(1,2)*ele%c_mat(2,1))
-
- call make_v_mats (ele, v_inv_mat = v_inv_mat)
- eta_vec = matmul (v_inv_mat, [ele%x%eta, ele%x%etap, ele%y%eta, ele%y%etap])
-
- ele%a%eta = eta_vec(1)
- ele%a%etap = eta_vec(2)
- ele%b%eta = eta_vec(3)
- ele%b%etap = eta_vec(4)
-
 !------------------
 ! Convert rbends to sbends and evaluate G if needed.
 ! Needed is the length and either: angle, G, or rho.