Start work on extending foil tracking to Lynch-Dahl scattering algorithm.

.github/workflows/dcs16_push.yml

@@ -13,6 +13,7 @@ env:
 
 # secrets.DCS_PUSH allows the ci.yml workflow to be triggered by this workflow. Documentation at:
 #		https://github.com/peter-evans/create-pull-request/blob/main/docs/concepts-guidelines.md#triggering-further-workflow-runs
+# DCS_PUSH is defined: in bmad-ecosystem -> Settings -> Secrets and variables -> Actions
 
 jobs:
   create-pull-request:

bmad/code/attribute_bookkeeper.f90

@@ -29,6 +29,7 @@ subroutine attribute_bookkeeper (ele, force_bookkeeping)
 use super_recipes_mod, only: super_brent
 use ptc_layout_mod, only: update_ele_from_fibre
 use taylor_mod, only: kill_taylor
+use particle_species_mod
 
 implicit none
 
@@ -40,6 +41,8 @@ subroutine attribute_bookkeeper (ele, force_bookkeeping)
 type (photon_element_struct), pointer :: ph
 type (wake_lr_mode_struct), pointer :: lr
 type (converter_prob_pc_r_struct), pointer :: ppcr
+type (molecular_component_struct), allocatable :: component(:)
+type (material_struct), pointer :: material
 
 real(rp) factor, e_factor, gc, f2, phase, E_tot, polarity, dval(num_ele_attrib$), time, beta
 real(rp) w_inv(3,3), len_old, f, dl, b_max, zmin, ky, kz
@@ -48,7 +51,7 @@ subroutine attribute_bookkeeper (ele, force_bookkeeping)
 real(rp) kick_magnitude, bend_factor, quad_factor, radius0, step_info(7), dz_dl_max_err
 real(rp) a_pole(0:n_pole_maxx), b_pole(0:n_pole_maxx)
 
-integer i, j, ix, ig, n, n_div, ixm, ix_pole_max, particle, geometry, i_max, status, material, z_material
+integer i, j, ix, ig, n, n_div, ixm, ix_pole_max, particle, geometry, i_max, status, z_material
 
 character(20) ::  r_name = 'attribute_bookkeeper'
 
@@ -449,26 +452,39 @@ subroutine attribute_bookkeeper (ele, force_bookkeeping)
 
 case (foil$)
 
-  material = species_id(ele%component_name)
-
-  if (ele%value(radiation_length$) == 0) then
-    ele%value(radiation_length_used$) = x0_radiation_length(material)
-  else
-    ele%value(radiation_length_used$) = ele%value(radiation_length$)
+  call molecular_components(ele%component_name, component)
+  n = size(component)
+  if (.not. allocated(ele%foil%material)) allocate(ele%foil%material(n))
+  if (n /= size(ele%foil%material)) then
+    call out_io(s_error$, r_name, 'NUMBER OF COMPONENTS IN: ' // quote(ele%component_name) // ' (' // int_str(n) // &
+                    ') IS NOT THE SAME AS OTHER PARAMETERS.')
+    if (global_com%exit_on_error) call err_exit
+    return
   endif
 
-  if (ele%value(density$) == 0) then
-    z_material = atomic_number(material)
-    ele%value(density_used$) = ElementDensity(z_material) * 1e3_rp  ! From xraylib. Convert to kg/m^3
-  else
-    ele%value(density_used$) = ele%value(density$)
-  endif
+  do ix = 1, n
+    material => ele%foil%material(ix)
+    material%species = species_id(component(ix)%atom)
+    z_material = atomic_number(material%species)
 
-  if (ele%value(thickness$) == 0) then
-    ele%value(area_density_used$) = ele%value(area_density$)
-  else
-    ele%value(area_density_used$) = ele%value(density_used$) * ele%value(thickness$) 
-  endif
+    if (material%radiation_length == 0) then
+      material%radiation_length_used = x0_radiation_length(material%species)
+    else
+      material%radiation_length_used = material%radiation_length
+    endif
+
+    if (material%density == 0) then
+      material%density_used = ElementDensity(z_material) * 1e3_rp  ! From xraylib. Convert to kg/m^3
+    else
+      material%density_used = material%density
+    endif
+
+    if (ele%value(thickness$) == 0) then
+      material%area_density_used = material%area_density
+    else
+      material%area_density_used = material%density_used * ele%value(thickness$)
+    endif
+  enddo
 
 ! Crystal
 

bmad/code/lat_sanity_check.f90

@@ -600,13 +600,6 @@ subroutine lat_sanity_check (lat, err_flag)
                       'WHICH DOES NOT HAVE AN ASSOCIATED ATOMIC NUMBER.')
         err_flag = .true.
       endif
