Fix Tao handling of mass/charge in expressions. (#1227)

bmad/code/closed_orbit_calc.f90

@@ -156,6 +156,7 @@ subroutine closed_orbit_calc (lat, closed_orb, i_dim, direction, ix_branch, err_
 bmad_com%aperture_limit_on = .false.
 bmad_com%spin_tracking_on = .false.
 bmad_com%spin_sokolov_ternov_flipping_on = .false.
+bmad_private%random_on = .false.
 
 n_ele = branch%n_ele_track
 betas = 1
@@ -580,6 +581,7 @@ subroutine end_cleanup(branch, reset_orb)
 !
 
 bmad_com = bmad_com_saved  ! Restore
+bmad_private%random_on = .true.
 
 if (n_dim == 4 .or. n_dim == 5) then
   call set_on_off (rfcavity$, branch%lat, restore_state$, ix_branch = branch%ix_branch, saved_values = on_off_state)

bmad/code/make_mat6.f90

@@ -91,6 +91,8 @@ recursive subroutine make_mat6 (ele, param, start_orb, end_orb, err_flag)
     if (global_com%exit_on_error) call err_exit
     return
   end select
+
+  if (ele%key == foil$) mat6_calc_method = tracking$
 endif
 
 ele%map_ref_orb_in = a_start_orb

bmad/doc/charged-tracking.tex

@@ -165,13 +165,13 @@ \section{Hamiltonian}
 \begin{equation}
   \wt m = \frac{m \, c^2}{c \, P_0} \qquad
   \left( a_x, a_y, \frac{a_s}{1+g \, x} \right) = \frac{q \, \Bf A}{P_0 \, c} \qquad 
-  \wt\psi(x,y,z) = \frac{q \, \psi}{P_0 \, c}
+  \wt\psi(x,y,z) = \frac{q \, \psi}{P_0}
   \label{mmccp}
 \end{equation}
 In terms of the normalized velocities $\beta_x$, $\beta_y$, the canonical momentum are
 \begin{equation}
-  p_x = \frac{m \, c^2}{P_0 \, c} \, \beta_x + a_x, \qquad 
-  p_y = \frac{m \, c^2}{P_0 \, c} \, \beta_y + a_y
+  p_x = \frac{m \, c^2}{P_0 \, c} \, \gamma \, \beta_x + a_x, \qquad 
+  p_y = \frac{m \, c^2}{P_0 \, c} \, \gamma \, \beta_y + a_y
   \label{pmc2pc}
 \end{equation}
 
@@ -204,7 +204,7 @@ \section{Hamiltonian}
 (\sref{s:phase.space}), \Eq{h1gx1} becomes
 \begin{equation}
   H = \frac{(p_x - a_x)^2}{2 (1 + p_z)} + \frac{(p_y - a_y)^2}{2 (1 + p_z)} - 
-  (1 + g \, x) \, (1 + p_z) - a_s +   \frac{1}{\beta_0} \, \sqrt{(1+p_z)^2 + \wt m^2}
+  (1 + p_z) - a_s +   \frac{1}{\beta_0} \, \sqrt{(1+p_z)^2 + \wt m^2}
   \label{hpapa}
 \end{equation}
 

bmad/doc/cover-page.tex

@@ -3,7 +3,7 @@
 
 \begin{flushright}
 \large
-  Revision: September 24, 2024 \\
+  Revision: October 10, 2024 \\
 \end{flushright}
 
 \pdfbookmark[0]{Preamble}{Preamble} 

bmad/low_level/make_mat6_tracking.f90

@@ -42,6 +42,7 @@ subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
 err_flag = .true.
 del_orb = bmad_com%d_orb
 abs_p = max(abs(start_orb%vec(2)) + abs(del_orb(2)), abs(start_orb%vec(4)) + abs(del_orb(4)), abs(del_orb(6)))
+bmad_private%random_on = .false.
 
 ! The factor of 1.01 is used to avoid roundoff problems.
 ! Note: init_coord is avoided since init_coord will make z and t consistent with the element's t_ref.
@@ -53,6 +54,7 @@ subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
 call track1 (start_orb0, ele, param, end_orb)
 if (end_orb%state /= alive$) then
   call out_io (s_error$, r_name, 'PARTICLE LOST IN TRACKING CENTRAL PARTICLE. MATRIX NOT CALCULATED FOR ELEMENT: ' // ele%name)
+  bmad_private%random_on = .true.
   return
 endif
 
@@ -66,6 +68,7 @@ subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
   call track1 (start, ele, param, end2)
   if (end2%state /= alive$) then
     call out_io (s_error$, r_name, 'PARTICLE LOST IN TRACKING (+). MATRIX NOT CALCULATED FOR ELEMENT: ' // ele%name)
+    bmad_private%random_on = .true.
     return
   endif
 
@@ -76,6 +79,7 @@ subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
   call track1 (start, ele, param, end1)
   if (end1%state /= alive$) then
     call out_io (s_error$, r_name, 'PARTICLE LOST IN TRACKING (-). MATRIX NOT CALCULATED FOR ELEMENT: ' // ele%name)
+    bmad_private%random_on = .true.
     return
   endif
 
@@ -86,6 +90,8 @@ subroutine make_mat6_tracking (ele, param, start_orb, end_orb, err_flag)
 
 ele%vec0 = end_orb%vec - matmul(mat6, start_orb%vec)
 ele%mat6 = mat6
+
+bmad_private%random_on = .true.
 err_flag = .false.
 
 !------------------------------------------------------

bmad/low_level/track_a_foil.f90

@@ -19,7 +19,7 @@
 !   mat6(6,6)  -- real(rp), optional: Transfer matrix through the element.
 !-
 
-subroutine track_a_foil (orbit, ele, param, mat6, make_matrix)
+recursive subroutine track_a_foil (orbit, ele, param, mat6, make_matrix)
 
 use bmad_interface, except_dummy => track_a_foil
 use random_mod
@@ -84,7 +84,9 @@ subroutine track_a_foil (orbit, ele, param, mat6, make_matrix)
       end select
     enddo
 
-    if (is_true(ele%value(scatter_test$))) then
+    if (.not. bmad_private%random_on) then
+      rnd = 0
+    elseif (is_true(ele%value(scatter_test$))) then
       rnd = 1
     else
       call ran_gauss(rnd)

bmad/modules/bmad_routine_interface.f90

@@ -2727,7 +2727,7 @@ subroutine track_a_sol_quad (orbit, ele, param, mat6, make_matrix)
   logical, optional :: make_matrix
 end subroutine
 
-subroutine track_a_foil (orbit, ele, param, mat6, make_matrix)
+recursive subroutine track_a_foil (orbit, ele, param, mat6, make_matrix)
   import
   implicit none
   type (coord_struct) orbit

bmad/modules/bmad_struct.f90

@@ -2263,10 +2263,11 @@ module bmad_struct
 type (bmad_common_struct), save, target :: bmad_com
 
 ! Bmad global private structure
-! For communication between Bmad and Bmad based programs.
+! For communication between Bmad routines and Bmad based programs.
 
 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 :: random_on = .true.       ! Temporarily turned off, for example, with the closed orbit calc.
 end type
 
 type (bmad_private_struct), save, target :: bmad_private

tao/code/tao_data_and_eval_mod.f90

@@ -1688,11 +1688,14 @@ subroutine tao_evaluate_stack (stack, n_size_in, use_good_user, value, err_flag,
     i2 = i2 + 1
     call value_transfer (stk2(i2)%value, stack(i)%value)
 
-  case (species_const$) 
+  case (species_const$) ! Something like "electron". Just push on stack.
     i2 = i2 + 1
     stk2(i2)%name = stack(i)%name
     call re_allocate(stk2(i2)%value, 1)
 
+  case (species$)
+    stk2(i2)%value = species_id(stk2(i2)%name)
+
   case (lat_num$, ele_num$)
     !!! This needs to be fixed to include default stuff
     !!! call tao_param_value_routine (stack(i)%name, '', stack(i), err_flag, print_err)

tao/code/tao_evaluate_expression.f90

@@ -83,8 +83,8 @@ subroutine tao_evaluate_expression (expression, n_size, use_good_user, value, er
 character(40) saved_prefix
 character(*), parameter :: r_name = "tao_evaluate_expression"
 
-logical delim_found, do_combine, use_good_user
-logical err_flag, err, wild, printit, found, species_here
+logical delim_found, do_combine, use_good_user, in_species_func
+logical err_flag, err, wild, printit, found
 logical, optional :: print_err
 
 ! Don't destroy the input expression
@@ -124,6 +124,7 @@ subroutine tao_evaluate_expression (expression, n_size, use_good_user, value, er
 n_func = 0
 i_lev = 0
 i_op = 0
+in_species_func = .false.
 
 allocate (stk(20), op(20))
 
@@ -133,7 +134,11 @@ subroutine tao_evaluate_expression (expression, n_size, use_good_user, value, er
 
   ! get a word
 
-  call word_read (phrase, '+-*/()^,}[ ', word, ix_word, delim, delim_found, phrase)
+  if (in_species_func) then
+    call word_read (phrase, ')', word, ix_word, delim, delim_found, phrase)
+  else
+    call word_read (phrase, '+-*/()^,}[ ', word, ix_word, delim, delim_found, phrase)
+  endif
 
   ! Args are counted counted at the beginning of the function and at each comma.
 
@@ -218,8 +223,6 @@ subroutine tao_evaluate_expression (expression, n_size, use_good_user, value, er
   do
     if (delim /= '*') exit
 
-
-
     ix0 = index(word, '::')
     ix4 = index(word, '|')
 
@@ -361,6 +364,12 @@ subroutine tao_evaluate_expression (expression, n_size, use_good_user, value, er
 
       call push_op_stack (op, i_op, l_func_parens$)
 
+      ! Parse function argument for functions that take a species.
+      select case (word2)
+      case ('mass_of', 'charge_of', 'species', 'antiparticle', 'anomalous_moment_of')
+        in_species_func = .true.
+      end select
+
     else
       call push_op_stack (op, i_op, l_parens$)
     endif
@@ -390,14 +399,7 @@ subroutine tao_evaluate_expression (expression, n_size, use_good_user, value, er
       endif
 
     else
-      species_here = .false.
-      if (i_op > 1) then
-        select case(op(i_op-1))   ! op(i_op) will be l_func_parens$
-        case (mass_of$, charge_of$, anomalous_moment_of$, antiparticle$, species$);  species_here = .true.
-        end select
-      endif
-
-      if (species_here) then
+      if (in_species_func) then
         call push_stack (stk, i_lev, species_const$)
         stk(i_lev)%name = word
       else

tao/doc/syntax.tex

@@ -153,15 +153,19 @@ \section{Arithmetic Expressions}
   \vn{ran_gauss}()            & Gaussian distributed random number with unit RMS \\
   \vn{ran_gauss}(sig_cut)     & Gaussian distributed random number truncated at sig_cut. \\
   \vn{int(x)}                 & Nearest integer with magnitude less then x \\
-  \vn{nint(x)}                & Nearest integer to x           \\
-  \vn{floor(x)}               & Nearest integer less than x    \\
-  \vn{ceiling(x)}             & Nearest integer greater than x \\
+  \vn{nint(x)}                & Nearest integer to x             \\
+  \vn{floor(x)}               & Nearest integer less than x      \\
+  \vn{ceiling(x)}             & Nearest integer greater than x   \\
   \vn{modulo(a, p)}           & a - floor(a/p) * p. Will be in range [0, p]. \\
-  \vn{average(arr)}           & Average value of an array      \\
-  \vn{rms(arr)}               & RMS value of an array          \\
-  \vn{sum(arr)}               & Sum of array values.           \\
-  \vn{min(arr)}               & Minimum of array values.       \\
-  \vn{max(arr)}               & Maximum of array values.
+  \vn{average(arr)}           & Average value of an array        \\
+  \vn{rms(arr)}               & RMS value of an array            \\
+  \vn{sum(arr)}               & Sum of array values.             \\
+  \vn{min(arr)}               & Minimum of array values.         \\
+  \vn{max(arr)}               & Maximum of array values.         \\
+  \vn{mass_of}(A)               & Mass of particle A             \\
+  \vn{charge_of}(A)             & Charge, in units of the elementary charge, of particle A \\
+  \vn{anomalous_moment_of}(A)   & Anomalous magnetic moment of particle A        \\
+  \vn{species}(A)               & Integer ID associated with species A
 \end{tabular} \newline
 Both \vn{ran} and \vn{ran_gauss} use a seeded random number generator.  Setting the seed is
 described in Section~\sref{s:globals}.

tao/version/tao_version_mod.f90

@@ -6,5 +6,5 @@
 !-
 
 module tao_version_mod
-character(*), parameter :: tao_version_date = "2024/10/08 16:23:17"
+character(*), parameter :: tao_version_date = "2024/10/09 03:31:12"
 end module