New dispersion derivative. (#789)

* Start work on cooler element.
* Changed name of element from cooler to feedback.
* New dispersion derivative.

bmad/code/attribute_set_bookkeeping.f90

@@ -35,6 +35,18 @@ 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 ('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 default
   if (.not. field_attribute_free(ele, attrib_name)) then
     ele%field_master = .not. ele%field_master

bmad/code/pointer_to_attribute.f90

@@ -487,6 +487,11 @@ subroutine pointer_to_attribute (ele, attrib_name, do_allocation, a_ptr, err_fla
 case ('ETAP_X');          a_ptr%r => ele%x%etap
 case ('ETAP_Y');          a_ptr%r => ele%y%etap
 case ('ETAP_Z');          a_ptr%r => ele%z%etap
+case ('DETA_A_DS');       a_ptr%r => ele%a%deta_ds
+case ('DETA_B_DS');       a_ptr%r => ele%b%deta_ds
+case ('DETA_X_DS');       a_ptr%r => ele%x%deta_ds
+case ('DETA_Y_DS');       a_ptr%r => ele%y%deta_ds
+case ('DETA_Z_DS');       a_ptr%r => ele%z%deta_ds
 case ('CMAT_11');         a_ptr%r => ele%c_mat(1,1)
 case ('CMAT_12');         a_ptr%r => ele%c_mat(1,2)
 case ('CMAT_21');         a_ptr%r => ele%c_mat(2,1)

bmad/code/twiss_at_element.f90

@@ -140,19 +140,21 @@ subroutine zero_ave (ave)
 
 type (ele_struct) ave
 
-ave%s          = 0;   ave%s_start    = 0
-ave%a%phi      = 0;   ave%b%phi      = 0
-ave%a%alpha    = 0;   ave%b%alpha    = 0
-ave%a%beta     = 0;   ave%b%beta     = 0
-ave%a%gamma    = 0;   ave%b%gamma    = 0
-ave%a%eta      = 0;   ave%b%eta      = 0
-ave%a%etap     = 0;   ave%b%etap     = 0
-ave%a%sigma    = 0;   ave%b%sigma    = 0
-ave%x%eta      = 0;   ave%y%eta      = 0
-ave%x%etap     = 0;   ave%y%etap     = 0
-ave%c_mat      = 0
-ave%gamma_c    = 0                                            
-ave%value(l$)  = 0
+ave%s         = 0;   ave%s_start   = 0
+ave%a%phi     = 0;   ave%b%phi     = 0
+ave%a%alpha   = 0;   ave%b%alpha   = 0
+ave%a%beta    = 0;   ave%b%beta    = 0
+ave%a%gamma   = 0;   ave%b%gamma   = 0
+ave%a%eta     = 0;   ave%b%eta     = 0
+ave%a%etap    = 0;   ave%b%etap    = 0
+ave%a%deta_ds = 0;   ave%b%deta_ds = 0
+ave%a%sigma   = 0;   ave%b%sigma   = 0
+ave%x%eta     = 0;   ave%y%eta     = 0
+ave%x%etap    = 0;   ave%y%etap    = 0
+ave%x%deta_ds = 0;   ave%y%deta_ds = 0
+ave%c_mat     = 0
+ave%gamma_c   = 0                                            
+ave%value(l$) = 0
 
 end subroutine
 
@@ -166,29 +168,33 @@ subroutine twiss_ave (ave, ele1, r)
 
 !
 
-ave%s          = ave%s          + r * ele1%s          
-ave%s_start    = ave%s_start    + r * ele1%s_start
-ave%c_mat      = ave%c_mat      + r * ele1%c_mat      
-ave%gamma_c    = ave%gamma_c    + r * ele1%gamma_c    
-ave%a%phi      = ave%a%phi      + r * ele1%a%phi      
-ave%a%alpha    = ave%a%alpha    + r * ele1%a%alpha    
-ave%a%beta     = ave%a%beta     + r * ele1%a%beta     
-ave%a%gamma    = ave%a%gamma    + r * ele1%a%gamma    
-ave%a%eta      = ave%a%eta      + r * ele1%a%eta      
-ave%a%etap     = ave%a%etap     + r * ele1%a%etap     
-ave%a%sigma    = ave%a%sigma    + r * ele1%a%sigma    
-ave%x%eta      = ave%x%eta      + r * ele1%x%eta  
-ave%x%etap     = ave%x%etap     + r * ele1%x%etap 
-ave%b%phi      = ave%b%phi      + r * ele1%b%phi      
-ave%b%alpha    = ave%b%alpha    + r * ele1%b%alpha    
-ave%b%beta     = ave%b%beta     + r * ele1%b%beta     
-ave%b%gamma    = ave%b%gamma    + r * ele1%b%gamma    
-ave%b%eta      = ave%b%eta      + r * ele1%b%eta      
-ave%b%etap     = ave%b%etap     + r * ele1%b%etap     
-ave%b%sigma    = ave%b%sigma    + r * ele1%b%sigma    
-ave%y%eta      = ave%y%eta      + r * ele1%y%eta  
-ave%y%etap     = ave%y%etap     + r * ele1%y%etap 
-ave%value(l$)  = ave%value(l$)  + ele1%value(l$) 
+ave%s         = ave%s         + r * ele1%s          
+ave%s_start   = ave%s_start   + r * ele1%s_start
+ave%c_mat     = ave%c_mat     + r * ele1%c_mat      
+ave%gamma_c   = ave%gamma_c   + r * ele1%gamma_c    
+ave%a%phi     = ave%a%phi     + r * ele1%a%phi      
+ave%a%alpha   = ave%a%alpha   + r * ele1%a%alpha    
+ave%a%beta    = ave%a%beta    + r * ele1%a%beta     
+ave%a%gamma   = ave%a%gamma   + r * ele1%a%gamma    
+ave%a%eta     = ave%a%eta     + r * ele1%a%eta      
+ave%a%etap    = ave%a%etap    + r * ele1%a%etap     
+ave%a%deta_ds = ave%a%deta_ds + r * ele1%a%deta_ds     
+ave%a%sigma   = ave%a%sigma   + r * ele1%a%sigma    
+ave%x%eta     = ave%x%eta     + r * ele1%x%eta  
+ave%x%etap    = ave%x%etap    + r * ele1%x%etap 
+ave%x%deta_ds = ave%x%deta_ds + r * ele1%x%deta_ds 
+ave%b%phi     = ave%b%phi     + r * ele1%b%phi      
+ave%b%alpha   = ave%b%alpha   + r * ele1%b%alpha    
+ave%b%beta    = ave%b%beta    + r * ele1%b%beta     
+ave%b%gamma   = ave%b%gamma   + r * ele1%b%gamma    
+ave%b%eta     = ave%b%eta     + r * ele1%b%eta      
+ave%b%etap    = ave%b%etap    + r * ele1%b%etap     
+ave%b%deta_ds = ave%b%deta_ds + r * ele1%b%deta_ds     
+ave%b%sigma   = ave%b%sigma   + r * ele1%b%sigma    
+ave%y%eta     = ave%y%eta     + r * ele1%y%eta  
+ave%y%etap    = ave%y%etap    + r * ele1%y%etap 
+ave%y%deta_ds = ave%y%deta_ds + r * ele1%y%deta_ds 
+ave%value(l$) = ave%value(l$) + ele1%value(l$) 
 
 end subroutine
 

bmad/code/twiss_from_mat6.f90

@@ -114,17 +114,34 @@ subroutine twiss_from_mat6 (mat6, orb0, ele, stable, growth_rate, status, type_o
 vec(2) = m6(2,6) + (m6(2,2) * orb0(2) + m6(2,4) * orb0(4) - orb0(2))
 vec(3) = m6(3,6) + (m6(3,2) * orb0(2) + m6(3,4) * orb0(4))
 vec(4) = m6(4,6) + (m6(4,2) * orb0(2) + m6(4,4) * orb0(4) - orb0(4))
+
 eta_vec = matmul(mat4, vec)
 
-eta_vec(2) = eta_vec(2)
-eta_vec(4) = eta_vec(4)
+ele%x%eta     = eta_vec(1)
+ele%x%deta_ds = eta_vec(2)
+ele%y%eta     = eta_vec(3)
+ele%y%deta_ds = eta_vec(4)
+
+ele%z%eta     = 0
+ele%z%deta_ds = 1
+
+eta_vec = matmul(mat_symp_conj(v), 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)
+
+!
+
+eta_vec = matmul(mat4, mat6(1:4,6))
+
+ele%x%eta   = eta_vec(1)
+ele%x%etap  = eta_vec(2)
+ele%y%eta   = eta_vec(3)
+ele%y%etap  = eta_vec(4)
 
-ele%x%eta  = eta_vec(1)
-ele%x%etap = eta_vec(2)
-ele%y%eta  = eta_vec(3)
-ele%y%etap = eta_vec(4)
-ele%z%eta  = 0
-ele%z%etap = 1
+ele%z%eta     = 0
+ele%z%etap    = 1
 
 eta_vec = matmul(mat_symp_conj(v), eta_vec)
 ele%a%eta  = eta_vec(1)

bmad/code/twiss_propagate1.f90

@@ -193,7 +193,7 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag)
 ! p_z2 is p_z at end of ele2 assuming p_z = 1 at end of ele1.
 ! This is just 1.0 (except for RF cavities).
 
-eta1_vec = [ele1%x%eta, ele1%x%etap * rel_p1, ele1%y%eta, ele1%y%etap * rel_p1, ele1%z%eta, 1.0_rp]
+eta1_vec = [ele1%x%eta, ele1%x%deta_ds * rel_p1, ele1%y%eta, ele1%y%deta_ds * rel_p1, ele1%z%eta, 1.0_rp]
 
 ! For a circular ring, defining the dependence of z on pz is problematical.
 ! With the RF off, z is not periodic so dz/dpz is dependent upon turn number.
@@ -231,6 +231,26 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag)
 eta_vec(2) = eta_vec(2) / rel_p2
 eta_vec(4) = eta_vec(4) / rel_p2
 
+ele2%x%eta     = eta_vec(1)
+ele2%x%deta_ds = eta_vec(2)
+ele2%y%eta     = eta_vec(3)
+ele2%y%deta_ds = eta_vec(4)
+ele2%z%eta     = eta_vec(5)
+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))
+
+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%x%eta  = eta_vec(1)
 ele2%x%etap = eta_vec(2)
 ele2%y%eta  = eta_vec(3)
