diff --git a/bmad/doc/bend-fiducial1.pdf b/bmad/doc/bend-fiducial1.pdf new file mode 100644 index 0000000000..686f12aea9 Binary files /dev/null and b/bmad/doc/bend-fiducial1.pdf differ diff --git a/bmad/doc/bend-fiducial2.pdf b/bmad/doc/bend-fiducial2.pdf new file mode 100644 index 0000000000..dcb6a578fb Binary files /dev/null and b/bmad/doc/bend-fiducial2.pdf differ diff --git a/bmad/doc/bend-vary1.pdf b/bmad/doc/bend-vary1.pdf new file mode 100644 index 0000000000..7e4d946374 Binary files /dev/null and b/bmad/doc/bend-vary1.pdf differ diff --git a/bmad/doc/bend-vary2.pdf b/bmad/doc/bend-vary2.pdf new file mode 100644 index 0000000000..07fb3b0343 Binary files /dev/null and b/bmad/doc/bend-vary2.pdf differ diff --git a/bmad/doc/charged-tracking.tex b/bmad/doc/charged-tracking.tex index 3244ff64b5..a940657fcb 100644 --- a/bmad/doc/charged-tracking.tex +++ b/bmad/doc/charged-tracking.tex @@ -368,7 +368,7 @@ \section{BeamBeam Tracking} \end{enumerate} There is an energy kick due to the motion of the strong beam. There are two parts to this $dp_z = -dp_{z,s} + dp_{z,h}$. One part, $dp_{z,s}$, is similar to the gravatational slingshot in orbital +dp_{z,s} + dp_{z,h}$. One part, $dp_{z,s}$, is similar to the gravitational slingshot in orbital mechanics. The slingshot energy kick is simply calculated using conservation of 4-momentum of the tracked particle and the strong beam where the mass of the strong beam is assumed to be large compared to the mass of the tracked particle.\footnote @@ -592,6 +592,146 @@ \section{Bend: Body Tracking with finite k1} m_{544} &= \frac{-1}{4 \, (1 + p_{z1})^2} \, (L + c_y \, s_y) \nonumber \end{alignat} +%--------------------------------------------------------------------------------- +%--------------------------------------------------------------------------------- +\section{Bend: Fiducial Point Calculations} +\label{s:bend.fiducial} + + +When the \vn{fiducial_pt} switch for a bend is set to something other than \vn{none}, changing one +of \vn{rho}, \vn{g}, \vn{b_field} or \vn{angle} in a program (that is, changing after the lattice +has been read in and the bend parameters calculated) involves adjustment to the other three +parameters along with adjustment to \vn{e1}, \vn{e2}, \vn{l}, \vn{l_chord}, and +\vn{l_rectangle}. This is done to keep the shape of the bend is invariant. Invariance is not +maintained with variation of any other parameter besides \vn{rho}, \vn{g}, \vn{b_field} and +\vn{angle}. For example, if \vn{e1} is varied the bend shape will vary. + +\begin{figure}[tb] + \centering + \hfill + \begin{subfigure}[b]{0.4\textwidth} + \includegraphics{bend-vary1.pdf} + \caption{ +With \vn{fiducial_pt} set to \vn{entrance_end}, $\bfr_1$ is the fiducial point at the entrance +end. By construction, the entrance point $\bfr_1$ and the +slope of the reference curve at $\bfr_1$ is invariant with the reference curve before (dashed line) +and after (solid line) being tangent to $\bfs_1$ where $\bfs_1$ being the perpendicular to $\bfx_1$. + } + \label{f:bend.fid1} + \end{subfigure} + \hfill + \begin{subfigure}[b]{0.4\textwidth} + \includegraphics{bend-vary2.pdf} + \caption{ +With \vn{fiducial_pt} set to \vn{center}, $\bfr_c$ is the fiducial point at the center. By +construction, the reference curve always goes through $\bfr_c$ and the tangent of the reference +curve at $\bfr_c$ is invariant.} + \label{f:bend.fid2} + \end{subfigure} + \hfill + \caption{ +Geometry with \vn{fiducial_pt} set to (a) \vn{entrance_end} and (b) \vn{center}. In both cases, +$\bfr_1$ and $\bfr_2$ are the entrance and exit reference points before and $\bfr'_1$ and $\bfr_2$ +are the entrance and exit points after variation of one of \vn{rho}, \vn{g}, \vn{b_field}, or +\vn{angle}. Similarly, $\rho$ and $\alpha$ are the bending radius and bending angle before +variation while $\rho'$ and $\alpha'$ are the bending radius afterwards. Finally, $e_1$ $e_2$ +are the face angles and rectangular length before variation, and $L'_r$ and $\bfr'_0$ are the +rectangular length and center of curvature after variation. + } + \label{f:bend.fid} +\end{figure} + +\fig{f:bend.fid} shows the situation when the \vn{fiducial_pt} is set to either \vn{entrance_end} or +\vn{center} (the situation for the \vn{exit_end} setting is analogous to the \vn{entrance_end} +setting and so is not discussed). For any one of the \vn{fiducial_pt} settings discussed there are +essentially two cases. One case is direct variation of the bend field via variation of \vn{rho}, +\vn{g}, or \vn{b_field}. This is called ``\vn{g}-variation''. The other type of variation is +variation of \vn{angle}. This is called ``angle-variation''. The discussion below shows how, with +\vn{g}-variation, \vn{l}, \vn{e1}, and \vn{e2} are calculated. With \vn{angle}-variation, \vn{l}, +\vn{g}, \vn{e1}, and \vn{e2} need to be calculated. Once \vn{l} and \vn{g} are know, the other +parameters \vn{l_chord}, \vn{l_rectangle}, \vn{l_sagitta} (and \vn{angle} for the \vn{g}-variation +case) can be readily computed. + +The \vn{entrance_end} analysis is as follows (\fig{f:bend.fid1}). The entrance end coordinates +around the point $\bfr_1$ are held fixed and as as a result $\bfr'_1 = \bfr_1$ and \vn{e1} does not +vary as well. $\bfr_2$ is the exit point before variation and $\bfr_3$ is the exit point after. The +position of $\bfr_3$ is calculated by first calculating the position of $\bfr_1$ in a coordinate +system centered at $\bfr_2$ and with axes parallel to the $(\bfs_1, \bfx_1)$ axes of the coordinate +system at $\bfr_1$ +\begin{equation} + \bar\bfr_1 = \left( -l_{\text{rectangle}}, \rho \, (1 - \cos\alpha) \right) +\end{equation} +Where the bar denotes that the coordinates are in the $(\bfs_1, \bfx_1)$ system. +The coordinates of $\bfr_1$ in the $(\bfe_s, \bfe_x)$ coordinate system +with origin at $\bfr_2$ and with $\bfe_x$ along the bend edge and $\bfe_s$ perpendicular to $\bfe_x$ +is a rotation $\bfR(\theta)$ +\begin{equation} + \bfr_1 = \bfR(\alpha - e_2) \, \bar\bfr_1 +\end{equation} +The angle $\theta_1$ of the vector $\bfs_1$, which is the invariant tangent of the reference curve +at the point $\bfr_1$, in the $(\bfe_s, \bfe_x)$ coordinate system (which is used from here on) is +\begin{equation} + \theta_1 = \alpha - e_1 +\end{equation} +The center of curvature after variation $\bfr'_0$ is +\begin{equation} + \bfr'_0 = \bfr_1 + \rho' \, \left( \sin\theta_1, -\cos\theta_1 \right) +\end{equation} +The reference trajectory after variation $\bfr'$ is a circular arc subject to the condition +\begin{equation} + \left| \bfr' - \bfr'_0 \right| = \rho'^2 + \label{rr0r} +\end{equation} +With \vn{g}-variation, the value of $\rho'$ is set (perhaps indirectly) by the User. To find the +point $\bfr'_2$, it is noted that in the $(e_s, e_x)$ coordinate system, the $s$ +coordinate of $\bfr'_2$, $r'_{2s}$ is zero. +so using this in\Eq{rr0r} and throwing away the unphysical root gives for the $x$ coordinate +\begin{equation} + r'_{2x} = r_{1x} + \frac{2 \, c}{-b - \sqrt{b^2-4 \, a \, c}} +\end{equation} +where +\begin{align} + a &= g' \CRNO + b &= 2 \, \cos\theta_1 \\ + c &= g' \, r_{1s}^2 + 2 \, r_{1s} \, \sin\theta_1 \nonumber +\end{align} +where $g' = 1 / \rho'$. +The rectangular length after variation $L'_r$ is then +\begin{equation} + L'_r = L_r + r'_{2x} * \sin\theta_1 +\end{equation} +where $L_r$ is the rectangular length before variation. +Finally, the length $L'$ after variation is +\begin{equation} + L' = \text{asinc} \left( g' \, L'_r \right) \, L'_r +\end{equation} +where $\text{asinc}$ is the function +\begin{equation} + \text{asinc}(\theta) = \frac{\sin^{-1}(\theta)}{\theta} +\end{equation} + +For \vn{fiducial_pt} set to \vn{entrance_end} and with \vn{angle}-variation, $\alpha'$ is know and +$g'$ can be computed via +\begin{equation} + g' = \frac{\sin(\alpha' - \theta_1) + \sin(\theta_1)}{r_{1s}} + \label{gatt} +\end{equation} +With this, all other parameters can be created. In both \vn{angle}-variation and \vn{g}-variation the new +face angle $e'_2$ is given by +\begin{equation} + e'_2 = e_2 + \alpha' - \alpha +\end{equation} + +For \vn{fiducial_pt} set to \vn{center}, The center point $\bfr_c$ (see \fig{f:bend.fid2}) is held +constant. Here the \vn{g}-variation analysis is similar to the \vn{g}-variation analysis with +\vn{fiducial_pt} set to \vn{entrance_end} (or \vn{exit_end} except in this case the reference orbit +to the right and left of $bfr_c$ are analyzed separately and the two lengths for each piece are +added together. For \vn{angle}-variation, the only situation where it is possible to keep $\bfr_c$ +fixed while varying the angle is when \vn{e1} and \vn{e2} are equal. In this instance, the +calculation is again similar to the \vn{angle}-variation analysis with the \vn{fiducial_pt} set to +either end. If \vn{e1} and \vn{e2} are not equal, a calculation is done that gives the desired angle +but the center point will shift. + %--------------------------------------------------------------------------------- %--------------------------------------------------------------------------------- \section{Converter Tracking} diff --git a/bmad/doc/elements.tex b/bmad/doc/elements.tex index 75795799dc..05f5eb1a8f 100644 --- a/bmad/doc/elements.tex +++ b/bmad/doc/elements.tex @@ -526,7 +526,7 @@ \section{Bends: Rbend and Sbend} \hfill \begin{subfigure}[b]{0.32\textwidth} \includegraphics{rbend-coords.pdf} - \caption{Rbend} + \caption{Rbend (with \vn{fiducial_pt} set to \vn{none} or \vn{center}).} \label{f:rbend} \end{subfigure} \hfill @@ -543,13 +543,18 @@ \section{Bends: Rbend and Sbend} \end{subfigure} \hfill \caption[Coordinate systems for rbend and sbend elements.] -{Coordinate systems for (a) normal (non-reversed)\vn{rbend}, (b) normal \vn{sbend}, and (c) + { +Coordinate systems for (a) normal (non-reversed)\vn{rbend}, (b) normal \vn{sbend}, and (c) \vn{reversed sbend} elements. The bends are viewed from ``above'' (viewed from positive $y$). Normal bends have \vn{g}, \vn{angle}, and \vn{rho} all positive. Reversed bends have \vn{g}, -\vn{angle}, and \vn{rho} all negative. For (a) and (b), as drawn, \vn{e1} and \vn{e2} are both -positive. For (c), as drawn, \vn{e1} and \vn{e2} are both negative. In all cases, \vn{L} is -positive. Notice that for reversed bends, the $x$-axis points towards the center of the bend while -for normal bends the $x$-axis points towards the outside.} +\vn{angle}, and \vn{rho} all negative. For (a) and (b), as drawn, the \vn{e1} and \vn{e2} face +angles are both positive. For (c), as drawn, \vn{e1} and \vn{e2} are both negative. In all cases, +\vn{L} is positive. Notice that for reversed bends, the $x$-axis points towards the center of the +bend while for normal bends the $x$-axis points towards the outside. Note: The geometry of the face +angles for \vn{rbend}s as shown in \fig{f:rbend} is modified if \vn{fiducial_pt} is set to +\vn{entrance_end} or \vn{exit_end}. See below for details. + + } \label{f:bend} \end{figure} @@ -575,8 +580,10 @@ \section{Bends: Rbend and Sbend} The rotation angle of the entrance pole face is \vn{e1} and at the exit face it is \vn{e2}. Zero \vn{e1} and \vn{e2} for an \vn{rbend} gives a rectangular magnet (\fig{f:rbend}). Zero \vn{e1} and \vn{e2} for an \vn{sbend} gives a wedge shaped magnet (\fig{f:sbend}). An \vn{sbend} with -an \vn{e1} = \vn{e2} = \vn{angle}/2 is equivalent to an \vn{rbend} with \vn{e1} = \vn{e2} = 0 (see -above). This formula holds for both positive and negative angles. +an \vn{e1} = \vn{e2} = \vn{angle}/2 is equivalent to an \vn{rbend} with \vn{e1} = \vn{e2} = 0. +This formula holds for both positive and negative angles. The geometry of the face +angles for \vn{rbend}s as shown in \fig{f:rbend} is modified if \vn{fiducial_pt} is set to +\vn{entrance_end} or \vn{exit_end}. See below for details. Note: The correspondence between \vn{e1} and \vn{e2} and the corresponding parameters used in the SAD program \cite{b:sad} is: @@ -756,7 +763,8 @@ \section{Bends: Rbend and Sbend} The attributes \vn{g}, \vn{angle}, and \vn{l} are mutually dependent. If any two are specified for an element \bmad will calculate the appropriate value for the third. After reading in a lattice, -\vn{angle} is considered a dependent variable (\sref{s:depend}). +\vn{angle} is considered the dependent variable so if \vn{l} or \vn{g} is veried, the value of +\vn{angle} will be set to \vn{g * l}. if \vn{theta} is varied, \vn{l} will be set accordingly. Since internally all \vn{rbend}s are converted to \vn{sbend}s, if one wants to vary the \vn{g} attribute of a bend and still keep the bend rectangular, an overlay (\sref{s:overlay}) can be diff --git a/bmad/modules/attribute_mod.f90 b/bmad/modules/attribute_mod.f90 index 91ad409006..e0a4ada5df 100644 --- a/bmad/modules/attribute_mod.f90 +++ b/bmad/modules/attribute_mod.f90 @@ -1487,7 +1487,7 @@ subroutine init_attribute_name_array () call init_attribute_name1 (rf_bend$, rho$, 'RHO', quasi_free$) call init_attribute_name1 (rf_bend$, l_chord$, 'L_CHORD', quasi_free$) call init_attribute_name1 (rf_bend$, l_sagitta$, 'L_SAGITTA', dependent$) -call init_attribute_name1 (rf_bend$, l_rectangle$, 'L_RECTANGLE', dependent$) +call init_attribute_name1 (rf_bend$, l_rectangle$, 'L_RECTANGLE', quasi_free$) call init_attribute_name1 (rf_bend$, fiducial_pt$, 'FIDUCIAL_PT') call init_attribute_name1 (rf_bend$, b_field$, 'B_FIELD', quasi_free$) call init_attribute_name1 (rf_bend$, init_needed$, 'init_needed', private$) @@ -1524,7 +1524,7 @@ subroutine init_attribute_name_array () call init_attribute_name1 (sbend$, init_needed$, 'init_needed', private$) call init_attribute_name1 (sbend$, l_chord$, 'L_CHORD', quasi_free$) call init_attribute_name1 (sbend$, l_sagitta$, 'L_SAGITTA', dependent$) -call init_attribute_name1 (sbend$, l_rectangle$, 'L_RECTANGLE', dependent$) +call init_attribute_name1 (sbend$, l_rectangle$, 'L_RECTANGLE', quasi_free$) call init_attribute_name1 (sbend$, fiducial_pt$, 'FIDUCIAL_PT') call init_attribute_name1 (sbend$, ptc_fringe_geometry$, 'PTC_FRINGE_GEOMETRY') call init_attribute_name1 (sbend$, b_field$, 'B_FIELD', quasi_free$) diff --git a/bmad/modules/changed_attribute_bookkeeper.f90 b/bmad/modules/changed_attribute_bookkeeper.f90 index 81294893db..442ec0366e 100644 --- a/bmad/modules/changed_attribute_bookkeeper.f90 +++ b/bmad/modules/changed_attribute_bookkeeper.f90 @@ -743,24 +743,7 @@ subroutine set_flags_for_changed_real_attribute (ele, attrib, set_dependent) ! is independent of the setting for %field_master. if (dep_set) then if (associated(a_ptr, ele%value(angle$))) then - if (ele%value(l$) /= 0) then - dangle = ele%value(angle$) - ele%value(g$) * ele%value(l$) - select case (nint(ele%value(fiducial_pt$))) - case (none_pt$) - ele%value(g$) = ele%value(angle$) / ele%value(l$) - case (center_pt$) - ele%value(g$) = ele%value(angle$) / ele%value(l$) - ele%value(e1$) = ele%value(e1$) + 0.5_rp * dangle - ele%value(e2$) = ele%value(e2$) + 0.5_rp * dangle - case (entrance_end$) - ele%value(g$) = sin(ele%value(angle$)) / ele%value(l_rectangle$) - ele%value(e2$) = ele%value(e2$) + dangle - case (exit_end$) - ele%value(g$) = sin(ele%value(angle$)) / ele%value(l_rectangle$) - ele%value(e1$) = ele%value(e1$) + dangle - end select - if (dep2_set) ele%value(b_field$) = ele%value(g$) * p0c_factor - endif + call bend_angle_fiducial_calc(ele, dep2_set, p0c_factor) elseif (associated(a_ptr, ele%value(l_rectangle$))) then select case (nint(ele%value(fiducial_pt$))) @@ -784,17 +767,17 @@ subroutine set_flags_for_changed_real_attribute (ele, attrib, set_dependent) elseif (associated(a_ptr, ele%value(rho$))) then if (ele%value(rho$) /= 0) then ele%value(g$) = 1.0_rp / ele%value(rho$) - call bend_fiducial_calc(ele, dep2_set, p0c_factor) + call bend_g_fiducial_calc(ele, dep2_set, p0c_factor) endif elseif (associated(a_ptr, ele%value(g$))) then ele%value(b_field$) = ele%value(g$) * p0c_factor - call bend_fiducial_calc(ele, dep2_set, p0c_factor) + call bend_g_fiducial_calc(ele, dep2_set, p0c_factor) elseif (dep2_set) then if (associated(a_ptr, ele%value(b_field$))) then ele%value(g$) = ele%value(b_field$) / p0c_factor - call bend_fiducial_calc(ele, dep2_set, p0c_factor) + call bend_g_fiducial_calc(ele, dep2_set, p0c_factor) elseif (associated(a_ptr, ele%value(db_field$))) then ele%value(dg$) = ele%value(db_field$) / p0c_factor @@ -862,31 +845,124 @@ subroutine set_flags_for_changed_real_attribute (ele, attrib, set_dependent) !-------------------------------------------------------------- contains -subroutine bend_fiducial_calc(ele, dep2_set, p0c_factor) +subroutine bend_g_fiducial_calc(ele, dep2_set, p0c_factor) + +type (ele_struct) ele +real(rp) p0c_factor, len1, len2 +logical dep2_set + +! There has been a change in g (and not any other parameters) so need to calculate shifts in L, e1, and e2 + +select case (nint(ele%value(fiducial_pt$))) +case (none_pt$) + +case (center_pt$) + call bend_g_fiducial_calc2(ele, 0.5_rp, ele%value(e1$), len1) + call bend_g_fiducial_calc2(ele, 0.5_rp, ele%value(e2$), len2) + ele%value(l$) = len1 + len2 + + +case (entrance_end$) + call bend_g_fiducial_calc2(ele, 1.0_rp, ele%value(e2$), ele%value(l$)) + +case (exit_end$) + call bend_g_fiducial_calc2(ele, 1.0_rp, ele%value(e1$), ele%value(l$)) +end select + +if (dep2_set) ele%value(b_field$) = ele%value(g$) * p0c_factor + +end subroutine bend_g_fiducial_calc + +!-------------------------------------------------------------- +! contains + +subroutine bend_g_fiducial_calc2(ele, factor, e_edge, len_out, angle0_in) + +type (ele_struct) ele +real(rp) factor, e_edge, len_out, g0, angle0, rp2x, r1(2), theta1, l_rect, length, a, b, c +real(rp), optional :: angle0_in + +! + +g0 = ele%value(l$) * ele%value(angle$) ! Old g +if (present(angle0_in)) then + angle0 = angle0_in ! Old angle +else + angle0 = factor * ele%value(angle$) ! Old angle +endif +theta1 = angle0 - e_edge + +r1 = factor * [-ele%value(l_rectangle$), -cos_one(angle0) * ele%value(l$)] +r1 = rot_2d(r1, theta1) + +a = ele%value(g$) +b = 2.0_rp * cos(theta1) +c = ele%value(g$) * r1(1)**2 + 2.0_rp * r1(1) * sin(theta1) +rp2x = r1(2) - 2.0_rp * c / (b + sqrt(b*b - 4.0_rp*a*c)) + +l_rect = factor * ele%value(l_rectangle$) + rp2x * sin(theta1) +length = asinc(ele%value(g$) * l_rect) * l_rect + +e_edge = e_edge + length * ele%value(g$) - angle0 +len_out = length + +end subroutine bend_g_fiducial_calc2 + +!-------------------------------------------------------------- +! contains + +subroutine bend_angle_fiducial_calc(ele, dep2_set, p0c_factor) type (ele_struct) ele -real(rp) p0c_factor, dangle +real(rp) p0c_factor, angle0, len1, len2 logical dep2_set -! There has been a change in g so need to calculate shifts in L, e1, and e2 +! There has been a change in angle (and not any other parameters) so need to calculate shifts in L, e1, and e2 -dangle = ele%value(g$) * ele%value(l$) - ele%value(angle$) +if (ele%value(l$) == 0) return select case (nint(ele%value(fiducial_pt$))) case (none_pt$) + ele%value(g$) = ele%value(angle$) / ele%value(l$) + case (center_pt$) - ele%value(e1$) = ele%value(e1$) + 0.5_rp * dangle - ele%value(e1$) = ele%value(e1$) + 0.5_rp * dangle + call bend_angle_fiducial_calc2(ele, 0.5_rp, ele%value(e2$), len1) + call bend_angle_fiducial_calc2(ele, 0.5_rp, ele%value(e2$), len2) + ele%value(l$) = len1 + len2 + ele%value(g$) = ele%value(angle$) / ele%value(l$) + case (entrance_end$) - ele%value(e2$) = ele%value(e2$) + dangle + call bend_angle_fiducial_calc2(ele, 1.0_rp, ele%value(e2$), ele%value(l$)) + case (exit_end$) - ele%value(e1$) = ele%value(e1$) + dangle + call bend_angle_fiducial_calc2(ele, 1.0_rp, ele%value(e1$), ele%value(l$)) end select if (dep2_set) ele%value(b_field$) = ele%value(g$) * p0c_factor +end subroutine bend_angle_fiducial_calc + +!-------------------------------------------------------------- +! contains + +subroutine bend_angle_fiducial_calc2(ele, factor, e_edge, len_out) + +type (ele_struct) ele +real(rp) factor, e_edge, angle0, theta1, r1(2), len_out + +! fiducial_pt = center is not considered here. + +angle0 = factor * ele%value(g$) * ele%value(l$) ! Old angle +theta1 = angle0 - e_edge + +r1 = factor * [-ele%value(l_rectangle$), -cos_one(angle0) * ele%value(l$)] +r1 = rot_2d(r1, theta1) + +ele%value(g$) = -(sin(theta1) + sin(ele%value(angle$) - theta1)) / r1(1) + +call bend_g_fiducial_calc2(ele, factor, e_edge, len_out, angle0) -end subroutine bend_fiducial_calc +end subroutine bend_angle_fiducial_calc2 end subroutine set_flags_for_changed_real_attribute diff --git a/bmad/output/type_ele.f90 b/bmad/output/type_ele.f90 index cf2e71a5c8..9e010c515f 100644 --- a/bmad/output/type_ele.f90 +++ b/bmad/output/type_ele.f90 @@ -1671,7 +1671,7 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c ! Temp until bend fiducial_pt code finished -if (ix2_attrib == fiducial_pt$ .or. ix2_attrib == l_rectangle$) ix2_attrib = -1 +!! if (ix2_attrib == fiducial_pt$ .or. ix2_attrib == l_rectangle$) ix2_attrib = -1 end function is_2nd_column_attribute diff --git a/bmad/output/write_bmad_lattice_file.f90 b/bmad/output/write_bmad_lattice_file.f90 index e5a92c726d..740be6184e 100644 --- a/bmad/output/write_bmad_lattice_file.f90 +++ b/bmad/output/write_bmad_lattice_file.f90 @@ -67,7 +67,7 @@ subroutine write_bmad_lattice_file (bmad_file, lat, err, output_form, orbit0) type (lat_ele_order_struct) order type (material_struct), pointer :: material -real(rp) s0, x_lim, y_lim, val, x, y +real(rp) s0, x_lim, y_lim, val, x, y, fid character(*) bmad_file character(4000) line @@ -774,8 +774,9 @@ subroutine write_bmad_lattice_file (bmad_file, lat, err, output_form, orbit0) if (y_lim /=0 .and. ele%value(y2_limit$) == y_lim) y_lim_good = .true. !---------------------------------------------------------------------------- - ! Print the element attributes. + ! Write the element attributes. + fid = nint(ele%value(fiducial_pt$)) attribute_loop: do j = 1, num_ele_attrib$ attrib = attribute_info(ele, j) val = ele%value(j) @@ -783,6 +784,9 @@ subroutine write_bmad_lattice_file (bmad_file, lat, err, output_form, orbit0) if (ele%key == sbend$) then if (j == fintx$ .and. ele%value(fintx$) == ele%value(fint$)) cycle if (j == hgapx$ .and. ele%value(hgapx$) == ele%value(hgap$)) cycle + + if (j == l$ .and. (fid == entrance_end$ .or. fid == entrance_end$)) cycle + if (j == l_rectangle$ .and. (fid == none_pt$ .or. fid == center_pt$)) cycle endif if (j == check_sum$) cycle if (x_lim_good .and. (j == x1_limit$ .or. j == x2_limit$)) cycle diff --git a/bmad/parsing/bmad_parser_mod.f90 b/bmad/parsing/bmad_parser_mod.f90 index 46c4680674..622566e886 100644 --- a/bmad/parsing/bmad_parser_mod.f90 +++ b/bmad/parsing/bmad_parser_mod.f90 @@ -1972,6 +1972,10 @@ subroutine parser_set_attribute (how, ele, delim, delim_found, err_flag, pele, c call get_switch (attrib_word, exact_multipoles_name(1:), ix, err_flag, ele, delim, delim_found); if (err_flag) return ele%value(exact_multipoles$) = ix +case ('FIDUCIAL_PT') + call get_switch (attrib_word, fiducial_pt_name(1:), ix, err_flag, ele, delim, delim_found); if (err_flag) return + ele%value(fiducial_pt$) = ix + case ('FIELD_CALC') call get_switch (attrib_word, field_calc_name(1:), ele%field_calc, err_flag, ele, delim, delim_found); if (err_flag) return @@ -7440,23 +7444,15 @@ subroutine settable_dep_var_bookkeeping (ele) if (ele%key == rbend$) then select case (nint(ele%value(fiducial_pt$))) case (none_pt$, center_pt$) - if (ele%value(e1$) == real_garbage$) ele%value(e1$) = 0 - if (ele%value(e2$) == real_garbage$) ele%value(e2$) = 0 + ele%value(e1$) = ele%value(e1$) + 0.5_rp * ele%value(angle$) + ele%value(e2$) = ele%value(e2$) + 0.5_rp * ele%value(angle$) case (entrance_end$) - if (ele%value(e1$) == real_garbage$) ele%value(e1$) = -0.5_rp * ele%value(angle$) - if (ele%value(e2$) == real_garbage$) ele%value(e2$) = 0.5_rp * ele%value(angle$) + ele%value(e2$) = ele%value(e2$) + ele%value(angle$) case (exit_end$) - if (ele%value(e1$) == real_garbage$) ele%value(e1$) = 0.5_rp * ele%value(angle$) - if (ele%value(e2$) == real_garbage$) ele%value(e2$) = -0.5_rp * ele%value(angle$) + ele%value(e1$) = ele%value(e1$) + ele%value(angle$) end select - ele%value(e1$) = ele%value(e1$) + 0.5_rp * ele%value(angle$) - ele%value(e2$) = ele%value(e2$) + 0.5_rp * ele%value(angle$) ele%key = sbend$ - - else - if (ele%value(e1$) == real_garbage$) ele%value(e1$) = 0 - if (ele%value(e2$) == real_garbage$) ele%value(e2$) = 0 endif ! diff --git a/bmad/parsing/set_ele_defaults.f90 b/bmad/parsing/set_ele_defaults.f90 index c826bc8536..f9bb684f00 100644 --- a/bmad/parsing/set_ele_defaults.f90 +++ b/bmad/parsing/set_ele_defaults.f90 @@ -292,8 +292,6 @@ subroutine set_ele_defaults (ele, do_allocate) ele%value(ptc_field_geometry$) = sector$ ele%value(fiducial_pt$) = none_pt$ ele%value(init_needed$) = true$ - ele%value(e1$) = real_garbage$ - ele%value(e2$) = real_garbage$ case (rcollimator$) ele%offset_moves_aperture = .true. diff --git a/regression_tests/CMakeLists.txt b/regression_tests/CMakeLists.txt index 80af69177a..f24b527b90 100644 --- a/regression_tests/CMakeLists.txt +++ b/regression_tests/CMakeLists.txt @@ -8,6 +8,7 @@ set (EXE_SPECS cmake_files/cmake.autoscale_test cmake_files/cmake.backwards_time_track_test cmake_files/cmake.beam_test + cmake_files/cmake.bend_test cmake_files/cmake.bbu_test cmake_files/cmake.cesr_test cmake_files/cmake.closed_orbit_test diff --git a/regression_tests/bend_test/bend_test.bmad b/regression_tests/bend_test/bend_test.bmad new file mode 100644 index 0000000000..a606234b16 --- /dev/null +++ b/regression_tests/bend_test/bend_test.bmad @@ -0,0 +1,20 @@ +f = 1 +beginning[e_tot] = 1e6*f +parameter[no_end_marker] = T + +parameter[particle] = electron +parameter[geometry] = open + +beginning[beta_a] = 10 +beginning[beta_b] = 10 +beginning[alpha_a] = 0 +beginning[alpha_b] = 0 +beginning[etap_x] = 0.001 + +b: rbend, l_rectangle = 2.0, angle = -0.1, fiducial_pt = entrance_end, e2 = -0 + +ln: line = (b) +use, ln + +no_digested + diff --git a/regression_tests/bend_test/bend_test.f90 b/regression_tests/bend_test/bend_test.f90 new file mode 100644 index 0000000000..ea08fc9bf5 --- /dev/null +++ b/regression_tests/bend_test/bend_test.f90 @@ -0,0 +1,47 @@ +program bend_test + +use bmad + +implicit none + +type (lat_struct), target :: lat +type (ele_struct), pointer :: ele +type (ele_struct) ele0 + +real(rp) dparam(5), pa(4), f +integer ie, ip, iz, ix_param(5), np + +! + +open (1, file = 'output.now', recl = 200) + +bmad_com%auto_bookkeeper = .false. +call bmad_parser('bend_test.bmad', lat) + +f = lat%ele(0)%value(p0c$) / c_light +dparam = [0.01_rp, 0.1_rp, 0.1_rp * f, 0.1_rp, 0.1_rp] +ix_param = [g$, angle$, b_field$, l$, l_rectangle$] + + +do ie = 1, lat%n_ele_track + ele0 = lat%ele(ie) + + do ip = 1, size(ix_param) + ele => lat%ele(ie) + ele = ele0 + np = ix_param(ip) + ele%value(np) = ele%value(np) + dparam(ip) + call set_flags_for_changed_attribute(ele, ele%value(np)) + call lattice_bookkeeper(lat) + + pa(1:3) = [ele%value(ix_param(1)), ele%value(ix_param(2)), ele%value(ix_param(3)) / f] + write (1, '(a, i0, 3a, 4f16.10)') '"', ie, '-', trim(attribute_name(ele, np)), ':g-ang-field" ABS 1E-10', pa(1:3) + + pa = [ele%value(ix_param(4)), ele%value(ix_param(5)), ele%value(e1$), ele%value(e2$)] + write (1, '(a, i0, 3a, 4f16.10)') '"', ie, '-', trim(attribute_name(ele, np)), ':l-lr-e1/2" ABS 1E-10', pa + enddo +enddo + +close (1) + +end program diff --git a/regression_tests/bend_test/output.correct b/regression_tests/bend_test/output.correct new file mode 100644 index 0000000000..9b56da4fc9 --- /dev/null +++ b/regression_tests/bend_test/output.correct @@ -0,0 +1,20 @@ +"1-p1-G" ABS 1E-10 0.0100000000 0.0200013336 -0.0100000000 +"1-p2-G" ABS 1E-10 2.0001333573 2.0000000000 0.0000000000 0.0200013336 +"1-p1-ANGLE" ABS 1E-10 0.0499167083 0.0998334166 -0.0499167083 +"1-p2-ANGLE" ABS 1E-10 2.0000000000 1.9966794182 0.0000000000 0.1000000000 +"1-p1-B_FIELD" ABS 1E-10 -0.1000000000 -0.2013579208 0.1000000000 +"1-p2-B_FIELD" ABS 1E-10 2.0135792079 2.0000000000 0.0000000000 -0.2013579208 +"1-p1-L" ABS 1E-10 0.0000000000 0.0000000000 0.0000000000 +"1-p2-L" ABS 1E-10 2.1000000000 2.1000000000 0.0000000000 0.0000000000 +"1-p1-L_RECTANGLE" ABS 1E-10 0.0000000000 0.0000000000 0.0000000000 +"1-p2-L_RECTANGLE" ABS 1E-10 2.1000000000 2.1000000000 0.0000000000 0.0000000000 +"2-p1-NOISE" ABS 1E-10 0.0100000000 0.0000000000 0.0000000000 +"2-p2-NOISE" ABS 1E-10 0.0000000000 0.0000000000 0.0000000000 0.0000000000 +"2-p1-!NULL" ABS 1E-10 0.0000000000 0.1000000000 0.0000000000 +"2-p2-!NULL" ABS 1E-10 0.0000000000 0.0000000000 0.0000000000 0.0000000000 +"2-p1-!NULL" ABS 1E-10 0.0000000000 0.0000000000 0.1000000000 +"2-p2-!NULL" ABS 1E-10 0.0000000000 0.0000000000 0.0000000000 0.0000000000 +"2-p1-L" ABS 1E-10 0.0000000000 0.0000000000 0.0000000000 +"2-p2-L" ABS 1E-10 0.1000000000 0.0000000000 0.0000000000 0.0000000000 +"2-p1-X_DISPERSION_ERR" ABS 1E-10 0.0000000000 0.0000000000 0.0000000000 +"2-p2-X_DISPERSION_ERR" ABS 1E-10 0.0000000000 0.1000000000 0.0000000000 0.0000000000 diff --git a/regression_tests/cmake_files/cmake.bend_test b/regression_tests/cmake_files/cmake.bend_test new file mode 100644 index 0000000000..80f2fdd54a --- /dev/null +++ b/regression_tests/cmake_files/cmake.bend_test @@ -0,0 +1,12 @@ +set (EXENAME bend_test) + +FILE (GLOB SRC_FILES "bend_test/*.f90") + +set (INC_DIRS +) + +set (LINK_LIBS + bmad + sim_utils + ${ACC_BMAD_LINK_LIBS} +) \ No newline at end of file diff --git a/sim_utils/interfaces/sim_utils_interface.f90 b/sim_utils/interfaces/sim_utils_interface.f90 index f57a7e0110..2337c1b780 100644 --- a/sim_utils/interfaces/sim_utils_interface.f90 +++ b/sim_utils/interfaces/sim_utils_interface.f90 @@ -665,6 +665,15 @@ function probability_funct(x) result (prob) real(rp) x end function +function quadratic_roots(coefs) result (root) + +use precision_def, only: rp + import + implicit none + real(rp) coefs(3) + complex(rp) root(2) +end function + subroutine query_string (query_str, upcase, return_str, ix, ios) implicit none character(*) return_str diff --git a/sim_utils/math/quadratic_roots.f90 b/sim_utils/math/quadratic_roots.f90 new file mode 100644 index 0000000000..a6cd8aa4ac --- /dev/null +++ b/sim_utils/math/quadratic_roots.f90 @@ -0,0 +1,62 @@ +!+ +! Function quadratic_roots(coefs) result (root) +! +! Routine to return the roots of the quadratic equation to machine precision. +! That is, avoiding the round-off error inherent in the standard formula. +! +! If coefs(3) is zero 0, root(2) is set to real_garbage$. +! If both coefs(2) and coefs(3) are zero, both roots are set to real_garbage$. +! +! When the roots are complex, the roots are ordered so that Imag(root(1)) < Imag(root(2)). +! When the roots are real, the roots are ordered so that root(1) <= root(2). +! +! Input: +! coefs(:) -- real(rp): Coefficients of the quadratic equation with +! 0 = coefs(1) + coefs(2) * x + coefs(3) * x^2 +! +! Output: +! root(2) -- complex(rp): Complex roots. +!- + +function quadratic_roots(coefs) result (root) + +use sim_utils_struct, only: rp, real_garbage$ +implicit none + +! + +real(rp) coefs(3), rad, den +complex(rp) root(2) + +! + +if (coefs(3) == 0) then + if (coefs(2) == 0) then + root = real_garbage$ + else + root(1) = -coefs(1) / coefs(2) + root(2) = real_garbage$ + endif + return +endif + +rad = coefs(2)**2 - 4.0_rp * coefs(1) * coefs(3) +den = 0.5_rp / coefs(3) + +if (rad < 0) then + rad = sqrt(-rad) + root(1) = den * cmplx(-coefs(2), -rad) + root(2) = den * cmplx(-coefs(2), rad) + +elseif (coefs(2) > 0) then + rad = sqrt(rad) + root(1) = den * (-coefs(2) - rad) + root(2) = 2.0_rp * coefs(1) / (-coefs(2) - rad) + +else + rad = sqrt(rad) + root(1) = den * (-coefs(2) + rad) + root(2) = 2.0_rp * coefs(1) / (-coefs(2) + rad) +endif + +end function diff --git a/tao/code/tao_set_mod.f90 b/tao/code/tao_set_mod.f90 index bdfacbc77a..ba5f023035 100644 --- a/tao/code/tao_set_mod.f90 +++ b/tao/code/tao_set_mod.f90 @@ -3428,7 +3428,7 @@ subroutine tao_set_drawing_cmd (drawing, component, value_str) case ('ele_name', 'name') needs_quotes = .true. component = 'ele_id%' // component(ix+1:) - case ('shape', 'color', 'label') + case ('shape', 'color', 'label', 'ele_id') needs_quotes = .true. end select diff --git a/tao/code/tao_show_this.f90 b/tao/code/tao_show_this.f90 index b23cea7e17..402a947bdb 100644 --- a/tao/code/tao_show_this.f90 +++ b/tao/code/tao_show_this.f90 @@ -4059,8 +4059,8 @@ subroutine tao_show_this (what, result_id, lines, nl) nl=nl+1; lines(nl) = '' nl=nl+1; lines(nl) = 'Element Shapes:' - nl=nl+1; lines(nl) = ' Shape Type Shape Multi Line_' - nl=nl+1; lines(nl) = ' Ele_ID Shape Color Size Label Draw Shape Width Offset' + nl=nl+1; lines(nl) = ' Line_' + nl=nl+1; lines(nl) = ' Ele_ID Shape Color Size Label Draw Multi Width Offset' nl=nl+1; lines(nl) = ' ------------------------------ ---------- ------- ---- ------ ----- ----- ----- ------' do i = 1, size(shapes) diff --git a/tao/version/tao_version_mod.f90 b/tao/version/tao_version_mod.f90 index 302d673bca..3d2f9d5705 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/13 12:55:46" +character(*), parameter :: tao_version_date = "2024/06/16 16:24:57" end module