Merge branch 'main' into OpenPMD/01

.gitignore

@@ -103,12 +103,12 @@ regression_tests/parse_test/z.bmad
 regression_tests/parse_test/z2.bmad
 regression_tests/regression.results
 regression_tests/sad_test/sler_1689.bmad
-regression_tests/twiss_track_test/digested.file
-regression_tests/twiss_track_test/rad_int_no_wig_cache.dat
-regression_tests/twiss_track_test/rad_int_no_wig_no_cache.dat
-regression_tests/twiss_track_test/rad_int_wig_cache.dat
-regression_tests/twiss_track_test/rad_int_wig_no_cache.dat
-regression_tests/twiss_track_test/twiss.out
+regression_tests/cesr_test/digested.file
+regression_tests/cesr_test/rad_int_no_wig_cache.dat
+regression_tests/cesr_test/rad_int_no_wig_no_cache.dat
+regression_tests/cesr_test/rad_int_wig_cache.dat
+regression_tests/cesr_test/rad_int_wig_no_cache.dat
+regression_tests/cesr_test/twiss.out
 regression_tests/long_term_tracking_test/sim1.dat
 regression_tests/long_term_tracking_test/beam-10.h5
 regression_tests/write_bmad_test/em_field1.gen_grad_map

bmad/code/attribute_set_bookkeeping.f90

@@ -35,17 +35,17 @@ subroutine attribute_set_bookkeeping (ele, attrib_name, err_flag, attrib_ptr)
 case ('GRADIENT')
     ele%value(voltage$) = ele%value(gradient$) * ele%value(l$)
 
-case ('ETAP_A'); ele%a%deta_ds = ele%a%etap
-case ('ETAP_B'); ele%b%deta_ds = ele%b%etap
-case ('ETAP_X'); ele%x%deta_ds = ele%x%etap
-case ('ETAP_Y'); ele%y%deta_ds = ele%y%etap
-case ('ETAP_Z'); ele%z%deta_ds = ele%z%etap
+case ('ETAP_A'); ele%a%deta_ds = real_garbage$
+case ('ETAP_B'); ele%b%deta_ds = real_garbage$
+case ('ETAP_X'); ele%x%deta_ds = real_garbage$
+case ('ETAP_Y'); ele%y%deta_ds = real_garbage$
+case ('ETAP_Z'); ele%z%deta_ds = real_garbage$
 
-case ('DETA_A_DS'); ele%a%etap = ele%a%deta_ds
-case ('DETA_B_DS'); ele%b%etap = ele%b%deta_ds
-case ('DETA_X_DS'); ele%x%etap = ele%x%deta_ds
-case ('DETA_Y_DS'); ele%y%etap = ele%y%deta_ds
-case ('DETA_Z_DS'); ele%z%etap = ele%z%deta_ds
+case ('DETA_A_DS'); ele%a%etap = real_garbage$
+case ('DETA_B_DS'); ele%b%etap = real_garbage$
+case ('DETA_X_DS'); ele%x%etap = real_garbage$
+case ('DETA_Y_DS'); ele%y%etap = real_garbage$
+case ('DETA_Z_DS'); ele%z%etap = real_garbage$
 
 case default
   if (.not. field_attribute_free(ele, attrib_name)) then

bmad/code/set_twiss.f90