@@ -246,6 +266,8 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag)
 ele2%b%eta  = eta_vec(3)
 ele2%b%etap = eta_vec(4)
 
+!
+
 if (present(err_flag)) err_flag = .false.
 
 end subroutine

bmad/code/twiss_propagate_all.f90

@@ -56,13 +56,21 @@ subroutine twiss_propagate_all (lat, ix_branch, err_flag, ie_start, ie_end, zero
 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/code/type_twiss.f90

@@ -5,10 +5,10 @@
 ! If the lines argument is not present, the element information is printed to the terminal.
 !
 ! The compact form looks like:
-!           Beta     Alpha     Gamma       Phi        Eta       Etap
-!            (m)       (-)     (1/m)     (rad)        (m)        (-)
-!  X:    29.8929    -2.953     0.325   11.9116     1.4442     0.1347
-!  Y:     1.3982     0.015     0.715   11.6300    -0.0006     0.0033     
+!           Beta     Alpha     Gamma       Phi        Eta       Etap    dEta/ds
+!            (m)       (-)     (1/m)     (rad)        (m)        (-)        (-)
+!  X:    29.8929    -2.953     0.325   11.9116     1.4442     0.1347     0.1347
+!  Y:     1.3982     0.015     0.715   11.6300    -0.0006     0.0033     0.0033
 !
 ! The default verbose form looks like:
 !                          A              B            Cbar                        C_mat
@@ -18,6 +18,7 @@
 !  Phi (rad)        0.00000000     0.00000000            X              Y              Z
 !  Eta (m)         -0.00212550     0.00142557    -0.00212346     0.00144302     0.00000000
 !  Etap            -0.03500656    -0.00125557    -0.03513890    -0.00102335     1.00000000
+!  dEta/ds         -0.03500656    -0.00125557    -0.03513890    -0.00102335     1.00000000
 !  Sigma            0.00034547     0.00005437     0.00034547     0.00005437
 !
 ! Input:
@@ -88,11 +89,11 @@ subroutine type_twiss (ele, frequency_units, compact_format, lines, n_lines)
 if (logic_option (.false., compact_format)) then
   write (li(1), '(10x, a)')  &
             'Beta     Alpha     Gamma       Phi        Eta       Etap'
-  li(2) = trim(str) // '        (m)        (-)'
+  li(2) = trim(str) // '        (m)        (-)        (-)'
   write (li(3), fmt) ' X:', ele%a%beta,  &
-          ele%a%alpha, ele%a%gamma, coef*ele%a%phi, ele%a%eta, ele%a%etap
+          ele%a%alpha, ele%a%gamma, coef*ele%a%phi, ele%a%eta, ele%a%etap, ele%a%deta_ds
   write (li(4), fmt) ' Y:', ele%b%beta,  &
-          ele%b%alpha, ele%b%gamma, coef*ele%b%phi, ele%b%eta, ele%b%etap
+          ele%b%alpha, ele%b%gamma, coef*ele%b%phi, ele%b%eta, ele%b%etap, ele%b%deta_ds
   nl = 4
 
 else
@@ -104,7 +105,8 @@ subroutine type_twiss (ele, frequency_units, compact_format, lines, n_lines)
   write (li(5), '(2x, a12, 2a, 12x, a, 3(14x, a))') freq_str, v(ele%a%phi*coef), v(ele%b%phi*coef), 'X', 'Y', 'Z'
   write (li(6), '(2x, a12, 5a)')               'Eta (m)     ', v(ele%a%eta),  v(ele%b%eta),  v(ele%x%eta),  v(ele%y%eta),  v(ele%z%eta)
   write (li(7), '(2x, a12, 5a)')               'Etap        ', v(ele%a%etap), v(ele%b%etap), v(ele%x%etap), v(ele%y%etap), v(ele%z%etap)
-  nl = 7
+  write (li(8), '(2x, a12, 5a)')               'dEta/ds     ', v(ele%a%deta_ds), v(ele%b%deta_ds), v(ele%x%deta_ds), v(ele%y%deta_ds), v(ele%z%deta_ds)
+  nl = 8
   if (ele%a%sigma /= 0 .or. ele%b%sigma /= 0) then
     nl=nl+1; write (li(nl), '(2x, a12, 5a)')   'Sigma       ', v(ele%a%sigma), v(ele%b%sigma), v(ele%x%sigma), v(ele%y%sigma)
   endif

bmad/doc/cover-page.tex

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

bmad/doc/linear-optics.tex

@@ -382,21 +382,38 @@ \section{Dispersion Calculation}
 \label{s:dispersion}
 \index{dispersion|hyperbf}
 
-The dispersion ($\eta$) and the dispersion derivative ($\eta'$) are defined by the equations
-\begin{align}
-  \eta_x(s) &\equiv \left. \frac{dx}{dp_z} \right|_s \comma \qquad
-    \eta'_x(s) \equiv \left. \frac{d\eta_x}{ds} \right|_s
+The dispersion $\eta$ is defined in the standard way
+\begin{equation}
+  \text{eta_x} = \eta_x(s) \equiv \left. \frac{dx}{dp_z} \right|_s \comma \qquad
+  \text{eta_y} = \eta_y(s) \equiv \left. \frac{dy}{dp_z} \right|_s
+  \label{eedxdpz}
+\end{equation}
+
+There are two dispersion derivatives that Bmad uses. One dispersion derivative corresponds to
+the definition used by the MAD and SAD programs
+\begin{equation}
+  \text{etap_x} = \equiv \left. \frac{dp_x}{dp_z} \right|_s \comma \qquad 
+  \text{etap_y} = \equiv \left. \frac{dp_y}{dp_z} \right|_s \comma \qquad 
+\end{equation}
+This dispersion derivative is useful when constructing particle bunch distributions and for
+various calculations like for calculating radiation integrals.
+
+The one drawback with the above dispersion derivative definition is that it is not always simply
+related to the longitudinal derivative of the dispersion $d\eta/ds$. This becomes a factor when
+designing lattices where, if some section of the lattice needs to be dispersion free, it is
+convienient to be able to optimize $d\eta/ds$ to zero. For this reason, $d\eta/ds$ is calculated
+by \bmad
+\begin{equation}
+  \text{deta_x_ds} \equiv \left. \frac{d\eta_x}{ds} \right|_s
     = \left. \frac{dx'}{dp_z} \right|_s \CRNO
-  \eta_y(s) &\equiv \left. \frac{dy}{dp_z} \right|_s \comma \qquad
-    \eta'_y(s) \equiv \left. \frac{d\eta_y}{ds} \right|_s
-    = \left. \frac{dy'}{dp_z} \right|_s \\
-  \eta_z(s) &\equiv \left. \frac{dz}{dp_z} \right|_s \nonumber
-\end{align}
+  \text{deta_y_ds} \equiv \left. \frac{d\eta_y}{ds} \right|_s
+    = \left. \frac{dy'}{dp_z} \right|_s 
+\end{equation}
 
-Given the dispersion at a given point, the dispersion at some other point is calculated as follows:
-Let $\Bf r = (x, p_x, y, p_y, z, p_z)$ be the reference orbit, around which the dispersion is to be
-calculated. Let $\bfV$ and $\bf M$ be the zeroth and first order components of the transfer map
-between two points labeled 1 and 2:
+For $d\eta/ds$, given the dispersion at a given point, the dispersion at some other point is
+calculated as follows: Let $\Bf r = (x, p_x, y, p_y, z, p_z)$ be the reference orbit, around which
+the dispersion is to be calculated. Let $\bfV$ and $\bf M$ be the zeroth and first order components
+of the transfer map between two points labeled 1 and 2:
 \begin{equation}
   \Bf r_2 = \bfM \, \Bf r_1 + \bfV
   \label{rmrv}
@@ -405,7 +422,7 @@ \section{Dispersion Calculation}
 \begin{equation}
   \bfeta = 
   \left( 
-    \eta_x, \eta'_x \, (1 + p_z), \eta_y, \eta'_y \, (1 + p_z), \eta_z, 1
+    \eta_x, d\eta_x/ds \, (1 + p_z), d\eta_y/ds, \eta'_y \, (1 + p_z), \eta_z, 1
   \right)
 \end{equation}
 Differentiating \Eq{rmrv} with respect to energy, the dispersion at point 2 in terms of the
@@ -470,3 +487,5 @@ \section{Dispersion Calculation}
 local dispersion, on the other hand, reflects the correlations between the particle energy and
 particle position within a beam.
 
+To calulate the normal mode dispersions, \Eq{avx} is used to transform from laboratory to normal mode
+coordinates.