Merge pull request #679 from bmad-sim/devel/68

Correct bookkeeping for atomic and molecular negatively charged particles.

bmad/code/attribute_bookkeeper.f90

@@ -25,6 +25,7 @@ subroutine attribute_bookkeeper (ele, force_bookkeeping)
 use s_fitting, only: check_bend
 use bookkeeper_mod, except_dummy => attribute_bookkeeper
 use xraylib_interface, except_dummy2 => attribute_bookkeeper
+use xraylib, dummy => r_e
 use super_recipes_mod, only: super_brent
 use ptc_layout_mod, only: update_ele_from_fibre
 use taylor_mod, only: kill_taylor
@@ -47,7 +48,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
+integer i, j, ix, ig, n, n_div, ixm, ix_pole_max, particle, geometry, i_max, status, material, z_material
 
 character(20) :: r_name = 'attribute_bookkeeper'
 
@@ -444,6 +445,31 @@ subroutine attribute_bookkeeper (ele, force_bookkeeping)
  val(rf_wavelength$) = 0
  endif
 
+! Foil
+
+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$)
+ 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
+
+ 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
+
 ! Crystal
 
 case (crystal$, multilayer_mirror$, mirror$, detector$, sample$, diffraction_plate$)

bmad/code/lat_sanity_check.f90

@@ -16,6 +16,7 @@
 subroutine lat_sanity_check (lat, err_flag)
 
 use bmad_interface, except_dummy => lat_sanity_check
+use xraylib, dummy => r_e
 
 implicit none
 
@@ -590,7 +591,25 @@ subroutine lat_sanity_check (lat, err_flag)
 
  endif
 