@@ -0,0 +1,81 @@
+!+
+! Subroutine set_twiss(branch, twiss_ele, ix_ele, err_flag, print_err)
+!
+! Routine to set the beginning Twiss of a lattice branch so that the Twiss parameters at branch%ele(ix_ele) are
+! the same as the Twiss parameters in twiss_ele.
+!
+! It is OK for twiss_ele to be an element in the branch (this routine makes a copy to be safe).
+!
+! Input:
+!   branch      -- branch_struct: Branch to modify.
+!   twiss_ele   -- ele_struct: Element with desired Twiss parameters.
+!   ix_ele      -- integer: Match branch%ele(ix_ele) Twiss to twiss_ele.
+!   err_flag    -- logical: Set True if there is an error. False otherwise.
+!   print_err   -- logical, optional: Print an error message if there is an error? Default is True.
+!-
+
+subroutine set_twiss(branch, twiss_ele, ix_ele, err_flag, print_err)
+
+use bmad, dummy => set_twiss
+
+implicit none
+
+type (branch_struct), target :: branch
+type (ele_struct) twiss_ele, t_ele, ele0, ele1
+type (ele_struct), pointer :: b_ele
+
+real(rp) xmat(6,6), xvec(6)
+
+integer ix_ele
+integer ie
+
+logical err_flag
+logical, optional :: print_err
+
+character(*), parameter :: r_name = 'set_twiss'
+
+!
+
+err_flag = .true.
+
+if (branch%param%geometry == closed$) then
+  if (logic_option(.true., print_err)) call out_io (s_error$, r_name, &
+                                             'TWISS CANNOT BE SET IN A BRANCH WITH AN OPEN GEOMETRY.')
+  return
+endif
+
+if (ix_ele < 0 .or. ix_ele > branch%n_ele_track) then
+  if (logic_option(.true., print_err)) then
+    call out_io (s_error$, r_name, 'IX_ELE ELEMENT INDEX (' // &
+                                int_str(ix_ele) // ') OUT OF RANGE: [0, ' // int_str(branch%n_ele_track) // '].')
+  endif
+  return
+endif
+
+!
+
+if (ix_ele == 0) then
+  call transfer_twiss(twiss_ele, branch%ele(0))
+  call twiss_propagate_all(branch%lat, branch%ix_branch, err_flag)
+  return
+endif
+
+!
+
+call transfer_twiss(twiss_ele, ele0)
+
+call transfer_matrix_calc(branch%lat, ele1%mat6, ele1%vec0, 0, ix_ele, branch%ix_branch)
+call mat_inverse(ele1%mat6, ele1%mat6)
+ele1%key = quadrupole$
+ele1%map_ref_orb_in = branch%ele(ix_ele)%map_ref_orb_out
+ele1%map_ref_orb_out = branch%ele(1)%map_ref_orb_in
+call twiss_propagate1(ele0, ele1, err_flag)
+if (err_flag) then
+  if (logic_option(.true., print_err)) call out_io (s_error$, r_name, 'CANNOT BACK PROPAGATE TWISS.')
+  return
+endif
+
+call transfer_twiss(ele1, branch%ele(0))
+call twiss_propagate_all(branch%lat, branch%ix_branch, err_flag)
+
+end subroutine set_twiss

bmad/code/twiss_propagate1.f90

@@ -32,14 +32,29 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag)
 real(rp) v_mat(4,4), v_inv_mat(4,4), 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), dpz2_dpz1, rel_p1, rel_p2
+real(rp) mat2(2,2), eta1_vec(6), eta_vec(6), vec(6), dpz2_dpz1, rel_p1, rel_p2
 real(rp) det_factor, deriv_rel, gamma2_c, df
 
 logical error
 logical, optional :: err_flag
 
 character(*), parameter :: r_name = 'twiss_propagate1'
 
+! Beginning element bookkeeping
+
+if (ele1%key == beginning_ele$) then
+  ele1%map_ref_orb_out = ele2%map_ref_orb_in
+  rel_p1 = 1 + ele1%map_ref_orb_out%vec(6)
+
+  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
+    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
+endif
+
 !
 
 if (present(err_flag)) err_flag = .true.
@@ -219,52 +234,36 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag)
   else
     dpz2_dpz1 = 1 / rel_p2
   endif
-  deriv_rel = dpz2_dpz1 * rel_p1
-  eta_vec(1) = (mat6(1,2) * orb%vec(2) + mat6(1,4) * orb%vec(4)) / deriv_rel
-  eta_vec(2) = (mat6(2,2) * orb%vec(2) + mat6(2,4) * orb%vec(4)) / deriv_rel - orb_out%vec(2) / rel_p2
-  eta_vec(3) = (mat6(3,2) * orb%vec(2) + mat6(3,4) * orb%vec(4)) / deriv_rel
-  eta_vec(4) = (mat6(4,2) * orb%vec(2) + mat6(4,4) * orb%vec(4)) / deriv_rel - orb_out%vec(4) / rel_p2
-  eta_vec(5) = (mat6(5,2) * orb%vec(2) + mat6(5,4) * orb%vec(4)) / deriv_rel
-  eta_vec(1:5) = eta_vec(1:5) + matmul (mat6(1:5,:), eta1_vec) / dpz2_dpz1
 endif
 
