Fix case where radiation setup fails due to extreme instability in th…

…e orbit.

bmad-doc/tutorial_bmad_tao/doc/tutorial_bmad_tao.tex

@@ -79,7 +79,7 @@
 
 \title{An Introduction and Tutorial to Bmad and Tao}
 \author{}
-\date{David Sagan and Chris Mayes \\ June 7, 2024}
+\date{David Sagan and Chris Mayes \\ July 18, 2024}
 
 \begin{document}
 
@@ -159,17 +159,21 @@ \subsection{How to use this Tutorial!}
 \Section{Orientation}
 \label{s:orientation}
 
+%------------------------------------------------------------------------------
+\subsection{Web Site}
+
 The \bmad web site is at:
 \begin{display}
  \url{\detokenize{www.classe.cornell.edu/bmad/}}
 \end{display}
+Here you will find documentation on the \bmad toolkit and associated programs as well as
+installation instructions.
 
 %------------------------------------------------------------------------------
 \subsection{Prerequisites}
 
-\bmad and \tao are generally used with Linux or Mac OS-X. Historically, while \bmad could be built
-under Windows, running \tao was problematical due to graphics issues. As present, these graphics issues have
-been resolved.
+\bmad and \tao can be run under Linux, Mac OS-X, or Windows running Windows Subsystem for Linux 2 (WSL2)
+which allows running in a Linux environment on Windows.
 
 Except when using a Python interface, \tao is accessed through the command line. Therefore you will
 need to run a terminal program to be able to run \tao.
@@ -236,12 +240,8 @@ \subsection{Tao Executable Installation}
 \subsection{Releases}
 \label{s:dist}
 
-A (GitHub) \vn{Release} is a set of files, including \bmad and \tao source files, which are used to
-build the \bmad, the \tao program, and various other simulation programs. A \vn{CLASSE Release} is like a
-GitHub \vn{Release} except that it is created on the Linux computer system at CLASSE (Cornell's
-Laboratory for Accelerator-based Sciences and Education). More information can be obtained from the
-\bmad web site. For the purpose of this tutorial, it is assumed that you are dealing with a
-GitHub \vn{Release}.
+A \vn{Release}, hosted on GitHub, is a set of files, including \bmad and \tao source files, which
+are used to build the \bmad, the \tao program, and various other simulation programs.
 
 The directories in a Release are:
 \begin{lstlisting}[mathescape]
@@ -3639,7 +3639,7 @@ \subsection{Reference Energy and the Lcavity and RFcavity Elements}
 
 For lattices with an \vn{open} geometry (\sref{s:lat.geom}), the reference energy/momentum at the
 beginning of the lattice is what is set in the lattice file. The reference energy/momentum for a
-downsteam element inherits the reference energy/momentum of the element just upstream. The exception
+downstream element inherits the reference energy/momentum of the element just upstream. The exception
 is for \vn{Lcavity} elements, which represent an RF cavity, where the reference energy/momentum at
 the exit end of the \vn{lcavity} is set so that a particle entering the cavity with zero phase space
 coordinates leaves with zero phase space coordinates and in particular phase space \vn{pz} will be
@@ -3728,16 +3728,16 @@ \subsection{Exercises [Answers in Section \sref{s:ans.phase}]}
 \item 
 Modify the \vn{lcavity} in the \vn{cavity.bmad} to have a small length (so the transit time is
 small), and set the beginning momentum small enough so the relativistic beta is significantly than
-one. Starting the particle with a finite $z$, calculate the ending $z$ after the cavity and verity
-that the cange in $z$ is consistent with the equation for $z$ given in \sref{s:phase.space.sub}.
+one. Starting the particle with a finite $z$, calculate the ending $z$ after the cavity and verify
+that the range in $z$ is consistent with the equation for $z$ given in \sref{s:phase.space.sub}.
 %
 \item
 \vn{Lcavity} elements have an attribute \vn{phi0_err} which varies the RF phase that a particle sees
 but does not change the reference energy. Add a finite \vn{phi0_err} to the \vn{cavity} element and
 verify that the reference energy does not change but that the phase space \vn{pz} of the particle
 (which is the normalized momentum deviation from the reference \sref{s:phase.space.sub}) does
 change. With the \vn{cavity.bmad} lattice, there is a range of values for \vn{phi0_err} where the
