More SR wake longitudinal table devel. (#1288)

bmad/code/pointer_to_attribute.f90

@@ -533,26 +533,21 @@ subroutine pointer_to_attribute (ele, attrib_name, do_allocation, a_ptr, err_fla
 case ('N_LORD');          a_ptr%i => ele%n_lord
 case ('LR_FREQ_SPREAD', 'LR_SELF_WAKE_ON', 'LR_WAKE%AMP_SCALE', 'LR_WAKE%TIME_SCALE', &
       'LR_WAKE%FREQ_SPREAD', 'LR_WAKE%SELF_WAKE_ON', &
-      'SR_WAKE%SCALE_WITH_LENGTH', 'SR_WAKE%AMP_SCALE', 'SR_WAKE%Z_SCALE')
+      'SR_WAKE%SCALE_WITH_LENGTH', 'SR_WAKE%AMP_SCALE', 'SR_WAKE%Z_SCALE', &
+      'SR_WAKE%Z_LONG%SMOOTHING_SIGMA')
   if (.not. associated(ele%wake)) then
     if (.not. do_allocation) goto 9100
     call init_wake (ele%wake, 0, 0, 0, 0, .true.)
   endif
   select case (a_name)
-  case ('LR_SELF_WAKE_ON', 'LR_WAKE%SELF_WAKE_ON')
-    a_ptr%l => ele%wake%lr%self_wake_on
-  case ('LR_WAKE%AMP_SCALE')
-    a_ptr%r => ele%wake%lr%amp_scale
-  case ('LR_WAKE%TIME_SCALE')
-    a_ptr%r => ele%wake%lr%time_scale
-  case ('LR_FREQ_SPREAD', 'LR_WAKE%FREQ_SPREAD')
-    a_ptr%r => ele%wake%lr%freq_spread
-  case ('SR_WAKE%AMP_SCALE')
-    a_ptr%r => ele%wake%sr%amp_scale
-  case ('SR_WAKE%Z_SCALE')
-    a_ptr%r => ele%wake%sr%z_scale
-  case ('SR_WAKE%SCALE_WITH_LENGTH')
-    a_ptr%l => ele%wake%sr%scale_with_length
+  case ('SR_WAKE%Z_LONG%SMOOTHING_SIGMA');             a_ptr%r => ele%wake%sr%z_long%smoothing_sigma
+  case ('SR_WAKE%AMP_SCALE');                          a_ptr%r => ele%wake%sr%amp_scale
+  case ('SR_WAKE%Z_SCALE');                            a_ptr%r => ele%wake%sr%z_scale
+  case ('SR_WAKE%SCALE_WITH_LENGTH');                  a_ptr%l => ele%wake%sr%scale_with_length
+  case ('LR_SELF_WAKE_ON', 'LR_WAKE%SELF_WAKE_ON');    a_ptr%l => ele%wake%lr%self_wake_on
+  case ('LR_WAKE%AMP_SCALE');                          a_ptr%r => ele%wake%lr%amp_scale
+  case ('LR_WAKE%TIME_SCALE');                         a_ptr%r => ele%wake%lr%time_scale
+  case ('LR_FREQ_SPREAD', 'LR_WAKE%FREQ_SPREAD');      a_ptr%r => ele%wake%lr%freq_spread
   end select
 
 case ('H_MISALIGN%ACTIVE');     a_ptr%l => ele%photon%h_misalign%active

bmad/doc/attributes.tex

@@ -2930,19 +2930,23 @@ \subsection{Short-Range Wakes}
 for short-range wakefields are given in \Sref{s:sr.wake.eq}.  The \vn{sr_wake} attribute is used to set
 wakefield parameters. The general form of this attribute is:
 \begin{example}
-  sr_wake = \{z_max = <real>, z_scale = <real>, amp_scale = <real>, 
+  sr_wake = \{
+      z_max = <real>, z_scale = <real>, amp_scale = <real>, 
       scale_with_length = <logical>,
       longitudinal = \{<amp>, <damp>, <k>, <phi>, <position_dependence>\},
       ...
       longitudinal = \{...\}, 
       transverse = \{<amp>, <damp>, <k>, <phi>, <polarization>, <particle_dependence>\},
       ...
-      transverse = \{...\} \}
+      transverse = \{...\} 
+      z_long = \{...\}
+    \}
 \end{example}
 The \vn{sr_wake} structure has optional components \vn{z_max}, \vn{z_scale}, \vn{amp_scale}, and
-\vn{scale_with_length} along with zero or more \vn{longitudinal} sub-structures each one specifying
-a single longitudinal mode, and zero or more \vn{transverse} sub-structures each one specifying a
-single transverse mode. Example:
+\vn{scale_with_length}. To specify wakes, zero or more \vn{longitudinal} sub-structures each one specifying
+a single longitudinal mode, zero or more \vn{transverse} sub-structures each one specifying a
+single transverse mode and zero or one \vn{z_long} sub-structures that can be used to specify the longitudinal
+wake as a function of z-position (\sref{s:sr.wake.z.long}). Example:
 \begin{example}
   cav9: lcavity, ..., sr_wake = \{z_max = 1.3e-3, scale_with_length = F, 
        longitudinal = \{3.23e14, 1.23e3, 3.62e3, 0.123, none\}, 
@@ -2974,8 +2978,8 @@ \subsection{Short-Range Wakes}
 Monopole modes are modes that are independent of transverse position and dipole modes are modes that
 are linear in the transverse position. 
 
-For the \vn{longitudinal} sub-structures, there is a 5\Th
-component which gives the transverse position dependence of the wake. Possible values are:
+For the \vn{longitudinal} sub-structures, there is a 5\Th component, called
+\vn{position_dependence}, gives the transverse position dependence of the wake. Possible values are:
 \begin{example}
   none          ! No position dependence
   x_leading     ! Linear in the leading particle x-position
@@ -3026,6 +3030,49 @@ \subsection{Short-Range Wakes}
 zero and the transverse wake has no terms independent of the transverse offsets nor terms that
 depend upon the trailing particle offset.
 
+%-----------------------------------------------------------------
+\subsection{Longitudinal Short-Range Wake With a Wake table}
+\label{s:sr.wake.z.long}
+
+The \vn{z_long} substructure of \vn{sr_wake} (\sref{s:sr.wake.syntax}) is an alternative to the
+\vn{longitudinal} substructure. With the \vn{z_long} substructure, the wake is specified using a
+table of equally spaced points in $z$-position or as a function of time. The syntax for \vn{z_long}
+is
+\begin{example}
+  z_long = \{
+    time_based = <T/F>,
+    smoothing_sigma = <sigma>,
+    position_dependence = <dep>,
+    w = \{
+      <z1>  <w1>,            ! Wake <w> table as a function of z-position (or time).
+      <z2>  <w2>,
+      ...
+    \}
+  \}
+\end{example}
+The \vn{w} component gives the single particle wake as a function of either z-position or time
+depending upon the setting of the optional \vn{time_based} logical. The default is \vn{False}. The z
+(or time) values in the \vn{w} table must be equally spaced and in increasing order. Unlike the
+\vn{longitudinal} and \vn{transverse} pseudo-mode description, the wake can be finite for both
+positive and negative z (or time). When the wake is parsed, if \vn{time_based} is \vn{True}, the
+table is converted to be z-based using the conversion $z = -c \, t$ where $c$ is the speed of
+light. If the table is not symmetric around $z = 0$ (if there are more points on one size of zero
+than the other), the table is extended to be symmetric (this is needed for tracking).  During
+tracking, the half-width of the wake (the magnitude of the maximum or minimum z of the symmeterized
+table) must be greater than the full width of the bunch (excluding 1\% outlier particles). If not,
+an error message is generated and tracking is stopped. 
+
+The optional \vn{position_dependence} component determines if the wake depends upon the transverse coordinates
+of the bunch particles. See \sref{s:sr.wake.syntax} for details. The default is \vn{none}. 
+
+The wake is applied to a bunch by convoluting the bunch distribution with the wake.  The optional
+\vn{smoothing_sigma} component, if set non-zero, applies a Gaussian smoothing filter to the
+convolution. The value of \vn{smoothing_sigma} is in meters if \vn{time_based} is False or seconds
+otherwise.
+
+The \vn{z_scale} and \vn{amp_scale} parameters of the \vn{sr_wake} structure (\sref{s:sr.wake.syntax})
+are used with \vn{z_long}.
+
 %-----------------------------------------------------------------
 \subsection{Short-Range Wakes --- Old Format}
 \label{s:sr.wake.old}

bmad/doc/cover-page.tex

@@ -3,7 +3,7 @@
 
 \begin{flushright}
 \large
-  Revision: November 9, 2024 \\
+  Revision: November 14, 2024 \\
 \end{flushright}
 
 \pdfbookmark[0]{Preamble}{Preamble} 

bmad/modules/bmad_struct.f90

@@ -601,9 +601,10 @@ module bmad_struct
   complex(rp), allocatable :: fbunch(:), w_out(:)      ! Scratch space.
   real(rp) :: dz = 0                                   ! Distance between points.
   real(rp) :: z0 = 0                                   ! Wake extent is [-z0, z0].
-  integer :: plane = not_set$                          ! x$, y$, xy$, z$.
-  integer :: position_dependence = not_set$            ! Transverse: leading$, trailing$, none$
+  real(rp) :: smoothing_sigma = 0                      ! 0 => No smoothing.
+  integer :: position_dependence = none$               ! Transverse: leading$, trailing$, none$
                                                        ! Longitudinal: x_leading$, ..., y_trailing$, none$
+  logical :: time_based = .false.                      ! Was input time based?
 end type
 
 type wake_sr_mode_struct    ! Psudo-mode Short-range wake struct 

bmad/modules/equality_mod.f90

@@ -448,10 +448,12 @@ elemental function eq_wake_sr_z_long (f1, f2) result (is_eq)
 is_eq = is_eq .and. (f1%dz == f2%dz)
 !! f_side.equality_test[real, 0, NOT]
 is_eq = is_eq .and. (f1%z0 == f2%z0)
-!! f_side.equality_test[integer, 0, NOT]
-is_eq = is_eq .and. (f1%plane == f2%plane)
+!! f_side.equality_test[real, 0, NOT]
+is_eq = is_eq .and. (f1%smoothing_sigma == f2%smoothing_sigma)
 !! f_side.equality_test[integer, 0, NOT]
 is_eq = is_eq .and. (f1%position_dependence == f2%position_dependence)
+!! f_side.equality_test[logical, 0, NOT]
+is_eq = is_eq .and. (f1%time_based .eqv. f2%time_based)
 
 end function eq_wake_sr_z_long
 

bmad/multiparticle/wake_mod.f90

@@ -482,22 +482,21 @@ end subroutine sr_transverse_wake_particle
 !--------------------------------------------------------------------------
 !--------------------------------------------------------------------------
 !+
-! Subroutine sr_z_long_wake (ele, bunch, z_min, z_max)
+! Subroutine sr_z_long_wake (ele, bunch, z_ave)
 !
 ! Subroutine to apply the short-range z-wake kick to a particle.
 !
 ! Input:
 !   ele         -- ele_struct: Element with wake.
 !   bunch       -- bunch_struct: Bunch before wake applied.
-!   z_min       -- real(rp): Minimum z-position of all live particles.
-!   z_max       -- real(rp): Maximum z-position of all live particles.
+!   z_ave       -- real(rp): Average z-position of all live particles.
 !
 ! Output:
 !   orbit   -- coord_struct: Ending particle coords.
 !   bunch       -- bunch_struct: Bunch before wake applied.
 !+
 
-subroutine sr_z_long_wake (ele, bunch, z_min, z_max)
+subroutine sr_z_long_wake (ele, bunch, z_ave)
 
 use spline_mod
 
@@ -507,10 +506,10 @@ subroutine sr_z_long_wake (ele, bunch, z_min, z_max)
 type (wake_sr_z_long_struct), pointer :: srz
 type (coord_struct), pointer :: p, orbit
 
-real(rp) x, f0, ff, f_add, kick, dz
-real(rp) z_min, z_max, z_ave
+real(rp) x, f0, ff, f_add, kick, dz, rz_rel, r1, r2
+real(rp) z_ave
 
-integer i, j, ix, n2
+integer i, j, ix1, ix2, n2, n_bad, nn
 logical ok
 
 character(*), parameter :: r_name = 'sr_z_long_wake'
@@ -523,37 +522,51 @@ subroutine sr_z_long_wake (ele, bunch, z_min, z_max)
 srz => sr%z_long
 if (srz%dz == 0) return
 
-if (sr%z_scale * (z_max - z_min) > srz%z0) then
-  call out_io (s_error$, r_name, &
-      'The bunch is longer than the sr-z wake can handle for element: ' // ele%name)
-  p%state = lost$
-  return
-endif
-
 f0 = sr%amp_scale * bunch%charge_live
 if (sr%scale_with_length) f0 = f0 * ele%value(l$) 
 
-! Compute wake
+! Compute binned bunch distribution and wake
 
-z_ave = 0.5_rp * (z_min + z_max)
-n2 = size(srz%w_out) / 2
+nn = size(srz%w_out)
+n2 = (nn - 1) / 2
 srz%w_out = 0
+n_bad = 0
 
 do i = 1, size(bunch%particle)
   p => bunch%particle(i)
-  ix = nint(sr%z_scale * (p%vec(5) - z_ave) / srz%dz) + n2
   if (p%state /= alive$) cycle
 
-  select case (srz%plane)
+  rz_rel = sr%z_scale * (p%vec(5) - z_ave) / srz%dz + n2 + 1 
+  ix1 = floor(rz_rel)
+  ix2 = ceiling(rz_rel)
+  if (ix1 < 1 .or. ix2 > nn) then
+    n_bad = n_bad + 1
+    cycle
+  endif
+
+  r1 = ix2 - rz_rel
+  r2 = rz_rel - ix1
+
+  select case (srz%position_dependence)
   case (none$, x_trailing$, y_trailing$)
-    srz%w_out(ix) = srz%w_out(ix) + 1
+    srz%w_out(ix1) = srz%w_out(ix1) + r1
+    srz%w_out(ix2) = srz%w_out(ix2) + r2
   case (x_leading$)
-    srz%w_out(ix) = srz%w_out(ix) + p%vec(1)
+    srz%w_out(ix1) = srz%w_out(ix1) + r1 * p%vec(1)
+    srz%w_out(ix2) = srz%w_out(ix2) + r2 * p%vec(1)
   case (y_leading$)
-    srz%w_out(ix) = srz%w_out(ix) + p%vec(3)
+    srz%w_out(ix1) = srz%w_out(ix2) + r1 * p%vec(3)
+    srz%w_out(ix2) = srz%w_out(ix2) + r2 * p%vec(3)
   end select
 enddo
 
+if (n_bad > 0.01 * size(bunch%particle)) then
+  call out_io (s_error$, r_name, &
+      'The bunch is longer than the sr-z wake can handle for element: ' // ele%name)
+  p%state = lost$
+  return
+endif
+
 call fft_1d(srz%w_out, -1)
 srz%w_out = srz%w_out * srz%fw * f0
 call fft_1d(srz%w_out, 1)
@@ -563,36 +576,21 @@ subroutine sr_z_long_wake (ele, bunch, z_min, z_max)
 do i = 1, size(bunch%particle)
   p => bunch%particle(i)
   if (p%state /= alive$) cycle
-  ix = nint((p%vec(5) - z_ave) / srz%dz) + n2
 
-  select case (srz%plane)
-  case (none$, x_leading$, y_leading$)
-    p%vec(6) = p%vec(6) - srz%w_out(ix)
-  case (x_trailing$)
-    p%vec(6) = p%vec(6) - srz%w_out(ix) * p%vec(1)
-  case (y_trailing$)
-    p%vec(6) = p%vec(6) - srz%w_out(ix) * p%vec(3)
-  end select
-enddo
+  rz_rel = sr%z_scale * (p%vec(5) - z_ave) / srz%dz + n2 + 1 
+  ix1 = floor(rz_rel)
+  ix2 = ceiling(rz_rel)
 
-call fft_1d(srz%w_out, -1)
-srz%w_out = srz%w_out * srz%fw * f0
-call fft_1d(srz%w_out, 1)
-
-! Apply wake
-
-do i = 1, size(bunch%particle)
-  p => bunch%particle(i)
-  if (p%state /= alive$) cycle
-  ix = nint(sr%z_scale * (p%vec(5) - z_ave) / srz%dz) + n2
+  r1 = ix2 - rz_rel
+  r2 = rz_rel - ix1
 
-  select case (srz%plane)
+  select case (srz%position_dependence)
   case (none$, x_leading$, y_leading$)
-    p%vec(6) = p%vec(6) - srz%w_out(ix)
+    p%vec(6) = p%vec(6) - (r1 * srz%w_out(ix1) + r2 * srz%w_out(ix2))
   case (x_trailing$)
-    p%vec(6) = p%vec(6) - srz%w_out(ix) * p%vec(1)
+    p%vec(6) = p%vec(6) - (r1 * srz%w_out(ix1) + r2 * srz%w_out(ix2)) * p%vec(1)
   case (y_trailing$)
-    p%vec(6) = p%vec(6) - srz%w_out(ix) * p%vec(3)
+    p%vec(6) = p%vec(6) - (r1 * srz%w_out(ix1) + r2 * srz%w_out(ix2)) * p%vec(3)
   end select
 enddo
 
@@ -751,7 +749,7 @@ subroutine track1_sr_wake (bunch, ele)
 
 ! Z-wake
 
-call sr_z_long_wake(ele, bunch, p(i1)%vec(5), p(i2)%vec(5))
+call sr_z_long_wake(ele, bunch, p((i1+i2)/2)%vec(5))
 
 ! Loop over all particles in the bunch and apply the mode wakes
 

bmad/output/type_ele.f90

@@ -1422,14 +1422,10 @@ subroutine type_ele (ele, type_zero_attrib, type_mat6, type_taylor, twiss_out, t
       call re_allocate (li, nl+size(ele%wake%sr%z_long%w)+100, .false.)
       nl=nl+1; li(nl) = '  Short-Range Z-dependent Longitudinal wake:'
       srz => ele%wake%sr%z_long
-      if (srz%plane == z$) then
-        nl=nl+1; li(nl) = '  plane = ' // trim(sr_z_plane_name(srz%plane))
-      else
-        nl=nl+1; li(nl) = '  plane = ' // trim(sr_z_plane_name(srz%plane)) // &
-                          ', position_dependence = ' // trim(sr_transverse_position_dep_name(srz%position_dependence))
-      endif
-      nl=nl+1; li(nl) = '  dz = ' // to_str(srz%dz)
-      nl=nl+1; li(nl) = '  +/- Wake range: ' // to_str(srz%z0)
+      nl=nl+1; li(nl) = '    smoothing_sigma     = ' // to_str(srz%smoothing_sigma)
+      nl=nl+1; li(nl) = '    position_dependence = ' // trim(sr_transverse_position_dep_name(srz%position_dependence))
+      nl=nl+1; li(nl) = '    dz = ' // to_str(srz%dz)
+      nl=nl+1; li(nl) = '    +/- Wake range: ' // to_str(srz%z0)
       nl=nl+1; write (li(nl), '(a, i0)') '  # wake points: ', size(srz%w)
 
     else