- ! Zero length cavity is a verboten
+ ! Foil
+
+ if (ele%key == foil$) then
+ if (atomic_number(ele%ref_species)== 0) then
+ call out_io (s_fatal$, r_name, &
+ 'ELEMENT: ' // ele_full_name(ele, '@N (&#)') // 'HAS REFERENCE SPECIES: ' // species_name(ele%ref_species), &
+ '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
 
  if (ele%key == lcavity$ .and. ele%value(l$) == 0) then
  call out_io (s_fatal$, r_name, &

bmad/code/type_ele.f90

@@ -1436,7 +1436,7 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c
 character(40) a_name, a2_name
 logical is_2nd_col_attrib
 
-character(41), parameter :: att_name(85) = [character(40):: 'X_PITCH', 'Y_PITCH', 'X_OFFSET', &
+character(41), parameter :: att_name(88) = [character(40):: 'X_PITCH', 'Y_PITCH', 'X_OFFSET', &
  'Y_OFFSET', 'Z_OFFSET', 'REF_TILT', 'TILT', 'ROLL', 'X1_LIMIT', 'Y1_LIMIT', &
  'FB1', 'FQ1', 'LORD_PAD1', 'HKICK', 'VKICK', 'KICK', 'FRINGE_TYPE', 'DS_STEP', 'R0_MAG', &
  'KS', 'K1', 'K2', 'G', 'DG', 'G_TOT', 'H1', 'E1', 'FINT', 'HGAP', &
@@ -1449,9 +1449,9 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c
  'ETA_Y0', 'ETAP_Y0', 'KICK0', 'X0', 'PX0', 'Y0', 'PY0', 'Z0', 'PZ0', &
  'C11_MAT0', 'C12_MAT0', 'C21_MAT0', 'C22_MAT0', 'HARMON', &
  'MODE_FLIP0', 'BETA_A_STRONG', 'BETA_B_STRONG', 'REF_TIME_START', &
- 'PX_KICK', 'PY_KICK', 'PZ_KICK']
+ 'PX_KICK', 'PY_KICK', 'PZ_KICK', 'DENSITY', 'RADIATION_LENGTH', 'AREA_DENSITY']
 
-character(41), parameter :: att2_name(85) = [character(40):: 'X_PITCH_TOT', 'Y_PITCH_TOT', 'X_OFFSET_TOT', &
+character(41), parameter :: att2_name(88) = [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', &
  'FB2', 'FQ2', 'LORD_PAD2', 'BL_HKICK', 'BL_VKICK', 'BL_KICK', 'FRINGE_AT', 'NUM_STEPS', 'R0_ELEC', &
  'BS_FIELD', 'B1_GRADIENT', 'B2_GRADIENT', 'B_FIELD', 'DB_FIELD', 'B_FIELD_TOT', 'H2', 'E2', 'FINTX', 'HGAPX', &
@@ -1464,7 +1464,7 @@ function is_2nd_column_attribute (ele, attrib_name, ix2_attrib) result (is_2nd_c
  'ETA_Y1', 'ETAP_Y1', 'MATRIX', 'X1', 'PX1', 'Y1', 'PY1', 'Z1', 'PZ1', &
  'C11_MAT1', 'C12_MAT1', 'C21_MAT1', 'C22_MAT1', 'HARMON_MASTER', &
  'MODE_FLIP1', 'ALPHA_A_STRONG', 'ALPHA_B_STRONG', 'DELTA_REF_TIME', &
- 'X_KICK', 'Y_KICK', 'Z_KICK']
+ 'X_KICK', 'Y_KICK', 'Z_KICK', 'DENSITY_USED', 'RADIATION_LENGTH_USED', 'AREA_DENSITY_USED']
 
 ! Exceptional cases
 

bmad/doc/beam-init.tex

@@ -73,7 +73,7 @@ \section{Beam_Init_Struct Structure}
 %
 \item[\%a_emit, \%b_emit, \%a_norm_emit, \%b_norm_emit] \Newline
 Normalized and unnormalized emittances. Either \vn{a_norm_emit} or \vn{a_emit} may be set but not
-both. similarly, either \vn{b_norm_emit} or \vn{b_emit} may be set but not both. 
+both. similarly, either \vn{b_norm_emit} or \vn{b_emit} may be set but not both.
 
 When simulating a ring, if any of these parameters is set negative, and if the \bmad based program
 being run has enabled it, the value of the parameter will be set the value as calculated from the

bmad/doc/cover-page.tex

@@ -3,7 +3,7 @@
 
 \begin{flushright}
 \large
- Revision: November 20, 2023 \\
+ Revision: December 6, 2023 \\
 \end{flushright}
 
 \pdfbookmark[0]{Preamble}{Preamble} 

bmad/doc/elements.tex

@@ -1133,7 +1133,7 @@ \section{Crab_Cavity}
  phi0_multipass = <Real> ! Phase (rad/2\(\pi\)) with respect to a multipass lord (\sref{s:multipass}).
  rf_frequency = <Real> ! RF frequency (Hz).
  harmon = <Real> ! Harmonic number
- harmon_master = <Logic> ! Is haromin or rf_frequency the dependent var with ref energy changes?
+ harmon_master = <Logic> ! Is harmon or rf_frequency the dependent var with ref energy changes?
  voltage ! Cavity voltage. Dependent attribute (\sref{s:depend}).
 \end{example}
 
@@ -2080,8 +2080,9 @@ \section{Foil}
 \label{s:foil}
 \index{foil|hyperbf}
 
-A \vn{foil} element represents a planar sheet which strips all the electrons from a particle.
-Additionally, there will be scattering of the trajectory of a particle going through a \vn{foil}.
+A \vn{foil} element represents a planar sheet of material which strips electrons from a particle.
+In conjunction, there will be scattering of the particle trajectory as well as an associated energy
+loss.
 
 General \vn{foil} attributes are:
 \begin{center}
@@ -2101,14 +2102,35 @@ \section{Foil}
 
 Attributes specific to a \vn{foil} element are:
 \begin{example}
- material_type = <String> ! Foil material.
- thickness = <Real> ! Foil thickness.
- final_charge = <Integer> ! Final charge state
-\end{example}
-
-The \vn{thickness} parameter along with \vn{material_type} is used to calculate the scattering
-sigma. Scattering is simulated to be Gaussian distributed with a sigma given by Eq.~(6) of Lynch and
-Dahl\cite{b:lynch}.
+ material_type = <String> ! Foil material.
+ thickness = <Real> ! Material thickness (m).
+ density = <Real> ! Input material density (kg/m^3).
+ density_used ! Density value used in tracking (kg/m^3). Dependent parameter.
+ radiation_length = <Real> ! Input material radiation length (m).
+ radiation_length_used ! Radiation length used in tracking (m). Dependent parameter.
+ area_density = <Real> ! Input material area density (kg/m^2).
+ area_density_used ! Area density used in tracking (kg/m^2).
+ final_charge = <Integer> ! Final charge state
+\end{example}
+
+Scattering is simulated to be Gaussian distributed with a sigma given by Eq.~(6) of Lynch and
+Dahl\cite{b:lynch}. Energy loss is calculated using the Bethe-Bloch formula.
+
+The radiation length used in the scattering calculation is given by the \vn{radiation_length_used}
+parameter. This value cannot be set directly. Rather, if the \vn{radiation_length} parameter is set
+non-zero, this value will be transferred to \vn{radiation_length_used}. If \vn{radiation_length} is
+zero (the default), the value of \vn{radiation_length_used} will be set by \bmad using the measured
+value from the published literature.
+
+Similarly, the \vn{area_density_used} (density of the material per unit of surface area) value
+needed for the calculation is not set directly but is set in one of two ways depending upon if the
+material \vn{thickness} is non-zero or not. If \vn{thickness} is non-zero, \vn{area_density_used} is
+set by the product of \vn{thickness} and \vn{density_used} while the value of \vn{density_used} is
+set by \bmad to be either the value \vn{density} if the \vn{density} is non-zero or by the measured
+density of the material as given in the published literature. If \vn{thickness} is zero (the
+default), the value of \vn{area_density_used} is set equal to the value of \vn{area_density}. Note that
+the value of \vn{density_used} is only used to set \vn{area_density_used} when \vn{thickness} is non-zero
+and otherwise does not affect the tracking calculation.
 
 In terms of element placement, The length of a \vn{foil} element (\Eq{l00l}) is considered to be
 zero. This is similar to the \vn{beambeam} element which is also considered to have zero length but
@@ -2119,11 +2141,10 @@ \section{Foil}
 
 Particles going through the \vn{foil} are stripped to have a final charge given by \vn{final_charge}.
 
-Energy loss is calculated using the Bethe-Bloch formula.
 
 Example:
 \begin{example}
- septum: foil, material_type = "Cu", thickness = 0.127
+ septum: foil, material_type = "Cu", thickness = 0.127, radiation_length = 
 \end{example}
 
 \newpage