From eb0cd7e04ca2ce952c8ae563bc182acf94a57aad Mon Sep 17 00:00:00 2001 From: David Sagan Date: Mon, 3 Jun 2024 17:50:17 -0500 Subject: [PATCH] More set_twiss and dispersion devel. (#990) * More set_twiss devel. * More dispersion devel. --- bmad/code/reverse_lat.f90 | 6 +- bmad/code/set_twiss.f90 | 40 ++++-- bmad/code/twiss_at_start.f90 | 70 ++-------- bmad/code/twiss_from_mat6.f90 | 49 +++---- bmad/code/twiss_propagate1.f90 | 54 ++++---- bmad/doc/linear-optics.tex | 45 +++++-- bmad/low_level/transfer_twiss.f90 | 2 + bmad/modules/bmad_routine_interface.f90 | 4 +- bmad/modules/bmad_struct.f90 | 1 - bmad/parsing/set_ele_defaults.f90 | 2 + regression_tests/cesr_test/output.correct | 138 ++++++++++---------- regression_tests/twiss_test/twiss_test.bmad | 100 ++++++++++---- regression_tests/twiss_test/twiss_test.f90 | 134 +++++++++++-------- tao/version/tao_version_mod.f90 | 2 +- 14 files changed, 350 insertions(+), 297 deletions(-) diff --git a/bmad/code/reverse_lat.f90 b/bmad/code/reverse_lat.f90 index feb533173d..8a0c5bcedb 100644 --- a/bmad/code/reverse_lat.f90 +++ b/bmad/code/reverse_lat.f90 @@ -71,11 +71,7 @@ subroutine reverse_lat (lat_in, lat_rev, track_antiparticle) ele_rev%s = 0 ele_rev%floor = branch_in%ele(nt)%floor - call transfer_twiss (branch_in%ele(nt), ele_rev) - ele_rev%a%alpha = -ele_rev%a%alpha - ele_rev%b%alpha = -ele_rev%b%alpha - ele_rev%c_mat(1,2) = -ele_rev%c_mat(1,2) - ele_rev%c_mat(2,1) = -ele_rev%c_mat(2,1) + call transfer_twiss (branch_in%ele(nt), ele_rev, .true.) ! All other elements diff --git a/bmad/code/set_twiss.f90 b/bmad/code/set_twiss.f90 index 54a65c259f..84f8f324f1 100644 --- a/bmad/code/set_twiss.f90 +++ b/bmad/code/set_twiss.f90 @@ -1,20 +1,27 @@ !+ -! Subroutine set_twiss(branch, twiss_ele, ix_ele, err_flag, print_err) +! Subroutine set_twiss(branch, twiss_ele, ix_ele, match_deta_ds, 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. +! For an open lattice branch, set the beginning Twiss, coupling, and dispersion of the branch so that the +! Twiss parameters at the "target" element 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). ! +! Since the the reference orbit at twiss_ele may be different from the reference orbit at the target ele, +! it is not possible, in general, to match both etap and deta_ds. Which is matched is set by match_deta_ds. +! +! Additionally, It is not possible to simultaneously match both (x,y) deta_ds and (a,b) normal mode deta_ds. +! This routine matches (x,y) deta_ds only (when match_deta_ds = True). +! ! 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. +! 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. +! match_deta_ds -- logical: If True, match deta_ds. If False, match etap. +! 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) +subroutine set_twiss(branch, twiss_ele, ix_ele, match_deta_ds, err_flag, print_err) use bmad, dummy => set_twiss @@ -24,12 +31,12 @@ subroutine set_twiss(branch, twiss_ele, ix_ele, err_flag, print_err) type (ele_struct) twiss_ele, ele0, ele1 type (ele_struct), pointer :: ele -real(rp) xmat(6,6), xvec(6) +real(rp) xmat(6,6), xvec(6), rel_p integer ix_ele integer ie -logical err_flag +logical match_deta_ds, err_flag logical, optional :: print_err character(*), parameter :: r_name = 'set_twiss' @@ -52,9 +59,18 @@ subroutine set_twiss(branch, twiss_ele, ix_ele, err_flag, print_err) return endif -! +! Note: twiss_propagate1 treats deta_ds as a dependent variable to etap call transfer_twiss(twiss_ele, ele0) +ele0%map_ref_orb_out = twiss_ele%map_ref_orb_out + +if (match_deta_ds) then + ele => branch%ele(ix_ele) + rel_p = 1.0_rp + ele%map_ref_orb_out%vec(6) + ele0%x%etap = ele0%x%deta_ds * rel_p + ele%map_ref_orb_out%vec(2) / rel_p + ele0%y%etap = ele0%y%deta_ds * rel_p + ele%map_ref_orb_out%vec(4) / rel_p +endif + do ie = ix_ele, 1, -1 ele => branch%ele(ie) call mat_inverse(ele%mat6, ele1%mat6) diff --git a/bmad/code/twiss_at_start.f90 b/bmad/code/twiss_at_start.f90 index ea3e3717f5..52f112fc0a 100644 --- a/bmad/code/twiss_at_start.f90 +++ b/bmad/code/twiss_at_start.f90 @@ -37,7 +37,7 @@ subroutine twiss_at_start (lat, status, ix_branch, type_out) type (ele_struct), target :: drift_ele type (branch_struct), pointer :: branch -real(rp) eta_vec(4), t0_4(4,4), mat6(6,6), map0(4), m56 +real(rp) eta_vec(5), t0_5(5,5), mat6(6,6), map0(5), m56 integer, optional, intent(in) :: ix_branch integer, optional, intent(out) ::status @@ -51,7 +51,7 @@ subroutine twiss_at_start (lat, status, ix_branch, type_out) ! Init one turn. T0 is the transverse part of the matrix -call mat_make_unit (t0_4) ! form unit matrix +call mat_make_unit (t0_5) ! form unit matrix eta_vec = 0 map0 = 0 m56 = 0 @@ -61,7 +61,7 @@ subroutine twiss_at_start (lat, status, ix_branch, type_out) ! Propagate the transfer map around branch. ! Since the RF is taken to be off we use a trick so we only have to multiply -! 4x4 matrices. +! 5x5 matrices. if (debug) then iu = lunget() @@ -80,23 +80,22 @@ subroutine twiss_at_start (lat, status, ix_branch, type_out) ele => drift_ele endif - m56 = m56 + ele%mat6(5,6) + dot_product(ele%mat6(5,1:4), eta_vec) - eta_vec = matmul (ele%mat6(1:4,1:4), eta_vec) + ele%mat6(1:4,6) - map0 = matmul (ele%mat6(1:4,1:4), map0) + ele%vec0(1:4) - t0_4 = matmul (ele%mat6(1:4,1:4), t0_4) + eta_vec = matmul (ele%mat6(1:5,1:5), eta_vec) + ele%mat6(1:5,6) + map0 = matmul (ele%mat6(1:5,1:5), map0) + ele%vec0(1:5) + t0_5 = matmul (ele%mat6(1:5,1:5), t0_5) if (debug) then write (iu, *) '!------------------------------------', n call type_ele (ele, .false., 0, .false., 0, lines = lines, n_lines = n_lines) do i = 1, n_lines write (iu, '(a)') trim(lines(i)) enddo - write (iu, *) 'Symplectic Check:', mat_symp_error(t0_4) - do i = 1, 4 - write (iu, '(5x, 4f14.8, 6x, 4f14.8)') t0_4(i,:), ele%mat6(i,1:4) + write (iu, *) 'Symplectic Check:', mat_symp_error(t0_5(1:4,1:4)) + do i = 1, 5 + write (iu, '(5x, 5f14.8, 6x, 5f14.8)') t0_5(i,:), ele%mat6(i,1:5) enddo write (iu, '(a, 6es18.10)') 'Map_ref_orb:', ele%map_ref_orb_in%vec - write (iu, '(a, 4es20.12)') 'Eta: ', eta_vec - write (iu, '(a, 4es20.12)') 'Map0:', map0 + write (iu, '(a, 5es20.12)') 'Eta: ', eta_vec + write (iu, '(a, 5es20.12)') 'Map0:', map0 endif enddo @@ -105,10 +104,8 @@ subroutine twiss_at_start (lat, status, ix_branch, type_out) ! Put 1-turn matrix into branch%param%t1_no_RF call mat_make_unit (mat6) -mat6(1:4,1:4) = t0_4 - -call mat6_dispersion (eta_vec, mat6) ! dispersion to %mat6 -mat6(5,6) = m56 +mat6(1:5,1:5) = t0_5 +mat6(1:5, 6) = eta_vec branch%param%t1_no_RF = mat6 ! Compute twiss parameters @@ -125,45 +122,4 @@ subroutine twiss_at_start (lat, status, ix_branch, type_out) branch%ele(0)%a%phi = 0 branch%ele(0)%b%phi = 0 -contains - -!-------------------------------------------------------------------------- -!-------------------------------------------------------------------------- -!-------------------------------------------------------------------------- -!+ -! Subroutine mat6_dispersion (m_i6, mat6) -! -! Subroutine to set the mat6(5, 1:4) terms given the vector mat6(1:4, 6) -! which is a measure of the dispersion. -! -! Input: -! m_i6(4) -- Real(rp): mat6(1:4, 6) components. -! mat6(6,6) -- Real(rp): Matrix with 4x4 x-y submatrix already made. -! -! Output: -! mat6(6,6) -- Real(rp): mat6(5, 1:4) components set. -!- - -subroutine mat6_dispersion (m_i6, mat6) - -implicit none - -real(rp), intent(inout) :: mat6(:,:) -real(rp), intent(in) :: m_i6(:) - -real(rp) vec4(4) - -! - -mat6(1:4, 6) = m_i6(1:4) - -vec4(1) = -m_i6(2) -vec4(2) = m_i6(1) -vec4(3) = -m_i6(4) -vec4(4) = m_i6(3) - -mat6(5,1:4) = matmul (vec4, mat6(1:4,1:4)) - -end subroutine mat6_dispersion - end subroutine diff --git a/bmad/code/twiss_from_mat6.f90 b/bmad/code/twiss_from_mat6.f90 index ae29a2e037..aaf06f5bff 100644 --- a/bmad/code/twiss_from_mat6.f90 +++ b/bmad/code/twiss_from_mat6.f90 @@ -47,12 +47,7 @@ subroutine twiss_from_mat6 (mat6, orb0, ele, stable, growth_rate, status, type_o stable = .false. m6 = mat6 - -if (bmad_private%normalize_twiss) then - rel_p = 1 + orb0(6) - m6(1:5:2, 2:6:2) = m6(1:5:2, 2:6:2) * rel_p - m6(2:6:2, 1:5:2) = m6(2:6:2, 1:5:2) / rel_p -endif +rel_p = 1 + orb0(6) mat4 = m6(1:4, 1:4) @@ -110,43 +105,31 @@ subroutine twiss_from_mat6 (mat6, orb0, ele, stable, growth_rate, status, type_o forall (i = 1:4) mat4(i,i) = mat4(i,i) + 1 call mat_inverse (mat4, mat4) -vec(1) = m6(1,6) + (m6(1,2) * orb0(2) + m6(1,4) * orb0(4)) -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 = matmul(mat4, mat6(1:4,6)) ele%x%eta = eta_vec(1) -ele%x%deta_ds = eta_vec(2) +ele%x%etap = eta_vec(2) +ele%x%deta_ds = ele%x%etap / rel_p - orb0(2) / rel_p**2 + ele%y%eta = eta_vec(3) -ele%y%deta_ds = eta_vec(4) +ele%y%etap = eta_vec(4) +ele%y%deta_ds = eta_vec(4) / rel_p - orb0(4) / rel_p*2 ele%z%eta = 0 +ele%z%etap = 1 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) +eta_vec = matmul(mat_symp_conj(v), eta_vec) +vec = matmul(mat_symp_conj(v), orb0(1:4)) -ele%z%eta = 0 -ele%z%etap = 1 +ele%a%eta = eta_vec(1) +ele%a%etap = eta_vec(2) +ele%a%deta_ds = ele%a%etap / rel_p - vec(2) / rel_p**2 -eta_vec = matmul(mat_symp_conj(v), 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) +ele%b%eta = eta_vec(3) +ele%b%etap = eta_vec(4) +ele%b%deta_ds = eta_vec(4) / rel_p - vec(4) / rel_p*2 end subroutine diff --git a/bmad/code/twiss_propagate1.f90 b/bmad/code/twiss_propagate1.f90 index 403fc3fbf6..e44456e8cd 100644 --- a/bmad/code/twiss_propagate1.f90 +++ b/bmad/code/twiss_propagate1.f90 @@ -25,14 +25,15 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag) type (ele_struct), target :: ele1, ele2 type (twiss_struct) twiss_a type (lat_param_struct) param -type (coord_struct), pointer :: orb, orb_out -integer key2 +type (coord_struct), pointer :: orb_in, orb_out + +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) 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_p2 +real(rp) mat2(2,2), eta1_vec(6), eta_vec(6), vec(6), dpz2_dpz1, rel_p1, rel_p21, rel_p2 real(rp) det_factor, deriv_rel, gamma2_c, df logical error @@ -65,6 +66,9 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag) if (error) return endif +geometry = open$ +if (associated(ele1%branch)) geometry = ele1%branch%param%geometry + !--------------------------------------------------------------------- ! init @@ -98,20 +102,12 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag) ! -orb => ele2%map_ref_orb_in +orb_in => ele2%map_ref_orb_in orb_out => ele2%map_ref_orb_out -rel_p1 = 1 + orb%vec(6) ! reference energy +rel_p1 = 1 + orb_in%vec(6) ! Map reference momentum rel_p2 = 1 + orb_out%vec(6) mat6 = ele2%mat6 -if (ele2%key /= e_gun$) then ! Energy change normalization is not applied to an e-gun - if (bmad_private%normalize_twiss) then - mat6(:, 2:6:2) = mat6(:, 2:6:2) * rel_p1 - mat6(2:6:2, :) = mat6(2:6:2, :) / rel_p2 - endif - rel_p2 = rel_p2 / rel_p1 - rel_p1 = 1 -endif !--------------------------------------------------------------------- ! det_factor is a renormalization factor to handle non-symplectic errors. @@ -208,34 +204,27 @@ 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, ele1%y%eta, ele1%y%etap, 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. ! With the RF on, dz/dpz, along with the other dispersion components, is not well defined. -! This being the case, eta1_vec(5) is just treated as zero for an rfcavity. +! This being the case, eta1_vec(5) is just treated as zero for a closed branch -if (key2 == rfcavity$) eta1_vec(5) = 0 +eta1_vec = [ele1%x%eta, ele1%x%etap, ele1%y%eta, ele1%y%etap, ele1%z%eta, 1.0_rp] +if (geometry == closed$) eta1_vec(5) = 0 ! Must avoid 0/0 divide at zero reference momentum. -! If rel_p1 = 0 then total momentum is zero and orb%vec(2) and orb%vec(4) must be zero. +! If rel_p1 = 0 then total momentum is zero and orb_in%vec(2) and orb_in%vec(4) must be zero. ! Also for elements with a static electric field, pz should include the potential energy and so the ! mat6(6,:) terms should be all zero. However Bmad is not using proper canonical coords so the ! mat6(6,:) terms are forced to zero. -if (rel_p1 == 0) then - dpz2_dpz1 = dot_product(mat6(6,:), eta1_vec) +if (rel_p1 == 0 .or. key2 == rfcavity$ .or. key2 == lcavity$) then + dpz2_dpz1 = dot_product(mat6(6,:), eta1_vec) else - if (key2 == rfcavity$ .or. key2 == lcavity$) then - dpz2_dpz1 = dot_product(mat6(6,:), eta1_vec) - else - dpz2_dpz1 = 1 / rel_p2 - endif + dpz2_dpz1 = rel_p1 / rel_p2 endif -! - eta_vec(1:5) = matmul (mat6(1:5,:), eta1_vec) / dpz2_dpz1 ele2%x%eta = eta_vec(1) @@ -246,9 +235,14 @@ subroutine twiss_propagate1 (ele1, ele2, err_flag) 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 +if (geometry == closed$) then + ele2%z%eta = 0 +else + ele2%z%eta = eta_vec(5) +endif + +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)) diff --git a/bmad/doc/linear-optics.tex b/bmad/doc/linear-optics.tex index 5819418090..d11a90eb89 100644 --- a/bmad/doc/linear-optics.tex +++ b/bmad/doc/linear-optics.tex @@ -391,8 +391,8 @@ \section{Dispersion Calculation} The associated momentum dispersion is: \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 + \text{etap_x} = \eta_{px} \equiv \left. \frac{dp_x}{dp_z} \right|_s \comma \qquad + \text{etap_y} = \eta_{py} \equiv \left. \frac{dp_y}{dp_z} \right|_s \comma \qquad \end{equation} The momentum dispersion is useful when constructing particle bunch distributions and for various calculations like for calculating radiation integrals. @@ -408,20 +408,30 @@ \section{Dispersion Calculation} \text{deta_x_ds} &\equiv \frac{d\eta_x}{ds} = \frac{d}{dp_z} \left( \frac{dx}{ds} \right) = \frac{d}{dp_z} \left( \frac{p_x}{1 + p_z} \right) - = \frac{1}{1 + p_z} \, \frac{dp_x}{dp_z} - \frac{p_x}{(1 + p_z)^2} \CRNO + = \frac{1}{1 + p_z} \, \eta_{px} - \frac{p_x}{(1 + p_z)^2} \CRNO \text{deta_y_ds} &\equiv \frac{d\eta_y}{ds} = \frac{d}{dp_z} \left( \frac{dy}{ds} \right) = \frac{d}{dp_z} \left( \frac{p_y}{1 + p_z} \right) - = \frac{1}{1 + p_z} \, \frac{dp_y}{dp_z} - \frac{p_y}{(1 + p_z)^2} + = \frac{1}{1 + p_z} \, \eta_{py} - \frac{p_y}{(1 + p_z)^2} \label{dexds} \end{align} -For a non-circular machine, there are two ways one can imagine defining the dispersion: Either with -respect to changes in energy at the beginning of the machine or with respect to the local change in -energy at the point of measurement. The former definition will be called ``non-local dispersion'' -and the latter definition will be called ``local dispersion'' which what -\bmad calculates. The non-local dispersion $\wt\bfeta(s_1)$ at some point -$s_1$ is related to the local dispersion $\bfeta(s_1)$ via +For a lattice branch with an open (non-circular) geometry, the dispersion of the $z$ phase space +coordinate, $\eta_z$ can be defined similar to the dispersion of the other coordinates. In this +case, the dispersion vector $\bfeta$ is defined by +\begin{equation} + \bfeta = (\eta_x, \eta_{px}, \eta_{y}, \eta_{py}, \eta_z, 1) +\end{equation} +and this vector is propagated via +\begin{equation} +\end{equation + +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 +with respect to the local change in energy at the point of measurement. The former definition will +be called ``non-local dispersion'' and the latter definition will be called ``local dispersion'' +which what \bmad calculates. The non-local dispersion $\wt\bfeta(s_1)$ at some point $s_1$ is +related to the local dispersion $\bfeta(s_1)$ via \begin{equation} \wt\bfeta(s_1) = \frac{dp_{z1}}{dp_{z0}} \, \bfeta(s_1) \end{equation} @@ -430,3 +440,18 @@ \section{Dispersion Calculation} dispersion, on the other hand, reflects the correlations between the particle energy and particle position within a beam. +For a closed geometry lattice branch, defining the dependence of $z$ on $p_z$ is problematical. +With the RF off, $z$ is not periodic so a closed orbit $z$ cannot be defined. With the RF on, the +dispersion of any of the phase space components is not well defined. This being the case, $\eta_z$ +is just treated as zero for a closed branch. + +Note: For a closed geometry branch with RF on, it is possible to define dispersions. If $\bfv$ is +the eigenvector of the eigenmode associated with longitudinal oscillations, the dispersion $\eta_x$ +can be defined by $\bfv(1) / \bfv(6)$ with similar definitions for the other dispersion components. +With this definition, the dispersion become complex. In the low RF limit, the dispersions $\eta_x$, +$\eta_{px}$, $\eta_y$, $\eta_{py}$ converge to the standard (real) values and $\eta_z$ diverges +to infinity.\footnote + { +This is assuming a linear system. In practice, the motion will become unstable due to +the finite size of the RF bucket. + } diff --git a/bmad/low_level/transfer_twiss.f90 b/bmad/low_level/transfer_twiss.f90 index e9989cc858..d3201f18d9 100644 --- a/bmad/low_level/transfer_twiss.f90 +++ b/bmad/low_level/transfer_twiss.f90 @@ -3,6 +3,7 @@ ! ! Routine to transfer the Twiss, coupling, and dispersion parameters from one element to another. ! +! Note: %map_ref_orb_out is transferred. ! ! Input: ! ele_in -- ele_struct: Element with existing Twiss parameters. @@ -32,6 +33,7 @@ subroutine transfer_twiss (ele_in, ele_out, reverse) ele_out%c_mat = ele_in%c_mat ele_out%gamma_c = ele_in%gamma_c ele_out%mode_flip = ele_in%mode_flip +ele_out%map_ref_orb_out = ele_in%map_ref_orb_out if (logic_option(.false., reverse)) then ele_out%x%etap = -ele_in%x%etap diff --git a/bmad/modules/bmad_routine_interface.f90 b/bmad/modules/bmad_routine_interface.f90 index f55370857c..8ef8e27325 100644 --- a/bmad/modules/bmad_routine_interface.f90 +++ b/bmad/modules/bmad_routine_interface.f90 @@ -2247,13 +2247,13 @@ subroutine set_status_flags (bookkeeping_state, stat) integer stat end subroutine -subroutine set_twiss(branch, twiss_ele, ix_ele, err_flag, print_err) +subroutine set_twiss(branch, twiss_ele, ix_ele, match_deta_ds, err_flag, print_err) import implicit none type (branch_struct), target :: branch type (ele_struct) twiss_ele integer ix_ele - logical err_flag + logical match_deta_ds, err_flag logical, optional :: print_err end subroutine diff --git a/bmad/modules/bmad_struct.f90 b/bmad/modules/bmad_struct.f90 index 6d23fd95e9..b3591b26aa 100644 --- a/bmad/modules/bmad_struct.f90 +++ b/bmad/modules/bmad_struct.f90 @@ -2241,7 +2241,6 @@ module bmad_struct type bmad_private_struct real(rp) :: rf_clock_period = 0 ! The RF clock is used by the long_term_tracking program to avoid time round-off errors. - logical :: normalize_twiss = .false. ! Experimental: Normalize Twiss calc by the ref pz about which the 1-turn matrix is calculated? end type type (bmad_private_struct), save, target :: bmad_private diff --git a/bmad/parsing/set_ele_defaults.f90 b/bmad/parsing/set_ele_defaults.f90 index c050fd9542..9585d06214 100644 --- a/bmad/parsing/set_ele_defaults.f90 +++ b/bmad/parsing/set_ele_defaults.f90 @@ -60,6 +60,8 @@ subroutine set_ele_defaults (ele, do_allocate) ele%value(p0c$) = -1 ele%value(inherit_from_fork$) = real_garbage$ ele%value(deta_ds_master$) = false$ + ele%z%etap = 1 + ele%z%deta_ds = 1 call mat_make_unit (ele%mat6) case (capillary$) diff --git a/regression_tests/cesr_test/output.correct b/regression_tests/cesr_test/output.correct index 70c6439697..baccea97e3 100644 --- a/regression_tests/cesr_test/output.correct +++ b/regression_tests/cesr_test/output.correct @@ -8,13 +8,13 @@ "Closed Orb 4 Del" ABS 1e-12 -6.9784E-15 -4.6018E-15 4.1413E-16 1.0235E-15 1.2878E-12 1.6941E-21 "First Descrip" STR "First" "Second Descrip" STR "Second" -"Cache Diff: I1-no_wig" ABS 1e-8 496 6.179560494E-04 +"Cache Diff: I1-no_wig" ABS 1e-8 496 6.179560500E-04 "Cache Diff: I2-no_wig" ABS 1e-8 10 5.551230207E-06 "Cache Diff: I3-no_wig" ABS 1e-8 10 1.523799473E-06 -"Cache Diff: I4a-no_wig" ABS 1e-8 156 1.405752967E-04 -"Cache Diff: I4b-no_wig" ABS 1e-8 428 2.912442458E-04 -"Cache Diff: I5a-no_wig" ABS 1e-8 824 8.203538367E-03 -"Cache Diff: I5b-no_wig" ABS 1e-8 402 9.804946516E-04 +"Cache Diff: I4a-no_wig" ABS 1e-8 156 1.405752968E-04 +"Cache Diff: I4b-no_wig" ABS 1e-8 428 2.912442463E-04 +"Cache Diff: I5a-no_wig" ABS 1e-8 824 8.203538350E-03 +"Cache Diff: I5b-no_wig" ABS 1e-8 402 9.804946504E-04 "Dif:Orb(1)" ABS 1.0E-07 -2.33685141549922E-14 "Dif:Orb(2)" ABS 1.0E-07 4.77123549003089E-13 "Dif:Orb(3)" ABS 1.0E-07 1.45020108558759E-15 @@ -25,29 +25,29 @@ "Dif:Beta_Y" ABS 1.0E-07 3.55004957722227E-09 "Dif:Alpha_X" ABS 1.0E-07 1.09135069879121E-05 "Dif:Alpha_Y" ABS 1.0E-07 1.48104998647092E-07 -"Dif:Eta_X " ABS 1.0E-07 4.28190161340315E-06 -"Dif:Eta_Y " ABS 1.0E-07 3.40306164975032E-09 -"Dif:Etap_X" ABS 1.0E-06 -3.04067295037319E-05 -"Dif:Etap_Y" ABS 1.0E-06 2.26378816960073E-06 -"Dif:Etap_X" ABS 1.0E-06 -3.03098643232719E-05 -"Dif:Etap_Y" ABS 1.0E-06 -9.73096473493506E-07 +"Dif:Eta_X " ABS 1.0E-07 4.28192300504812E-06 +"Dif:Eta_Y " ABS 1.0E-07 3.40292745522851E-09 +"Dif:Etap_X" ABS 1.0E-06 -3.04071629210270E-05 +"Dif:Etap_Y" ABS 1.0E-06 2.26380020518007E-06 +"Dif:Etap_X" ABS 1.0E-06 -3.04073952120190E-05 +"Dif:Etap_Y" ABS 1.0E-06 2.16063616251174E-06 "Lat1:Beta_a" ABS 1.0E-06 5.74755728175682E+00 "Lat1:Alpha_a" ABS 1.0E-06 6.17183084954029E-01 -"Lat1:Eta_a" ABS 1.0E-06 3.81079580272918E-01 -"Lat1:Etap_a" ABS 1.0E-06 -1.50985407173454E-01 -"Lat1:Deta_a_Ds" ABS 1.0E-06 -1.49500566944722E-01 -"Lat1:Eta_x" ABS 1.0E-06 3.81072825907900E-01 -"Lat1:Etap_x" ABS 1.0E-06 -1.50984512302200E-01 -"Lat1:Deta_x_Ds" ABS 1.0E-06 -1.49499675027797E-01 +"Lat1:Eta_a" ABS 1.0E-06 3.81079580219623E-01 +"Lat1:Etap_a" ABS 1.0E-06 -1.50985407074940E-01 +"Lat1:Deta_a_Ds" ABS 1.0E-06 -1.49491430805345E-01 +"Lat1:Eta_x" ABS 1.0E-06 3.81072825854607E-01 +"Lat1:Etap_x" ABS 1.0E-06 -1.50984512203686E-01 +"Lat1:Deta_x_Ds" ABS 1.0E-06 -1.49490538943290E-01 "Lat1:Phi_a" ABS 1.0E-06 6.67148580318659E+00 "Lat1:Beta_b" ABS 1.0E-06 3.00630062638158E+01 "Lat1:Alpha_y" ABS 1.0E-06 -1.70228379556182E+00 -"Lat1:Eta_b" ABS 1.0E-06 -2.69791050930001E-02 -"Lat1:Etap_b" ABS 1.0E-06 -1.86491344564010E-03 -"Lat1:Deta_b_Ds" ABS 1.0E-06 -1.85833101000457E-03 -"Lat1:Eta_y" ABS 1.0E-06 -2.26187288392922E-02 -"Lat1:Etap_y" ABS 1.0E-06 -1.53834524911109E-03 -"Lat1:Deta_y_Ds" ABS 1.0E-06 -1.53513786617504E-03 +"Lat1:Eta_b" ABS 1.0E-06 -2.69791050794704E-02 +"Lat1:Etap_b" ABS 1.0E-06 -1.86491344469377E-03 +"Lat1:Deta_b_Ds" ABS 1.0E-06 -1.85821671279631E-03 +"Lat1:Eta_y" ABS 1.0E-06 -2.26187288294449E-02 +"Lat1:Etap_y" ABS 1.0E-06 -1.53834524844234E-03 +"Lat1:Deta_y_Ds" ABS 1.0E-06 -1.53504330951996E-03 "Lat1:Phi_y" ABS 1.0E-06 5.93989686093169E+00 "Lat1:Orb X" ABS 1.0E-10 3.25325012796147E-03 "Lat1:Orb P_X" ABS 1.0E-10 -1.48509969948590E-03 @@ -57,66 +57,66 @@ "Lat1:Orb P_Z" ABS 1.0E-10 5.99542531502537E-05 "Lat1:Chrom_x" ABS 1.0E-05 -1.17355450638179E+00 "Lat1:Chrom_y" ABS 1.0E-05 -9.41786942512124E-01 -"Lat1:Synch_int(1)" ABS 1.0E-06 8.84194890491963E+00 +"Lat1:Synch_int(1)" ABS 1.0E-06 8.84194890492006E+00 "Lat1:Synch_int(2)" ABS 1.0E-06 1.04424863370109E-01 "Lat1:Synch_int(3)" ABS 1.0E-06 2.30420629537313E-03 -"Lat1:Sige_e" ABS 1.0E-10 6.75242117495030E-04 -"Lat1:Sig_z" ABS 1.0E-08 1.85933501254008E-02 +"Lat1:Sige_e" ABS 1.0E-10 6.75242117421361E-04 +"Lat1:Sig_z" ABS 1.0E-08 1.85933501233728E-02 "Lat1:E_loss" ABS 1.0E-01 1.15047899411770E+06 -"Lat1:A%Emittance" ABS 1.0E-12 2.02129364313202E-07 -"Lat1:B%Emittance" ABS 1.0E-14 1.79689109249269E-11 -"Lat1:Z%Emittance" ABS 1.0E-11 1.25550131100021E-05 -"Lat1:A%Synch_int(4)" ABS 1.0E-07 -1.35181830461404E-03 -"Lat1:A%Synch_int(5)" ABS 1.0E-07 5.20827232373484E-04 -"Lat1:B%Synch_int(4)" ABS 1.0E-07 -4.09638242913126E-05 -"Lat1:B%Synch_int(5)" ABS 1.0E-11 4.56385534427079E-08 -"Lat1:Z%Synch_int(4)" ABS 1.0E-08 -1.39278212890535E-03 +"Lat1:A%Emittance" ABS 1.0E-12 2.02129364314750E-07 +"Lat1:B%Emittance" ABS 1.0E-14 1.79689109250504E-11 +"Lat1:Z%Emittance" ABS 1.0E-11 1.25550131072629E-05 +"Lat1:A%Synch_int(4)" ABS 1.0E-07 -1.35181825937264E-03 +"Lat1:A%Synch_int(5)" ABS 1.0E-07 5.20827232154712E-04 +"Lat1:B%Synch_int(4)" ABS 1.0E-07 -4.09638242656544E-05 +"Lat1:B%Synch_int(5)" ABS 1.0E-11 4.56385534430111E-08 +"Lat1:Z%Synch_int(4)" ABS 1.0E-08 -1.39278208363830E-03 "Cache Diff: I1-wig" ABS 1e-8 398 6.427705296E-04 "Cache Diff: I2-wig" ABS 1e-8 917 1.796173843E-06 -"Cache Diff: I3-wig" ABS 1e-8 917 2.004272990E-07 -"Cache Diff: I4a-wig" ABS 1e-8 56 5.063753845E-05 +"Cache Diff: I3-wig" ABS 1e-8 917 2.004272991E-07 +"Cache Diff: I4a-wig" ABS 1e-8 56 5.063753843E-05 "Cache Diff: I4b-wig" ABS 1e-8 917 1.234030878E-04 -"Cache Diff: I5a-wig" ABS 1e-8 50 2.867657460E-05 -"Cache Diff: I5b-wig" ABS 1e-8 862 1.389303588E-05 -"Lat2:Beta_a" ABS 1.0E-05 1.63241327696990E+01 -"Lat2:Alpha_a" ABS 1.0E-05 -1.51916901683308E+00 -"Lat2:Eta_a" ABS 1.0E-05 2.25296842978666E+00 -"Lat2:Etap_a" ABS 1.0E-05 1.86023795976717E-01 -"Lat2:Deta_a_Ds" ABS 1.0E-05 1.85471529599840E-01 -"Lat2:Eta_x" ABS 1.0E-05 2.25293014970256E+00 -"Lat2:Etap_x" ABS 1.0E-05 1.86006226079607E-01 -"Lat2:Deta_x_Ds" ABS 1.0E-05 1.85453998341411E-01 -"Lat2:Phi_a" ABS 1.0E-05 5.68599666766945E+00 -"Lat2:Beta_b" ABS 1.0E-05 2.65318213916511E+01 -"Lat2:Alpha_b" ABS 1.0E-05 1.51859630159153E+00 -"Lat2:Eta_b" ABS 1.0E-05 -8.81652263108863E-03 -"Lat2:Etap_b" ABS 1.0E-05 -7.50700635275607E-04 -"Lat2:Deta_b_Ds" ABS 1.0E-05 -7.46691029978985E-04 -"Lat2:Eta_y" ABS 1.0E-05 -2.03201128939941E-02 -"Lat2:Etap_y" ABS 1.0E-05 1.95191415148415E-03 -"Lat2:Deta_y_Ds" ABS 1.0E-05 1.94624367073931E-03 -"Lat2:Phi_b" ABS 1.0E-05 4.76647532786524E+00 +"Cache Diff: I5a-wig" ABS 1e-8 50 2.867657461E-05 +"Cache Diff: I5b-wig" ABS 1e-8 862 1.389303589E-05 +"Lat2:Beta_a" ABS 1.0E-05 1.63241327697001E+01 +"Lat2:Alpha_a" ABS 1.0E-05 -1.51916901683317E+00 +"Lat2:Eta_a" ABS 1.0E-05 2.25296842978664E+00 +"Lat2:Etap_a" ABS 1.0E-05 1.86023795976713E-01 +"Lat2:Deta_a_Ds" ABS 1.0E-05 1.85471761338328E-01 +"Lat2:Eta_x" ABS 1.0E-05 2.25293014970255E+00 +"Lat2:Etap_x" ABS 1.0E-05 1.86006226079604E-01 +"Lat2:Deta_x_Ds" ABS 1.0E-05 1.85454230057977E-01 +"Lat2:Phi_a" ABS 1.0E-05 5.68599666766948E+00 +"Lat2:Beta_b" ABS 1.0E-05 2.65318213916517E+01 +"Lat2:Alpha_b" ABS 1.0E-05 1.51859630159158E+00 +"Lat2:Eta_b" ABS 1.0E-05 -8.81652263112590E-03 +"Lat2:Etap_b" ABS 1.0E-05 -7.50700635274169E-04 +"Lat2:Deta_b_Ds" ABS 1.0E-05 -7.46691960696562E-04 +"Lat2:Eta_y" ABS 1.0E-05 -2.03201128939990E-02 +"Lat2:Etap_y" ABS 1.0E-05 1.95191415148456E-03 +"Lat2:Deta_y_Ds" ABS 1.0E-05 1.94624610264159E-03 +"Lat2:Phi_b" ABS 1.0E-05 4.76647532786521E+00 "Lat2:Orb X" ABS 1.0E-07 1.18493311706587E-02 -"Lat2:Orb P_X" ABS 1.0E-07 5.52227738196111E-04 -"Lat2:Orb Y" ABS 1.0E-07 -2.65182968854781E-05 -"Lat2:Orb P_Y" ABS 1.0E-07 5.67048074484464E-06 +"Lat2:Orb P_X" ABS 1.0E-07 5.52227738196148E-04 +"Lat2:Orb Y" ABS 1.0E-07 -2.65182968855430E-05 +"Lat2:Orb P_Y" ABS 1.0E-07 5.67048074485485E-06 "Lat2:Orb Z" ABS 1.0E-07 -1.49778902918623E-03 "Lat2:Orb P_Z" ABS 1.0E-07 -1.25318583885772E-06 -"Lat2:Chrom_x" ABS 1.0E-04 6.41525629895057E+00 -"Lat2:Chrom_y" ABS 1.0E-04 5.70528459498121E+00 -"Lat2:Synch_int(1)" ABS 1.0E-06 8.66018423912541E+00 +"Lat2:Chrom_x" ABS 1.0E-04 6.41525629885287E+00 +"Lat2:Chrom_y" ABS 1.0E-04 5.70528459500785E+00 +"Lat2:Synch_int(1)" ABS 1.0E-06 8.66018423912542E+00 "Lat2:Synch_int(2)" ABS 1.0E-06 9.97053316812563E-01 "Lat2:Synch_int(3)" ABS 1.0E-06 2.62658509891663E-01 "Lat2:Sige_e" ABS 1.0E-10 8.41149646134634E-04 "Lat2:Sig_z" ABS 1.0E-08 2.27386275458037E-02 "Lat2:E_loss" ABS 1.0E-01 1.77232886436639E+05 -"Lat2:A%Emittance" ABS 1.0E-12 1.23659335790085E-07 -"Lat2:B%Emittance" ABS 1.0E-14 9.54067929108944E-11 +"Lat2:A%Emittance" ABS 1.0E-12 1.23659335790078E-07 +"Lat2:B%Emittance" ABS 1.0E-14 9.54067929109133E-11 "Lat2:Z%Emittance" ABS 1.0E-11 1.91265885137401E-05 -"Lat2:A%Synch_int(4)" ABS 1.0E-07 -5.90562374954835E-02 -"Lat2:A%Synch_int(5)" ABS 1.0E-07 2.50457845985694E-02 -"Lat2:B%Synch_int(4)" ABS 1.0E-07 6.85852825837135E-04 -"Lat2:B%Synch_int(5)" ABS 1.0E-11 1.81627016056473E-05 -"Lat2:Z%Synch_int(4)" ABS 1.0E-08 -5.83703846696463E-02 +"Lat2:A%Synch_int(4)" ABS 1.0E-07 -5.90562374954820E-02 +"Lat2:A%Synch_int(5)" ABS 1.0E-07 2.50457845985681E-02 +"Lat2:B%Synch_int(4)" ABS 1.0E-07 6.85852825836967E-04 +"Lat2:B%Synch_int(5)" ABS 1.0E-11 1.81627016056510E-05 +"Lat2:Z%Synch_int(4)" ABS 1.0E-08 -5.83703846696451E-02 "Lat2:Lat" STR "T" diff --git a/regression_tests/twiss_test/twiss_test.bmad b/regression_tests/twiss_test/twiss_test.bmad index a64e8773cb..780ad5ec8f 100644 --- a/regression_tests/twiss_test/twiss_test.bmad +++ b/regression_tests/twiss_test/twiss_test.bmad @@ -10,48 +10,100 @@ beginning[beta_b] = 10 beginning[alpha_a] = +0.1 beginning[alpha_b] = -0.2 + +!------------------------- + +beginning[eta_x] = 0. +beginning[etap_x] = 0. +beginning[eta_y] = 0. +beginning[etap_y] = -0. +beginning[eta_z] = 0. + +beginning[cmat_11] = 0. +beginning[cmat_12] = 0. +beginning[cmat_21] = 0. +beginning[cmat_22] = 0. + +particle_start[x] = 0.00 +particle_start[px] = 0.00 +particle_start[y] = 0.00 +particle_start[py] = 0.00 +particle_start[z] = 0. +particle_start[pz] = 0. + + +!------------------------- + beginning[eta_x] = 0.1 beginning[etap_x] = 0.2 -beginning[eta_y] = 0.1 -beginning[etap_y] = 0.2 +beginning[eta_y] = 0.3 +beginning[etap_y] = -0.2 +beginning[eta_z] = 0.2 -!beginning[cmat_11] = 0.1 -!beginning[cmat_12] = 0.2 -!beginning[cmat_21] = 0.3 -!beginning[cmat_22] = 0.4 +beginning[cmat_11] = 0.1 +beginning[cmat_12] = 0.2 +beginning[cmat_21] = 0.3 +beginning[cmat_22] = 0.4 -!particle_start[x] = 0.001 +particle_start[x] = 0.001 particle_start[px] = 0.002 -!particle_start[y] = 0.003 -!particle_start[py] = 0.004 -!particle_start[z] = 0.1 -!particle_start[pz] = 0.5 +particle_start[y] = 0.003 +particle_start[py] = 0.004 +particle_start[z] = 0.1 +particle_start[pz] = 0.5 + +!------------------------- + +beginning[eta_x] = 0. +beginning[etap_x] = 0. +beginning[eta_y] = 0. +beginning[etap_y] = -0. +beginning[eta_z] = 0. + +beginning[cmat_11] = 0. +beginning[cmat_12] = 0. +beginning[cmat_21] = 0. +beginning[cmat_22] = 0. + +particle_start[x] = 0.001 +particle_start[px] = 0.00 +particle_start[y] = 0.00 +particle_start[py] = 0.00 +particle_start[z] = 0. +particle_start[pz] = 0 + + +!------------------------- a: quadrupole, l = 1, k1 = 0.3, tilt = 1 -z: quadrupole, l = 1, k1 = 0.1, tilt = 2 +z: quadrupole, l = 1, k1 = -0.1, tilt = 2 -ln: line = (a, z) +ln0: line = (a, z) !------------------------- -rf1: rfcavity, l = 1, rf_frequency = 1e9, tracking_method = bmad_standard, voltage = 1e10, phi0 = pi/2 +rf1: rfcavity, l = 1, rf_frequency = 1e9, tracking_method = runge_kutta, voltage = 4e10, phi0 = 0.3 eg: e_gun, l = 1, tracking_method = runge_kutta, voltage = 0e5 -z2: quadrupole, l = 0, k1 = 0.1, tilt = 2 +z2: quadrupole, l = 1, k1 = 0.1, tilt = 2 -lneg: line = (rf1, z2) -lneg[geometry] = open -lneg[particle] = electron -lneg[p0c] = 1e11 -lneg[beta_a] = 10 -lneg[beta_b] = 10 -lneg[alpha_a] = +0.1 -lneg[alpha_b] = -0.2 +lnrf: line = (rf1, z2) +lnrf[geometry] = open +lnrf[particle] = electron +lnrf[p0c] = 1e11 +lnrf[beta_a] = 10 +lnrf[beta_b] = 10 +lnrf[alpha_a] = +0.1 +lnrf[alpha_b] = -0.2 !------------------------- -use, lneg +use, lnrf no_digested +! To Do: +! e_gun +! electric multipole +! lcavity diff --git a/regression_tests/twiss_test/twiss_test.f90 b/regression_tests/twiss_test/twiss_test.f90 index 65c0afdcf7..26e365d0f8 100644 --- a/regression_tests/twiss_test/twiss_test.f90 +++ b/regression_tests/twiss_test/twiss_test.f90 @@ -6,75 +6,103 @@ program twiss_test type (lat_struct), target :: lat type (branch_struct), pointer :: branch -type (ele_struct) ele0 -type (ele_struct), pointer :: elen -type (coord_struct), allocatable :: orbit(:) +type (ele_struct) t_ele, ds_ele +type (ele_struct), pointer :: ele0, elen +type (coord_struct), allocatable :: orb0(:), orb1(:) type (twiss_struct) dt, tt type (xy_disp_struct) dxy, xy -real(rp) dc(2,2), max_diff +real(rp) dc(2,2), max_diff, dpz, eta_vec(6), dorb(6) integer n, ib, status logical err character(200) lat_file +character(3) b_str -! Idea is to check set_twiss +! Check twiss_propagate vs true dispersion calculated via tracking. open (1, file = 'output.now') +dpz = 1d-8 lat_file = 'twiss_test.bmad' call bmad_parser(lat_file, lat) do ib = 0, ubound(lat%branch, 1) branch => lat%branch(ib) - call reallocate_coord(orbit, branch%n_ele_max) - call twiss_and_track(lat, orbit, status, ib, .true.) + b_str = 'br' // int_str(ib) + ele0 => branch%ele(0) + call reallocate_coord(orb0, branch%n_ele_max) + call twiss_and_track(lat, orb0, status, ib, .true.) n = branch%n_ele_track - elen => branch%ele(n-1) - call transfer_twiss(branch%ele(n), ele0) - call set_twiss(branch, ele0, n-1, err) + eta_vec = [ele0%x%eta, ele0%x%etap, ele0%y%eta, ele0%y%etap, ele0%z%eta, 1.0_rp] + + call reallocate_coord(orb1, branch%n_ele_max) + call init_coord(orb1(0), orb0(0)%vec + dpz * eta_vec, ele0, downstream_end$) + call track_all(lat, orb1, ib) + + dorb = (orb1(n)%vec - orb0(n)%vec) / (orb1(n)%vec(6) - orb0(n)%vec(6)) + elen => branch%ele(n) + eta_vec = [elen%x%eta, elen%x%etap, elen%y%eta, elen%y%etap, elen%z%eta, 1.0_rp] - max_diff = 0 if (ib /= 0) write (1, '(a)') + write (1, '(2a, 6es14.6)') quote(b_str // '-Fwd-Orb'), ' ABS 1e-7', dorb + write (1, '(2a, 6es14.6)') quote(b_str // '-Fwd-dOrb'), ' ABS 1e-14', dorb - eta_vec + max_diff = maxval(abs(dorb - eta_vec)) + write (1, '(2a, es14.6)') quote(b_str // '-dMax-Fwd'), ' ABS 1e-7', max_diff - tt = ele0%a - dt = twiss_diff(ele0%a, elen%a, max_diff) - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-a-Twiss'), ' REL 1e-7', tt%beta, tt%alpha, tt%gamma - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-a-dTwiss'), ' ABS 1E-14', dt%beta, dt%alpha, dt%gamma - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-a-Disp'), ' REL 1e-7', tt%eta, tt%etap, tt%deta_ds - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-a-dDisp'), ' ABS 1E-14', dt%eta, dt%etap, dt%deta_ds - - tt = ele0%b - dt = twiss_diff(ele0%b, elen%b, max_diff) - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-b-Twiss'), ' REL 1e-7', tt%beta, tt%alpha, tt%gamma - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-b-dTwiss'), ' ABS 1E-14', dt%beta, dt%alpha, dt%gamma - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-b-Disp'), ' REL 1e-7', tt%eta, tt%etap, tt%deta_ds - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-b-dDisp'), ' ABS 1E-14', dt%eta, dt%etap, dt%deta_ds - - tt = ele0%z - dt = twiss_diff(ele0%z, elen%z, max_diff) - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-z-Twiss'), ' REL 1e-7', tt%beta, tt%alpha, tt%gamma - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-z-dTwiss'), ' ABS 1E-14', dt%beta, dt%alpha, dt%gamma - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-z-Disp'), ' REL 1e-7', tt%eta, tt%etap, tt%deta_ds - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-z-dDisp'), ' ABS 1E-14', dt%eta, dt%etap, dt%deta_ds - - dc = ele0%c_mat - elen%c_mat - max_diff = max(max_diff, maxval(abs(dc))) - write (1, '(2a, 4es14.6)') quote('br' // int_str(ib) // '-C_mat'), ' ABS 1e-7', ele0%c_mat - write (1, '(2a, 4es14.6)') quote('br' // int_str(ib) // '-dC_mat'), ' ABS 1e-14', dc + ! Check of set_twiss by setting Twiss at ele n-2 with the Twiss at ele n-1. + ! Note: (a,b) %deta_ds are not matched to. - xy = ele0%x - dxy = xy_disp_diff(ele0%x, elen%x, max_diff) - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-x-Disp'), ' REL 1e-7', xy%eta, xy%etap, xy%deta_ds - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-x-dDisp'), ' ABS 1E-14', dxy%eta, dxy%etap, dxy%deta_ds + call transfer_twiss(branch%ele(n-1), t_ele) + call set_twiss(branch, t_ele, n-2, .true., err) + elen => branch%ele(n-2) + call transfer_twiss(elen, ds_ele) + call set_twiss(branch, t_ele, n-2, .false., err) - xy = ele0%y - dxy = xy_disp_diff(ele0%y, elen%y, max_diff) - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-y-Disp'), ' REL 1e-7', xy%eta, xy%etap, xy%deta_ds - write (1, '(2a, 3es14.6)') quote('br' // int_str(ib) // '-y-dDisp'), ' ABS 1E-14', dxy%eta, dxy%etap, dxy%deta_ds + max_diff = 0 + tt = t_ele%a + dt = twiss_diff(t_ele%a, ds_ele%a, elen%a, max_diff) + write (1, '(a)') + write (1, '(2a, 3es14.6)') quote(b_str // '-a-Twiss'), ' REL 1e-7', tt%beta, tt%alpha, tt%gamma + write (1, '(2a, 3es14.6)') quote(b_str // '-a-dTwiss'), ' ABS 1E-14', dt%beta, dt%alpha, dt%gamma + write (1, '(2a, 3es14.6)') quote(b_str // '-a-Disp'), ' REL 1e-7', tt%eta, tt%etap, tt%deta_ds + write (1, '(2a, 3es14.6)') quote(b_str // '-a-dDisp'), ' ABS 1E-14', dt%eta, dt%etap ! , dt%deta_ds + + tt = t_ele%b + dt = twiss_diff(t_ele%b, ds_ele%b, elen%b, max_diff) + write (1, '(a)') + write (1, '(2a, 3es14.6)') quote(b_str // '-b-Twiss'), ' REL 1e-7', tt%beta, tt%alpha, tt%gamma + write (1, '(2a, 3es14.6)') quote(b_str // '-b-dTwiss'), ' ABS 1E-14', dt%beta, dt%alpha, dt%gamma + write (1, '(2a, 3es14.6)') quote(b_str // '-b-Disp'), ' REL 1e-7', tt%eta, tt%etap, tt%deta_ds + write (1, '(2a, 3es14.6)') quote(b_str // '-b-dDisp'), ' ABS 1E-14', dt%eta, dt%etap ! , dt%deta_ds + + tt = t_ele%z + dt = twiss_diff(t_ele%z, ds_ele%z, elen%z, max_diff) + write (1, '(a)') + write (1, '(2a, 3es14.6)') quote(b_str // '-z-Twiss'), ' REL 1e-7', tt%beta, tt%alpha, tt%gamma + write (1, '(2a, 3es14.6)') quote(b_str // '-z-dTwiss'), ' ABS 1E-14', dt%beta, dt%alpha, dt%gamma + write (1, '(2a, 3es14.6)') quote(b_str // '-z-Disp'), ' REL 1e-7', tt%eta, tt%etap, tt%deta_ds + write (1, '(2a, 3es14.6)') quote(b_str // '-z-dDisp'), ' ABS 1E-14', dt%eta, dt%etap, dt%deta_ds + + xy = t_ele%x + dxy = xy_disp_diff(t_ele%x, ds_ele%x, elen%x, max_diff) + write (1, '(a)') + write (1, '(2a, 3es14.6)') quote(b_str // '-x-Disp'), ' REL 1e-7', xy%eta, xy%etap, xy%deta_ds + write (1, '(2a, 3es14.6)') quote(b_str // '-x-dDisp'), ' ABS 1E-14', dxy%eta, dxy%etap, dxy%deta_ds + + xy = t_ele%y + dxy = xy_disp_diff(t_ele%y, ds_ele%y, elen%y, max_diff) + write (1, '(2a, 3es14.6)') quote(b_str // '-y-Disp'), ' REL 1e-7', xy%eta, xy%etap, xy%deta_ds + write (1, '(2a, 3es14.6)') quote(b_str // '-y-dDisp'), ' ABS 1E-14', dxy%eta, dxy%etap, dxy%deta_ds + + dc = t_ele%c_mat - elen%c_mat + max_diff = max(max_diff, maxval(abs(dc))) + write (1, '(a)') + write (1, '(2a, 4es14.6)') quote(b_str // '-C_mat'), ' ABS 1e-7', t_ele%c_mat + write (1, '(2a, 4es14.6)') quote(b_str // '-dC_mat'), ' ABS 1e-14', dc - write (1, '(2a, 2l3, a)') quote('br' // int_str(ib) // '-Flib'), ' STR "', ele0%mode_flip, ele0%mode_flip .eqv. elen%mode_flip, '"' - write (1, '(2a, es14.6)') quote('br' // int_str(ib) // '-dMax'), ' ABS 1e-7', max_diff + write (1, '(2a, 2l1, a)') quote(b_str // '-flip'), ' STR "', t_ele%mode_flip, t_ele%mode_flip .eqv. elen%mode_flip, '"' + write (1, '(2a, es14.6)') quote(b_str // '-dMax-Rev'), ' ABS 1e-7', max_diff enddo @@ -84,9 +112,9 @@ program twiss_test !-------------------------------------------------------------------------------------------- contains -function twiss_diff(twiss1, twiss2, max_diff) result (twiss_d) +function twiss_diff(twiss1, ds_twiss, twiss2, max_diff) result (twiss_d) -type (twiss_struct) twiss1, twiss2, twiss_d +type (twiss_struct) twiss1, ds_twiss, twiss2, twiss_d real(rp) max_diff ! @@ -97,28 +125,28 @@ function twiss_diff(twiss1, twiss2, max_diff) result (twiss_d) twiss_d%phi = twiss1%phi - twiss2%phi twiss_d%eta = twiss1%eta - twiss2%eta twiss_d%etap = twiss1%etap - twiss2%etap -twiss_d%deta_ds = twiss1%deta_ds - twiss2%deta_ds +twiss_d%deta_ds = twiss1%deta_ds - ds_twiss%deta_ds twiss_d%sigma = twiss1%sigma - twiss2%sigma twiss_d%sigma_p = twiss1%sigma_p - twiss2%sigma_p twiss_d%emit = twiss1%emit - twiss2%emit twiss_d%norm_emit = twiss1%norm_emit - twiss2%norm_emit max_diff = max(max_diff, abs(twiss_d%beta), abs(twiss_d%alpha), abs(twiss_d%gamma), & - abs(twiss_d%eta), abs(twiss_d%etap), abs(twiss_d%deta_ds)) + abs(twiss_d%eta), abs(twiss_d%etap)) ! , abs(twiss_d%deta_ds)) end function twiss_diff !-------------------------------------------------------------------------------------------- ! contains -function xy_disp_diff(xy_disp1, xy_disp2, max_diff) result (xy_disp_d) +function xy_disp_diff(xy_disp1, ds_xy, xy_disp2, max_diff) result (xy_disp_d) -type (xy_disp_struct) xy_disp1, xy_disp2, xy_disp_d +type (xy_disp_struct) xy_disp1, ds_xy, xy_disp2, xy_disp_d real(rp) max_diff xy_disp_d%eta = xy_disp1%eta - xy_disp2%eta xy_disp_d%etap = xy_disp1%etap - xy_disp2%etap -xy_disp_d%deta_ds = xy_disp1%deta_ds - xy_disp2%deta_ds +xy_disp_d%deta_ds = xy_disp1%deta_ds - ds_xy%deta_ds xy_disp_d%sigma = xy_disp1%sigma - xy_disp2%sigma max_diff = max(max_diff, abs(xy_disp_d%eta), abs(xy_disp_d%etap), abs(xy_disp_d%deta_ds)) diff --git a/tao/version/tao_version_mod.f90 b/tao/version/tao_version_mod.f90 index ff4a5fe005..1a28f65d14 100644 --- a/tao/version/tao_version_mod.f90 +++ b/tao/version/tao_version_mod.f90 @@ -6,5 +6,5 @@ !- module tao_version_mod -character(*), parameter :: tao_version_date = "2024/06/02 11:46:56" +character(*), parameter :: tao_version_date = "2024/06/03 14:07:52" end module