-eta_vec(2) = eta_vec(2) / rel_p2
-eta_vec(4) = eta_vec(4) / rel_p2
+!
+
+eta1_vec = [ele1%x%eta, ele1%x%etap, ele1%y%eta, ele1%y%etap, ele1%z%eta, 1.0_rp]
+eta_vec(1:5) = matmul (mat6(1:5,:), eta1_vec) / dpz2_dpz1
 
 ele2%x%eta     = eta_vec(1)
-ele2%x%deta_ds = eta_vec(2)
+ele2%x%etap    = eta_vec(2)
+ele2%x%deta_ds = eta_vec(2) / rel_p2 - orb_out%vec(2) / rel_p2**2
+
 ele2%y%eta     = eta_vec(3)
-ele2%y%deta_ds = eta_vec(4)
+ele2%y%etap    = eta_vec(4)
+ele2%y%deta_ds = eta_vec(4) / rel_p2 - orb_out%vec(4) / rel_p2**2
+
 ele2%z%eta     = eta_vec(5)
+ele2%z%etap    = ele1%z%etap * dpz2_dpz1
 ele2%z%deta_ds = ele1%z%deta_ds * dpz2_dpz1
 
 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%deta_ds = eta_vec(2)
-ele2%b%eta     = eta_vec(3)
-ele2%b%deta_ds = eta_vec(4)
-
-!
-
-eta1_vec = [ele1%x%eta, ele1%x%etap, ele1%y%eta, ele1%y%etap, ele1%z%eta, 1.0_rp]
-eta_vec(1:5) = matmul (mat6(1:5,:), eta1_vec) / dpz2_dpz1
+ele2%a%etap    = eta_vec(2)
+ele2%a%deta_ds = eta_vec(2) / rel_p2 - vec(2) / rel_p2**2
 
-ele2%x%eta  = eta_vec(1)
-ele2%x%etap = eta_vec(2)
-ele2%y%eta  = eta_vec(3)
-ele2%y%etap = eta_vec(4)
-ele2%z%eta  = eta_vec(5)
-ele2%z%etap = ele1%z%etap * dpz2_dpz1
-
-call make_v_mats (ele2, v_mat, v_inv_mat)
-eta_vec(1:4) = matmul (v_inv_mat, eta_vec(1:4))
-
-ele2%a%eta  = eta_vec(1)
-ele2%a%etap = eta_vec(2)
-ele2%b%eta  = eta_vec(3)
-ele2%b%etap = eta_vec(4)
+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
 
 !
 
@@ -326,7 +325,7 @@ subroutine twiss1_propagate (twiss1, mat2, ele_key, length, twiss2, err)
 
 det = determinant (mat2)
 
-if (det == 0 .or. twiss1%beta == 0) return
+if (det == 0 .or. twiss1%beta == 0 .or. twiss1%beta > 1d100) return  ! Limit max beta to prevent numerical overflow.
 
 m11 = mat2(1,1)
 m12 = mat2(1,2)

bmad/code/twiss_propagate_all.f90

@@ -54,22 +54,6 @@ subroutine twiss_propagate_all (lat, ix_branch, err_flag, ie_start, ie_end)
 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
 
-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)
-ele%a%eta  = eta_vec(1)
-ele%a%etap = eta_vec(2)
-ele%b%eta  = eta_vec(3)
-ele%b%etap = eta_vec(4)
-
-eta_vec = [ele%x%eta, ele%x%deta_ds, ele%y%eta, ele%y%deta_ds]
-eta_vec = matmul (v_inv_mat, eta_vec)
-ele%a%eta     = eta_vec(1)
-ele%a%deta_ds = eta_vec(2)
-ele%b%eta     = eta_vec(3)
-ele%b%deta_ds = eta_vec(4)
-
 ! Propagate twiss
 
 if (present(err_flag)) err_flag = .true.

bmad/doc/cover-page.tex

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

bmad/doc/ele-struct.tex

@@ -187,7 +187,7 @@ \section{Element Attribute Bookkeeping}
 programs to be slow if attributes are continually being changed. In order to keep track what
 bookkeeping has been done, the \vn{ele%status} component is used by the appropriate bookkeeping
 routines for making sure the bookkeeping overhead is keep to a minimum. However, ``intelligent''
-bookkeeping is only done if explicitly enabled in a program. See \sref{s:lat.bookkeeping} for more
+bookkeeping is only done if explicitly enabled in a program. See \sref{s:lat.book} for more
 details.
 
 %--------------------------------------------------------------------------

bmad/doc/elements.tex

