From 0a0d1e496c4aff494dc2a108e09cb46c4d96ce8f Mon Sep 17 00:00:00 2001 From: David Sagan Date: Thu, 10 Oct 2024 03:18:47 +0200 Subject: [PATCH] Fix Tao handling of mass/charge in expressions. (#1227) --- bmad/code/closed_orbit_calc.f90 | 2 ++ bmad/code/make_mat6.f90 | 2 ++ bmad/doc/charged-tracking.tex | 8 +++---- bmad/doc/cover-page.tex | 2 +- bmad/low_level/make_mat6_tracking.f90 | 6 ++++++ bmad/low_level/track_a_foil.f90 | 6 ++++-- bmad/modules/bmad_routine_interface.f90 | 2 +- bmad/modules/bmad_struct.f90 | 3 ++- tao/code/tao_data_and_eval_mod.f90 | 5 ++++- tao/code/tao_evaluate_expression.f90 | 28 +++++++++++++------------ tao/doc/syntax.tex | 20 +++++++++++------- tao/version/tao_version_mod.f90 | 2 +- 12 files changed, 54 insertions(+), 32 deletions(-) diff --git a/bmad/code/closed_orbit_calc.f90 b/bmad/code/closed_orbit_calc.f90 index 47bc00d56c..03030d0fb1 100644 --- a/bmad/code/closed_orbit_calc.f90 +++ b/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) diff --git a/bmad/code/make_mat6.f90 b/bmad/code/make_mat6.f90 index 07cff674ea..4f1c84f32d 100644 --- a/bmad/code/make_mat6.f90 +++ b/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 diff --git a/bmad/doc/charged-tracking.tex b/bmad/doc/charged-tracking.tex index dcadbf7f15..097b7bed9e 100644 --- a/bmad/doc/charged-tracking.tex +++ b/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} diff --git a/bmad/doc/cover-page.tex b/bmad/doc/cover-page.tex index a9cf2a56f3..2f117f15d6 100644 --- a/bmad/doc/cover-page.tex +++ b/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} diff --git a/bmad/low_level/make_mat6_tracking.f90 b/bmad/low_level/make_mat6_tracking.f90 index 71bb9e73de..b2db0bd098 100644 --- a/bmad/low_level/make_mat6_tracking.f90 +++ b/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. !------------------------------------------------------ diff --git a/bmad/low_level/track_a_foil.f90 b/bmad/low_level/track_a_foil.f90 index 2ab60c0855..ea36b63157 100644 --- a/bmad/low_level/track_a_foil.f90 +++ b/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) diff --git a/bmad/modules/bmad_routine_interface.f90 b/bmad/modules/bmad_routine_interface.f90 index 4b4ecdcf8d..2e711bb3fb 100644 --- a/bmad/modules/bmad_routine_interface.f90 +++ b/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 diff --git a/bmad/modules/bmad_struct.f90 b/bmad/modules/bmad_struct.f90 index 78848f3034..7b6d895b7d 100644 --- a/bmad/modules/bmad_struct.f90 +++ b/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 diff --git a/tao/code/tao_data_and_eval_mod.f90 b/tao/code/tao_data_and_eval_mod.f90 index 06fe0c686b..06490d2ce1 100644 --- a/tao/code/tao_data_and_eval_mod.f90 +++ b/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) diff --git a/tao/code/tao_evaluate_expression.f90 b/tao/code/tao_evaluate_expression.f90 index 613f2ed930..05c8d1ecfd 100644 --- a/tao/code/tao_evaluate_expression.f90 +++ b/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 diff --git a/tao/doc/syntax.tex b/tao/doc/syntax.tex index 30f4daa954..00ecd16d79 100644 --- a/tao/doc/syntax.tex +++ b/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}. diff --git a/tao/version/tao_version_mod.f90 b/tao/version/tao_version_mod.f90 index e09d4b8655..2cb4528b9a 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/10/08 16:23:17" +character(*), parameter :: tao_version_date = "2024/10/09 03:31:12" end module