-cavity will decellerate the particle enough so that the particle will turn around and not make
+cavity will decelerate the particle enough so that the particle will turn around and not make
 it through the cavity. Approximately what is this range?
 \end{enumerate}
 
@@ -3981,7 +3981,7 @@ \subsection{Exercises [Answers in Section \sref{s:ans.super}]}
 
 \begin{enumerate}[label=\thesection.\arabic{enumi}]
 \item
-Create a lattice with a sextupole superimpsed in the middle of a quadrupole. Make the length of the 
+Create a lattice with a sextupole superimposed in the middle of a quadrupole. Make the length of the 
 sextupole shorter than the length of the quadrupole. Make a second lattice similar to the first but
 this time use \vn{jumbo} superposition. Compare both lattices.
 \end{enumerate}
@@ -5326,7 +5326,7 @@ \subsection{Exercises}
  = actual_measured + (design_orbit - model_orbit)
  = actual_measured + (design_orbit - 
  (actual_measured - (golden_orbit - design_orbit)))
- = gloden_orbit
+ = golden_orbit
 \end{code}
 %
 \item

bmad/code/chrom_calc.f90

@@ -58,8 +58,6 @@ subroutine chrom_calc (lat, delta_e, chrom_x, chrom_y, err_flag, &
 
 ! Init setup
 
-call cpu_time(time0)
-
 ix_br = integer_option(0, ix_branch)
 branch => lat%branch(ix_br)
 ele => branch%ele(0)
@@ -183,7 +181,5 @@ subroutine chrom_calc (lat, delta_e, chrom_x, chrom_y, err_flag, &
 
 if (present(err_flag)) err_flag = .false.
 
-call cpu_time(time1)
-if (bmad_com%debug) print '(a, f12.2)', 'chrom_calc execution time:', time1 - time0
 
 end subroutine

bmad/code/radiation_integrals.f90

@@ -108,7 +108,6 @@ subroutine radiation_integrals (lat, orbit, mode, ix_cache, ix_branch, rad_int_b
 real(rp) v(4,4), v_inv(4,4), z_here, z_start, mc2, gamma, gamma4, gamma6
 real(rp) kz, fac, c, s, factor, g2, g_x0, dz_small, z1, const_q, vec0(6), mat6(6,6)
 real(rp), parameter :: const_q_factor = 55 * h_bar_planck * c_light / (32 * sqrt_3) ! Cf: Sands Eq 5.46 pg 124.
-real time0, time1
 
 integer, optional :: ix_cache, ix_branch
 integer i, j, k, n, ix, ixe, ib, ip, ir, key2, n_step, ix_pole_max
@@ -121,8 +120,6 @@ subroutine radiation_integrals (lat, orbit, mode, ix_cache, ix_branch, rad_int_b
 !---------------------------------------------------------------------
 ! Allocate rad_int_by_ele
 
-call cpu_time(time0)
-
 if (present(rad_int_by_ele)) then
  if (allocated(rad_int_by_ele%branch)) then
  if (ubound(rad_int_by_ele%branch, 1) /= ubound(lat%branch, 1)) deallocate (rad_int_by_ele%branch)
@@ -682,9 +679,6 @@ subroutine radiation_integrals (lat, orbit, mode, ix_cache, ix_branch, rad_int_b
  enddo
 endif
 
-call cpu_time(time1)
-if (bmad_com%debug) print '(a, f12.2)', 'radiation_integrals execution time:', time1 - time0
-
 !------------------------------------------------------------------------
 contains
 

bmad/code/tracking_rad_map_setup.f90 → bmad/low_level/tracking_rad_map_setup.f90

@@ -1,5 +1,5 @@
 !+
-! Subroutine tracking_rad_map_setup (ele, tollerance, ref_edge, rad_map)
+! Subroutine tracking_rad_map_setup (ele, tollerance, ref_edge, rad_map, err_flag)
 !
 ! Routine to setup the radiation damping and excitation matrices for an element.
 ! These matrices can then be used when tracking to simulate radiation effects.
@@ -13,9 +13,10 @@
 !
 ! Output:
 ! rad_map -- rad_map_struct: Structure holding the matrices.
+! err_flag -- logical: Set True if there is an error. False otherwise.
 !-
 
-subroutine tracking_rad_map_setup (ele, tollerance, ref_edge, rad_map)
+subroutine tracking_rad_map_setup (ele, tollerance, ref_edge, rad_map, err_flag)
 
 use rad_6d_mod, dummy => tracking_rad_map_setup
 use f95_lapack, only: dpotrf_f95
@@ -33,10 +34,12 @@ subroutine tracking_rad_map_setup (ele, tollerance, ref_edge, rad_map)
 integer ref_edge
 integer i, info
 
-logical err
+logical err_flag, err
 
 !
 
+err_flag = .true.
+
 branch => pointer_to_branch(ele)
 tol = tollerance / branch%param%total_length
 orb0 = ele%map_ref_orb_in
@@ -49,7 +52,8 @@ subroutine tracking_rad_map_setup (ele, tollerance, ref_edge, rad_map)
 !
 
 call rad1_damp_and_stoc_mats (ele, .true., orb0, orb1, rad_map, &
- tol*branch%param%g2_integral, tol*branch%param%g3_integral)
+ tol*branch%param%g2_integral, tol*branch%param%g3_integral, err)
+if (err) return
 
 select case (ref_edge)
 case (upstream_end$)
@@ -67,12 +71,15 @@ subroutine tracking_rad_map_setup (ele, tollerance, ref_edge, rad_map)
 
 call dpotrf_f95 (rad_map%stoc_mat, 'L', info = info)
 if (info /= 0) then
- rad_map%stoc_mat = 0 ! Cholesky failed
+ rad_map%stoc_mat = 0 ! Cholesky failed.
+ err_flag = .false. ! But this is OK.
  return
 endif
 
 do i = 2, 6
  rad_map%stoc_mat(1:i-1, i) = 0
 enddo
 
+err_flag = .false.
+
 end subroutine

bmad/modules/bmad_routine_interface.f90

@@ -2543,6 +2543,13 @@ subroutine tilt_mat6 (mat6, tilt)
  real(rp) tilt, mat6(6,6)
 end subroutine
 
+subroutine to_surface_coords (lab_orbit, ele, surface_orbit)
+ import
+ implicit none
+ type (coord_struct) lab_orbit, surface_orbit
+ type (ele_struct) ele
+end subroutine
+
 subroutine track_a_beambeam (orbit, ele, param, track, mat6, make_matrix)
  import
  implicit none
@@ -2912,13 +2919,14 @@ subroutine track1_time_runge_kutta (orbit, ele, param, err_flag, track, t_end, d
  type (track_struct), optional :: track
 end subroutine
 
-subroutine tracking_rad_map_setup (ele, tollerance, ref_edge, rad_map)
+subroutine tracking_rad_map_setup (ele, tollerance, ref_edge, rad_map, err_flag)
  import
  implicit none
  type (ele_struct), target :: ele
  type (rad_map_struct) rad_map
  real(rp) tollerance
  integer ref_edge
+ logical err_flag
 end subroutine
 
 subroutine transfer_ac_kick (ac_kick_in, ac_kick_out)

bmad/modules/changed_attribute_bookkeeper.f90

@@ -743,7 +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
- call bend_angle_fiducial_calc(ele, dep2_set, p0c_factor)
+ call bend_angle_fiducial(ele, dep2_set, p0c_factor)
 
  elseif (associated(a_ptr, ele%value(l_rectangle$))) then
  select case (nint(ele%value(fiducial_pt$)))
@@ -767,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_g_fiducial_calc(ele, dep2_set, p0c_factor)
+ call bend_g_fiducial(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_g_fiducial_calc(ele, dep2_set, p0c_factor)
+ call bend_g_fiducial(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_g_fiducial_calc(ele, dep2_set, p0c_factor)
+ call bend_g_fiducial(ele, dep2_set, p0c_factor)
 
  elseif (associated(a_ptr, ele%value(db_field$))) then
  ele%value(dg$) = ele%value(db_field$) / p0c_factor
@@ -845,7 +845,7 @@ subroutine set_flags_for_changed_real_attribute (ele, attrib, set_dependent)
 !--------------------------------------------------------------
 contains
 
-subroutine bend_g_fiducial_calc(ele, dep2_set, p0c_factor)
+subroutine bend_g_fiducial(ele, dep2_set, p0c_factor)
 
 type (ele_struct) ele
 real(rp) p0c_factor, len1, len2
@@ -857,26 +857,26 @@ subroutine bend_g_fiducial_calc(ele, dep2_set, p0c_factor)
 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)
+ call bend_g_fiducial2(ele, 0.5_rp, ele%value(e1$), len1)
+ call bend_g_fiducial2(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$))
+ call bend_g_fiducial2(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$))
+ call bend_g_fiducial2(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
+end subroutine bend_g_fiducial
 
 !--------------------------------------------------------------
 ! contains
 
-subroutine bend_g_fiducial_calc2(ele, factor, e_edge, len_out, angle0_in)
+subroutine bend_g_fiducial2(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
@@ -906,12 +906,12 @@ subroutine bend_g_fiducial_calc2(ele, factor, e_edge, len_out, angle0_in)
 e_edge = e_edge + length * ele%value(g$) - angle0
 len_out = length
 
-end subroutine bend_g_fiducial_calc2
+end subroutine bend_g_fiducial2
 
 !--------------------------------------------------------------
 ! contains
 
-subroutine bend_angle_fiducial_calc(ele, dep2_set, p0c_factor)
+subroutine bend_angle_fiducial(ele, dep2_set, p0c_factor)
 
 type (ele_struct) ele
 real(rp) p0c_factor, angle0, len1, len2
@@ -926,26 +926,26 @@ subroutine bend_angle_fiducial_calc(ele, dep2_set, p0c_factor)
  ele%value(g$) = ele%value(angle$) / ele%value(l$)
 
 case (center_pt$)
- 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)
+ call bend_angle_fiducial2(ele, 0.5_rp, ele%value(e2$), len1)
+ call bend_angle_fiducial2(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$)
- call bend_angle_fiducial_calc2(ele, 1.0_rp, ele%value(e2$), ele%value(l$))
+ call bend_angle_fiducial2(ele, 1.0_rp, ele%value(e2$), ele%value(l$))
 
 case (exit_end$)
- call bend_angle_fiducial_calc2(ele, 1.0_rp, ele%value(e1$), ele%value(l$))
+ call bend_angle_fiducial2(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
+end subroutine bend_angle_fiducial
 
 !--------------------------------------------------------------
 ! contains
 
-subroutine bend_angle_fiducial_calc2(ele, factor, e_edge, len_out)
+subroutine bend_angle_fiducial2(ele, factor, e_edge, len_out)
 
 type (ele_struct) ele
 real(rp) factor, e_edge, angle0, theta1, r1(2), len_out
@@ -960,9 +960,9 @@ subroutine bend_angle_fiducial_calc2(ele, factor, e_edge, len_out)
 
 ele%value(g$) = -(sin(theta1) + sin(ele%value(angle$) - theta1)) / r1(1)
 
-call bend_g_fiducial_calc2(ele, factor, e_edge, len_out, angle0)
+call bend_g_fiducial2(ele, factor, e_edge, len_out, angle0)
 
-end subroutine bend_angle_fiducial_calc2
+end subroutine bend_angle_fiducial2
 
 end subroutine set_flags_for_changed_real_attribute