@@ -2154,8 +2154,9 @@ \section{Foil}
   \caption[Foil geometry.]
 {Foil geometry. Parameters \vn{thickness}, \vn{dthickness_dx}, \vn{x1_edge}, \vn{x2_edge},
 \vn{y1_edge}, and \vn{y2_edge} determine the foil geometry in body coordinates (see \fig{f:coords}).
-Orientation parameters like \vn{x_offset} (\sref{s:offset}) orient the foil with respect to
-laboratory coordinates but do not change the foil shape.}
+For the drawing, \vn{x1_edge} is negative (the associated edge is to the left of $x = 0$) and the
+\vn{x2_edge} is positive.  Orientation parameters like \vn{x_offset} (\sref{s:offset}) orient the
+foil with respect to laboratory coordinates but do not change the foil shape.}
   \label{f:foil}
 \end{figure}
 
@@ -2194,10 +2195,10 @@ \section{Foil}
   scatter_test        = <Logic>   ! For testing scattering. Default: False.
   scatter_method      = <Switch>  ! Scattering algorithm. Default: highland.
   dthickness_dx       = <Real>    ! Wedge slope when the foil is wedge shaped.
-  x1_edge             = <Real>    ! Foil edge in the +x direction. Default: -99 m.
-  x2_edge             = <Real>    ! Foil edge in the -x direction. Default:  99 m.
-  y1_edge             = <Real>    ! Foil edge in the +y direction. Default: -99 m.
-  y2_edge             = <Real>    ! Foil edge in the -y direction. Default:  99 m.
+  x1_edge             = <Real>    ! Foil edge in the x-direction. Default: -99 m.
+  x2_edge             = <Real>    ! Foil edge in the x-direction. Default:  99 m.
+  y1_edge             = <Real>    ! Foil edge in the y-direction. Default: -99 m.
+  y2_edge             = <Real>    ! Foil edge in the y-direction. Default:  99 m.
 \end{example}
 
 Scattering is simulated to be Gaussian distributed with a sigma calculated in one of two methods.
@@ -3388,7 +3389,7 @@ \section{Match}
   alpha_a1, alpha_b1 = <Real>            ! Exit alphas
   eta_x0, eta_y0     = <Real>            ! Entrance etas 
   eta_x1, eta_y1     = <Real>            ! Exit etas 
-  etap_x0, etap_y0   = <Real>            ! Entrance eta' 
+  etap_x0, etap_y0   = <Real>            ! Entrance momentum dispersion
   etap_x1, etap_y1   = <Real>            ! Exit eta'
   c11_mat0, c12_mat0, c21_mat0, c22_mat0 = <Real>  ! Entrance coupling.
   c11_mat1, c12_mat1, c21_mat1, c22_mat1 = <Real>  ! Exit coupling.

bmad/doc/example-program.tex

@@ -100,7 +100,7 @@ \section{Explanation of the Simple_Bmad_Program}
 subroutines. In particular, the \vn{lat} variable (line 5), which is of type \vn{lat_struct}
 (\sref{s:lat.struct}), holds all of the lattice information: The list of elements, their attributes,
 etc.  The manditory setting of \vn{bmad_com%auto_bookkeeper} to \vn{False} in line 13 enables the
-``intelligent'' bookkeeping of lattice attributes as discussed in \sref{s:lat.bookkeeping}).  The
+``intelligent'' bookkeeping of lattice attributes as discussed in \sref{s:lat.book}).  The
 call to \Hyperref{r:bmad.parser}{bmad_parser} (line 17) causes the lattice file \vn{lat.bmad} to be
 parsed and the lattice information is stored the \vn{lat} variable. Note: To get a listing of the
 \vn{lat_struct} components or to find out more about \vn{bmad_parser} use the \vn{getf} command as
@@ -118,7 +118,7 @@ \section{Explanation of the Simple_Bmad_Program}
 changed, the call to \Hyperref{r:set.flags.for.changed.attribute}{set_flags_for_changed_attribute}
 in line 23 is needed for \bmad to inform \bmad that this attribute has changed and the call to
 \Hyperref{r:lattice.bookkeeper}{lattice_bookkeeper} does the necessary lattice bookkeeping
-(\sref{s:lat.bookkeeping}).
+(\sref{s:lat.book}).
 
 The call to \Hyperref{r:lat.make.mat6}{lat_make_mat6} in line 25 recalculates the linear transfer
 matrix for the \vn{CLEO_SOL} element.