-
-      if (ele%value(radiation_length_used$) == real_garbage$) then
-        call out_io(s_fatal$, r_name, &
-            'ELEMENT: ' // ele_full_name(ele, '@N (&#)'), & 
-            'CANNOT HANDLE NON-ATOMIC MATERIAL_TYPE: ' // ele%component_name)
-        err_flag = .true.
-      endif
     endif
 
     ! Zero length cavity is verboten

bmad/code/pointer_to_attribute.f90

@@ -722,7 +722,6 @@ subroutine pointer_to_attribute (ele, attrib_name, do_allocation, a_ptr, err_fla
 select case (a_name)
 ! attrib_type = is_real$
 ! attrib_type = is_logical$
-case ('SCATTER');                        a_ptr%r => ele%value(scatter$)
 case ('MATRIX');                         a_ptr%r => ele%value(matrix$)
 case ('KICK0');                          a_ptr%r => ele%value(kick0$)
 case ('FLEXIBLE');                       a_ptr%r => ele%value(flexible$)
@@ -771,6 +770,7 @@ subroutine pointer_to_attribute (ele, attrib_name, do_allocation, a_ptr, err_fla
 case ('PTC_FIELD_GEOMETRY');             a_ptr%r => ele%value(ptc_field_geometry$)
 case ('REF_ORBIT_FOLLOWS');              a_ptr%r => ele%value(ref_orbit_follows$)
 case ('REF_COORDS');                     a_ptr%r => ele%value(ref_coords$)
+case ('SCATTER_METHOD');                 a_ptr%r => ele%value(scatter_method$)
 case ('SPACE_CHARGE_METHOD');            a_ptr%i => ele%space_charge_method
 case ('SPIN_TRACKING_METHOD');           a_ptr%i => ele%spin_tracking_method
 case ('TRACKING_METHOD');                a_ptr%i => ele%tracking_method

bmad/code/type_ele.f90

@@ -80,6 +80,7 @@ subroutine type_ele (ele, type_zero_attrib, type_mat6, type_taylor, twiss_out, t
 type (str_index_struct) str_index
 type (rad_map_struct), pointer :: rm0, rm1
 type (photon_reflect_table_struct), pointer :: rt
+type (material_struct), pointer :: matter
 
 integer, optional, intent(in) :: type_mat6, twiss_out, type_field
 integer, optional, intent(out) :: n_lines
@@ -401,14 +402,35 @@ subroutine type_ele (ele, type_zero_attrib, type_mat6, type_taylor, twiss_out, t
   call encode_2nd_column_parameter (li, nl2, nl, 'REF_SPECIES', str_val = species_name(ele%ref_species))
 endif
 
+! Foil
+
+if (associated(ele%foil)) then
+    nl=nl+1; li(nl) = ''
+    nl=nl+1; li(nl) = 'Material_type: ' // ele%component_name
+
+    do ix = 1, size(ele%foil%material)
+      matter = ele%foil%material(ix)
+      if (size(ele%foil%material) > 1) then
+        nl=nl+1; li(nl) = ''
+        nl=nl+1; li(nl) = 'Component: ' // species_name(matter%species)
+      endif
+
+      nl=nl+1; write(li(nl), '(3(a, es14.6))') 'Density      =', matter%density, &
+                  'Area_Density      =', matter%area_density,      'Radiation_Length      =', matter%radiation_length
+      nl=nl+1; write(li(nl), '(3(a, es14.6))') 'Density_Used =', matter%density_used, &
+                  'Area_Density_Used =', matter%area_density_used, 'Radiation_Length_Used =', matter%radiation_length_used
+      nl=nl+1; 
+    enddo
+endif
+
 ! Converter
 
 if (associated(ele%converter)) then
   do im = 1, size(ele%converter%dist)
   enddo
 endif
 
-! Cartesian map
+! Cartesian map. The type_field logical is useful since field tables can be very large.
 
 if (associated(ele%cartesian_map)) then
   if (integer_option(no$, type_field) == no$) then
@@ -1453,10 +1475,11 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c
                 'PZ_APERTURE_WIDTH2', 'Z_APERTURE_WIDTH2', 'CMAT_11', 'CMAT_21', 'X_DISPERSION_ERR', &
                 'X_DISPERSION_CALIB', 'K1X', 'RF_FREQUENCY', 'UPSTREAM_ELE_DIR', 'SIG_X', &
                 'BETA_A0', 'BETA_B0', 'ALPHA_A0', 'ALPHA_B0', 'ETA_X0', 'ETAP_X0', 'X1_EDGE', 'Y1_EDGE', &
-                'ETA_Y0', 'ETAP_Y0', 'KICK0', 'X0', 'PX0', 'Y0', 'PY0', 'Z0', 'PZ0', &
+                'ETA_Y0', 'ETAP_Y0', 'KICK0', 'X0', 'PX0', 'Y0', 'PY0', 'Z0', 'PZ0', 'ATOMIC_WEIGHT', &
                 'C11_MAT0', 'C12_MAT0', 'C21_MAT0', 'C22_MAT0', 'HARMON', 'FINAL_CHARGE', &
-                'MODE_FLIP0', 'BETA_A_STRONG', 'BETA_B_STRONG', 'REF_TIME_START', &
-                'PX_KICK', 'PY_KICK', 'PZ_KICK', 'DENSITY', 'RADIATION_LENGTH', 'AREA_DENSITY']
+                'MODE_FLIP0', 'BETA_A_STRONG', 'BETA_B_STRONG', 'REF_TIME_START', 'THICKNESS', &
+                'PX_KICK', 'PY_KICK', 'PZ_KICK', &
+                'F_FACTOR']
 
 character(41), parameter :: att2_name(91) = [character(40):: 'X_PITCH_TOT', 'Y_PITCH_TOT', 'X_OFFSET_TOT', &
                 'Y_OFFSET_TOT', 'Z_OFFSET_TOT', 'REF_TILT_TOT', 'TILT_TOT', 'ROLL_TOT', 'X2_LIMIT', 'Y2_LIMIT', &
@@ -1468,10 +1491,11 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c
                 'PZ_APERTURE_CENTER', 'Z_APERTURE_CENTER', 'CMAT_12', 'CMAT_22', 'Y_DISPERSION_ERR', &
                 'Y_DISPERSION_CALIB', 'K1Y', 'RF_WAVELENGTH', 'DOWNSTREAM_ELE_DIR', 'SIG_Y', &
                 'BETA_A1', 'BETA_B1', 'ALPHA_A1', 'ALPHA_B1', 'ETA_X1', 'ETAP_X1', 'X2_EDGE', 'Y2_EDGE', &
-                'ETA_Y1', 'ETAP_Y1', 'MATRIX', 'X1', 'PX1', 'Y1', 'PY1', 'Z1', 'PZ1', &
+                'ETA_Y1', 'ETAP_Y1', 'MATRIX', 'X1', 'PX1', 'Y1', 'PY1', 'Z1', 'PZ1', 'Z_CHARGE', &
                 'C11_MAT1', 'C12_MAT1', 'C21_MAT1', 'C22_MAT1', 'HARMON_MASTER', 'SCATTER', &
-                'MODE_FLIP1', 'ALPHA_A_STRONG', 'ALPHA_B_STRONG', 'DELTA_REF_TIME', &
-                'X_KICK', 'Y_KICK', 'Z_KICK', 'DENSITY_USED', 'RADIATION_LENGTH_USED', 'AREA_DENSITY_USED']
+                'MODE_FLIP1', 'ALPHA_A_STRONG', 'ALPHA_B_STRONG', 'DELTA_REF_TIME', 'DREL_THICKNESS_DX', &
+                'X_KICK', 'Y_KICK', 'Z_KICK', &
+                'SCATTER_METHOD']
 
 ! Exceptional cases
 
@@ -1496,9 +1520,11 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c
 end select
 
 ! Is a 2nd column attribute if corresponding first column attribute exists
+! Note: There are rare cases where the corresponding 1st column parameter does not exist for the particular element type.
 
 call match_word (attrib_name, att2_name, ix, .true., .false.)
 if (ix > 0) then
+  if (.not. has_attribute(ele, att_name(ix))) return
   ia = attribute_index(ele, att_name(ix))
   is_2nd_col_attrib = (ia > 0)
   return
@@ -1507,7 +1533,10 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c
 ! If the attribute has a corresponding 2nd column attribute, set ix2_attrib accordingly.
 
 call match_word (attrib_name, att_name, ix, .true., .false.)
-if (ix > 0) ix2_attrib = attribute_index(ele, att2_name(ix))
+if (ix > 0) then
+  if (.not. has_attribute(ele, att2_name(ix))) return
+  ix2_attrib = attribute_index(ele, att2_name(ix))
+endif
 
 end function is_2nd_column_attribute
 

bmad/code/write_bmad_lattice_file.f90

@@ -79,6 +79,7 @@ subroutine write_bmad_lattice_file (bmad_file, lat, err, output_form, orbit0)
 type (expression_atom_struct), pointer :: stack(:)
 type (str_index_struct) str_index
 type (lat_ele_order_struct) order
+type (material_struct), pointer :: material
 
 real(rp) s0, x_lim, y_lim, val, x, y
 
@@ -534,6 +535,48 @@ subroutine write_bmad_lattice_file (bmad_file, lat, err, output_form, orbit0)
       line = trim(line) // '}'
     endif
 
+    ! Foil
+
+    if (associated(ele%foil)) then
+      if (size(ele%foil%material) > 1) then
+        if (any(ele%foil%material%density /= 0)) then
+          line = trim(line) // ', density = ('
+          do n = 1, size(ele%foil%material)
+            if (n == 1) then; line = trim(line) // re_str(ele%foil%material(n)%density)
+            else;             line = trim(line) // ', ' // re_str(ele%foil%material(n)%density)
+            endif
+          enddo
+          line = trim(line) // ')'
+        endif
+
+        if (any(ele%foil%material%area_density /= 0)) then
+          line = trim(line) // ', area_density = ('
+          do n = 1, size(ele%foil%material)
+            if (n == 1) then; line = trim(line) // re_str(ele%foil%material(n)%area_density)
+            else;             line = trim(line) // ', ' // re_str(ele%foil%material(n)%area_density)
+            endif
+          enddo
+          line = trim(line) // ')'
+        endif
+
+        if (any(ele%foil%material%radiation_length /= 0)) then
+          line = trim(line) // ', radiation_length = ('
+          do n = 1, size(ele%foil%material)
+            if (n == 1) then; line = trim(line) // re_str(ele%foil%material(n)%radiation_length)
+            else;             line = trim(line) // ', ' // re_str(ele%foil%material(n)%radiation_length)
+            endif
+          enddo
+          line = trim(line) // ')'
+        endif
+
+      else
+        material => ele%foil%material(1)
+        if (material%density /= 0)          line = trim(line) // ', density = ' // re_str(material%density)
+        if (material%area_density /= 0)     line = trim(line) // ', area_density = ' // re_str(material%area_density)
+        if (material%radiation_length /= 0) line = trim(line) // ', radiation_length = ' // re_str(material%radiation_length)
+      endif
+    endif
+
     ! Cartesian_map.
 
     if (associated(ele%cartesian_map)) then

bmad/doc/charged-tracking.tex

@@ -763,16 +763,6 @@ \section{Drift Tracking}
   p_l = \sqrt{1 - \frac{p_x^2 + p_y^2}{(1 + p_z)^2}}
 \end{equation}
 
-%---------------------------------------------------------------------------------
-%---------------------------------------------------------------------------------
-%\section{FFT Grid Tracking with Space Charge}
-%\label{s:egun.sc}
-%\index{e_gun}
-%\index{space charge}
-
-
-
-
 %---------------------------------------------------------------------------------
 %---------------------------------------------------------------------------------
 \section{ElSeparator Tracking}
@@ -849,7 +839,99 @@ \section{ElSeparator Tracking}
 
 %---------------------------------------------------------------------------------
 %---------------------------------------------------------------------------------
-\section{Kicker, Hkicker, and Vkicker, Tracking}
+\section{Foil Tracking}
+\label{s:foil.std}
+\index{foil}
+
+A particle going through a \vn{foil} element is scattered (has an angular change in direction of
+propagation), has an energy loss, and the charge of the particle may be affected. The following two
+subsections give the formulas used for scattering and energy loss.
+
+%---------------------------------------------------------------------------------
+\subsection{Scattering in a Foil}
+\label{s:foil.scatter}
+
+For the scattering part of the simulation, the user can select between one of two algorithms, both
+of which are given in the paper by Lynch and Dahl\cite{b:lynch}. Both methods scatter using the
+equation
+\begin{equation}
+  (dp_x, dp_y) = \frac{p \, \sigma}{P_0} \, (r_1, r_2)
+  \label{dpr1s}
+\end{equation}
+where $dp_x$ and $dp_y$ are the change of the $p_x$ and $p_y$ phase space coordinates, $p$ is the
+particle momentum, $P_0$ is the reference momentum, $r_1$ and $r_2$ are Gaussian random numbers with
+unit sigma and zero mean, and $\sigma$ is the sigma of the angular scattering distribution. The
+factor of $p/P_0$ is due to a translation between change in angle and change in phase space momenta
+(see \Eq{xpa1p}). The actual scattering distribution has $1/\theta^4$ tails ($\theta$ is the
+scattering angle) due to single event large angle scattering (Rutherford scattering). By assuming a
+Gaussian disribution, these tails are not simulated.
+
+The \vn{Highland} algorithm uses Eq~(6) of Lynch and Dahl to calculate the scattering sigma:
+\begin{equation}
+  \sigma = \frac{S_2 \, z}{p \, \beta} \sqrt{\frac{X}{X_0}} \, \left[
+  1 + \epsilon \, \log_{10} \left( \frac{X \, z^2}{X_0 \, \beta^2} \right) \right]
+  \label{sszpb}
+\end{equation}
+where the constant $S_2$ has a value of $13.6 \cdot 10^6$~{eV}, $X$ is the foil thickness, $X_0$ is
+the material radiation length, $z$ is the particle charge, $\beta$ is the particle relativistic
+beta, and $\epsilon$ is a constant with value 0.88. The form of \Eq{sszpb} is slightly different
+than Lynch and Dahl Eq~(6) due to the use by \bmad of SI units (except for energy which uses units
+of eV). In terms of accuracy, Lynch and Dahl write:
+\begin{quote}
+This form (Highland equation) takes into account the $p$ and $z$ dependence quite well at small
+$Z$, but for large $Z$ and small $X$ the $P$-dependence is not taken into account very well. For
+example, for $10^{-3}$ radiation lengths of lead and $p$ = 0.1 (MeV), it overestimates Moliere
+scattering angle by 25\% for singly charged particles.
+\end{quote}
+
+The \vn{Lynch_Dahl} algormithm uses Eq~(7) of Lynch and Dahl to calculate the scattering sigma:
+\begin{equation}
+  \sigma = \frac{\chi_c^2}{1 + F^2} \left[ \frac{1 + \nu}{\nu} \ln (1 + \nu) - 1 \right]
+\end{equation}
+where 
+\begin{equation}
+  \nu = \frac{0.5 \, \chi_c^2}{\chi_\alpha^2 (1 - F)}
+  \label{ncc1f}
+\end{equation}
+with
+\begin{align}
+  \chi_c^2 &= 0.157 \cdot 10^{11} \, \frac{Z (Z+1) X}{A} \left[ \frac{z}{p \, \beta} \right]^2 \\
+  \chi_\alpha &= 2.007 \cdot 10^{-5} \, \frac{Z^{2/3}}{p^2} 
+  \left[1 + 3.34 \cdot 10^6 \,  \left( \frac{Z \, z \, \alpha}{\beta} \right)^2 \right]
+\end{align}
+and $A$ is the atomic weight, $p$ is the particle momentum, $\alpha$ is the fine structure constant,
+and $F$ is a fit parameter.
+Quoting Lynch and Dahl:
+\begin{quote}
+Much of the difficulty in approximating multiple Coulomb scattering in terms of the radiation length
+(Highland formula) is that the number of radiation lengths is a poor measure of the scattering. We
+can get a much better simple expression for the scattering if we do not use the radiation length. An
+expression that does much better than the previous ones is (and gives the Lynch-Dahl formula)...
+
+...This form (Lynch-Dahl), which is not exact, was motivated from a calculation of the RMS angle of the
+screened Rutherford cross section ... The constant 0.5 in (Eq.~\ref{ncc1f} above) was determined
+empirically. For $F$ anywhere in the range of 90\% to 99.5\% this expression represents Moliere
+scattering to better than 2\%
+\end{quote}
+
+%---------------------------------------------------------------------------------
+\subsection{Energy Loss in a Foil}
+\label{s:foil.eloss}
+
+The particle energy loss per unit length $dE/dx$ through a foil is calculated using the \vn{Bethe-Bloch} formula
+\begin{equation}
+  - \left\langle\frac{dE}{dx}\right\rangle = 
+  \frac{4 \pi}{m_e c^2} \cdot \frac{nz^2}{\beta^2} \cdot \left(\frac{e^2}{4\pi\varepsilon_0}\right)^2 \cdot 
+  \left[\ln \left(\frac{2m_e c^2 \beta^2}{I \cdot (1-\beta^2)}\right) - \beta^2\right]
+  \label{ex4pmc}
+\end{equation}
+where $n$ is the material electron density, $I$ is the mean excitation energy, $z$ is the particle
+charge, $c$ is the speed of light, $\epsilon_0$ is the vacuum permittivity, $\beta = {v}/{c}$, is
+the normalized velocity, and $e$ and $m_e$ the electron charge and rest mass respectively.
+
+%---------------------------------------------------------------------------------
+%---------------------------------------------------------------------------------
+\section{Kicker, Hkicker, and Vkicker Tracking}
 \label{s:kicker.std}
 \index{kicker}
 \index{hkicker}