Merge pull request #528 from bmad-sim/devel/step16

Devel/step16

bmad/code/attribute_bookkeeper.f90

@@ -362,6 +362,10 @@ subroutine attribute_bookkeeper (ele, force_bookkeeping)
  if (val(ds_step$) <= 0) then
  if (val(num_steps$) <= 0 .or. abs(val(l$)) == 0) then
  val(ds_step$) = bmad_com%default_ds_step
+ if (ele%key == wiggler$ .or. ele%key == undulator$) call out_io(s_warn$, r_name, &
+ 'Warning! Element: ' // trim(ele%name) // ' which is a ' // key_name(ele%key), &
+ 'is using the bmad_com%default_ds_step since ds_step is not set in the element.', &
+ 'this may be very inaccurate.')
  else
  val(ds_step$) = abs(val(l$)) / val(num_steps$)
  endif

bmad/code/lat_compute_reference_energy.f90

@@ -543,7 +543,8 @@ subroutine ele_compute_ref_energy_and_time (ele0, ele, param, err_flag)
 if (key == converter$) then
  ele%ref_species = ele%converter%species_out
 elseif (key == foil$) then
- n = atomic_number(ele0%ref_species)
+ if (ele0%value(final_charge$) == real_garbage$) ele0%value(final_charge$) = atomic_number(ele0%ref_species)
+ n = nint(ele0%value(final_charge$))
  ele%ref_species = set_species_charge(ele0%ref_species, n)
 else
  ele%ref_species = ele0%ref_species

bmad/doc/cover-page.tex

@@ -3,7 +3,7 @@
 
 \begin{flushright}
 \large
- Revision: September 13, 2023 \\
+ Revision: September 28, 2023 \\
 \end{flushright}
 
 \pdfbookmark[0]{Preamble}{Preamble} 

bmad/doc/elements.tex

@@ -2056,6 +2056,7 @@ \section{Foil}
 \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
@@ -2069,8 +2070,9 @@ \section{Foil}
 The 6x6 transfer matrix of a \vn{foil} element is the unit matrix. That is, scattering does
 not affect the transfer matrix.
 
-Currently, particles going through the \vn{foil} are considered to be fully stripped. Additionally,
-energy losses are ignored.
+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}

bmad/low_level/track_a_foil.f90

@@ -32,10 +32,10 @@ subroutine track_a_foil (orbit, ele, param, mat6, make_matrix)
 type (lat_param_struct) :: param
 
 real(rp), optional :: mat6(6,6)
-real(rp) x0, xx0, sigma, z, rnd(2), density
+real(rp) x0, xx0, sigma, rnd(2), density, I_excite, mass_material, dE_dx, E0, E1, pc, p2, elec_density
 real(rp), parameter :: S2 = 13.6e6_dp, epsilon = 0.088 ! Factor of 1e6 is due to original formula using MeV/c for momentum
 
-integer material, atomic_num
+integer i, material, z_material, z_part, n_step
 
 logical, optional :: make_matrix
 
@@ -51,21 +51,36 @@ subroutine track_a_foil (orbit, ele, param, mat6, make_matrix)
  return
 endif
 
-atomic_num = atomic_number(material)
-density = ElementDensity(atomic_num) * 1e3_rp ! Convert to kg/m^3
+z_material = atomic_number(material)
+density = ElementDensity(z_material) * 1e3_rp ! From xraylib. Convert to kg/m^3
 
-!
+! Angle scatter
 
-z = atomic_number(orbit%species)
+z_part = atomic_number(orbit%species)
 xx0 = ele%value(thickness$) / (x0 / density)
-sigma = S2 * z * sqrt(xx0) / (orbit%p0c * orbit%beta) * (1.0_rp + epsilon * log10(xx0*z*z/orbit%beta**2))
-
-!
+sigma = S2 * z_part * sqrt(xx0) / (orbit%p0c * orbit%beta) * (1.0_rp + epsilon * log10(xx0*z_part**2/orbit%beta**2))
 
 call ran_gauss(rnd)
 orbit%vec(2) = orbit%vec(2) + rnd(1) * sigma
 orbit%vec(4) = orbit%vec(4) + rnd(2) * sigma
 
+! Energy loss
+
+I_excite = mean_excitation_energy_over_z(z_material) * z_material
+mass_material = mass_of(material) / atomic_mass_unit
+elec_density = N_avogadro * density * z_material / (1.0e-3_rp * mass_material) ! number_electrons / m^3
+
+n_step = nint(ele%value(num_steps$))
+do i = 1, n_step
+ p2 = (orbit%p0c * (1.0_rp + orbit%vec(6)) / mass_of(orbit%species))**2 ! beta^2 / (1 - beta^2) term
+ dE_dx = -fourpi * mass_of(electron$) * elec_density * (z_part * r_e / orbit%beta)**2 * &
+ (log(2 * mass_of(electron$) * p2 / I_excite) - orbit%beta**2)
+ E0 = orbit%p0c * (1.0_rp + orbit%vec(6)) / orbit%beta
+ E1 = E0 - dE_dx * ele%value(thickness$) / n_step
+ call convert_total_energy_to(E1, orbit%species, beta = orbit%beta, pc = pc)
+ orbit%vec(6) = (pc - orbit%p0c) / orbit%p0c
+enddo
+
 !
 
 if (logic_option(.false., make_matrix)) then

bmad/modules/attribute_mod.f90

@@ -689,7 +689,7 @@ subroutine init_attribute_name_array ()
  call init_attribute_name1 (i, l$, 'L')
 
  if (i == converter$) cycle
- if (i == foil$) cycle
+ if (i == foil$)  cycle
 
  call init_attribute_name1 (i, symplectify$, 'SYMPLECTIFY')
  call init_attribute_name1 (i, taylor_map_includes_offsets$, 'TAYLOR_MAP_INCLUDES_OFFSETS')
@@ -936,8 +936,9 @@ subroutine init_attribute_name_array ()
 call init_attribute_name1 (converter$, angle_out_max$, 'ANGLE_OUT_MAX')
 call init_attribute_name1 (converter$, species_out$, 'SPECIES_OUT')
 
-call init_attribute_name1 (foil$, thickness$, 'THICKNESS')
-call init_attribute_name1 (foil$, material_type$, 'MATERIAL_TYPE')
+call init_attribute_name1 (foil$, thickness$, 'THICKNESS')
+call init_attribute_name1 (foil$, material_type$, 'MATERIAL_TYPE')
+call init_attribute_name1 (foil$, final_charge$, 'FINAL_CHARGE')
 
 call init_attribute_name1 (lens$, l$, 'L')
 call init_attribute_name1 (lens$, radius$, 'RADIUS')
@@ -1911,7 +1912,7 @@ function attribute_type (attrib_name, ele) result (attrib_type)
 
 case ('CARTESIAN_MAP', 'CYLINDRICAL_MAP', 'FIELD_OVERLAPS', 'GEN_GRAD_MAP', 'GRID_FIELD', 'REF_ORBIT', &
  'SUPERIMPOSE', 'H_MISALIGN', 'DISPLACEMENT', 'SEGMENTED', 'PIXEL', 'TERM', 'ENERGY_PROBABILITY_CURVE', &
- 'VAR', 'WALL', 'AMP_VS_TIME', 'FREQUENCIES', 'X_KNOT', 'SR_WAKE', 'LR_WAKE', 'CURVATURE')
+ 'VAR', 'WALL', 'AMP_VS_TIME', 'FREQUENCIES', 'X_KNOT', 'SR_WAKE', 'LR_WAKE', 'CURVATURE', 'REFLECTIVITY_TABLE')
  attrib_type = is_struct$
 
 case default
@@ -2043,7 +2044,7 @@ function attribute_units (attrib_name, unrecognized_units) result (attrib_units)
  'PC_OUT_MIN', 'PC_OUT_MAX', 'E_TOT_STRONG')
  attrib_units = 'eV'
 
-case ('DELTA_REF_TIME', 'REF_TIME', 'T', 'T_OFFSET', 'DELTA_TIME', 'DT_MAX')
+case ('DELTA_REF_TIME', 'REF_TIME', 'REF_TIME_START', 'T', 'T_OFFSET', 'DELTA_TIME', 'DT_MAX')
  attrib_units = 'sec'
 
 case ('EMITTANCE_A', 'EMITTANCE_B', 'EMITTANCE_Z')

bmad/modules/bmad_struct.f90

@@ -1627,7 +1627,7 @@ module bmad_struct
 integer, parameter :: tilt$ = 2, roll$ = 2, n_part$ = 2, inherit_from_fork$ = 2 ! Important: tilt$ = roll$
 integer, parameter :: ref_tilt$ = 3, direction$ = 3, repetition_frequency$ = 3
 integer, parameter :: kick$ = 3, x_gain_err$ = 3, taylor_order$ = 3, r_solenoid$ = 3
-integer, parameter :: k1$ = 4, kx$ = 4, harmon$ = 4, h_displace$ = 4, y_gain_err$ = 4
+integer, parameter :: k1$ = 4, kx$ = 4, harmon$ = 4, h_displace$ = 4, y_gain_err$ = 4, final_charge$ = 4
 integer, parameter :: critical_angle_factor$ = 4, tilt_corr$ = 4, ref_coords$ = 4, dt_max$ = 4
 integer, parameter :: graze_angle$ = 5, k2$ = 5, b_max$ = 5, v_displace$ = 5, gradient_tot$ = 5
 integer, parameter :: ks$ = 5, flexible$ = 5, crunch$ = 5, ref_orbit_follows$ = 5, pc_out_min$ = 5

bmad/parsing/set_ele_defaults.f90

@@ -153,6 +153,10 @@ subroutine set_ele_defaults (ele, do_allocate)
  ele%value(upstream_coord_dir$) = 1
  ele%value(downstream_coord_dir$) = 1
 
+case (foil$)
+ ele%value(num_steps$) = 10
+ ele%value(final_charge$) = real_garbage$
+
 case (girder$)
  ele%value(origin_ele_ref_pt$) = center_pt$
 

regression_tests/scripts/svn_regression.py

@@ -7,7 +7,6 @@
 
 import subprocess
 svn_repo = {
- "https://accserv.classe.cornell.edu/svn/trunk/src/open_spacecharge" : "50675",
  "https://accserv.classe.cornell.edu/svn/packages/PGPLOT" : "55094", # Change from 50204
  "https://accserv.classe.cornell.edu/svn/packages/fftw" : "54912", # Change from 54142
  "https://accserv.classe.cornell.edu/svn/packages/fgsl" : "55230", # Change from 54914
@@ -19,8 +18,6 @@
  "https://accserv.classe.cornell.edu/svn/packages/openmpi" : "51988",
  "https://accserv.classe.cornell.edu/svn/packages/plplot" : "55075", # Change from 54925
  "https://accserv.classe.cornell.edu/svn/packages/xraylib" : "54916", # Change from 53181
- "https://accserv.classe.cornell.edu/svn/trunk/util" : "55221", # Change from 55177
- "https://accserv.classe.cornell.edu/svn/trunk/build_system" : "54615",
 }
 
 all_pass = True

sim_utils/interfaces/particle_species_mod.f90

@@ -140,6 +140,39 @@ module particle_species_mod
  6.4176, 6.3688, 6.2899, 6.1907, 6.0651, 6.2833, 6.1868, 6.1477, 6.0560, 6.0726, &
  5.9319, 5.9990]
 
+! Mean excitation energy (in eV) normalized by atomic Z. 
+! This is used in the Bethe-Bloch formula for calculating dE/dx of charged particles through matter.
+! Data from: https://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html
+!
+! This is suplemented with values from:
+! Hans Bichsel,
+! "Stopping Power of Fast Charged Particles in Heavy Elements",
+! NIST NISTIR 4550, April 1991.
+! Values in this paper are from the b_e column in table 1 (pg 34) and are for 19 element in the range Z = 57 to 92.
+
+real(rp), parameter :: mean_excitation_energy_over_z(92) = [ &
+ 19.20, 20.90, 13.33, 15.93, 15.20, 13.00, 11.71, 11.88, 12.78, 13.70, &
+ 13.55, 13.00, 12.77, 12.36, 11.53, 11.25, 10.24, 10.44, 10.00, 9.55, &
+ 10.29, 10.59, 10.65, 10.71, 10.88, 11.00, 11.00, 11.11, 11.10, 11.00, &
+ 10.77, 10.94, 10.52, 10.24, 9.80, 9.78, 9.81, 9.63, 9.72, 9.82, &
+ 10.17, 10.10, 9.95, 10.02, 9.98, 10.22, 10.00, 9.77, 9.96, 9.76, &
+ 9.55, 9.33, 9.26, 8.93, 8.87, 8.77, 8.32, 8.76, 8.64, 9.10, &
+ 9.18, 9.05, 9.21, 8.83, 9.45, 9.17, 9.55, 9.56, 9.77, 9.66, &
+ 9.77, 9.32, 10.05, 10.53, 9.81, 9.82, 10.23, 10.08, 10.00, 10.00, &
+ 10.00, 9.50, 8.98, 9.88, 9.71, 9.23, 9.51, 9.39, 9.45, 8.51, &
+ 9.65, 9.09]
+
+! Table from NIST before Bichsel substitution
+! 19.20, 20.90, 13.33, 15.93, 15.20, 13.00, 11.71, 11.88, 12.78, 13.70, &
+! 13.55, 13.00, 12.77, 12.36, 11.53, 11.25, 10.24, 10.44, 10.00, 9.55, &
+! 10.29, 10.59, 10.65, 10.71, 10.88, 11.00, 11.00, 11.11, 11.10, 11.00, &
+! 10.77, 10.94, 10.52, 10.24, 9.80, 9.78, 9.81, 9.63, 9.72, 9.82, &
+! 10.17, 10.10, 9.95, 10.02, 9.98, 10.22, 10.00, 9.77, 9.96, 9.76, &
+! 9.55, 9.33, 9.26, 8.93, 8.87, 8.77, 8.79, 9.02, 9.07, 9.10, &
+! 9.18, 9.26, 9.21, 9.23, 9.45, 9.52, 9.70, 9.68, 9.77, 9.77, &
+! 9.77, 9.79, 9.84, 9.82, 9.81, 9.82, 9.83, 10.13, 10.00, 10.00, &
+! 10.00, 10.04, 9.92, 9.88, 9.71, 9.23, 9.51, 9.39, 9.45, 9.41, &
+! 9.65, 9.67]
 
 ! Isotopes from NIST data 2020/Jan.
 ! The first number in a row is the average mass for natually occuring isotope mixtures.

tao/code/tao_data_and_eval_mod.f90

@@ -3148,7 +3148,7 @@ recursive subroutine tao_evaluate_a_datum (datum, u, tao_lat, datum_value, valid
  endif
 
  j = index('abc', data_type(10:10))
- if (j == 0 .or. len_trim(data_type) > 11) then
+ if (j == 0) then
  call tao_set_invalid (datum, 'DATA_TYPE = "' // trim(data_type) // '" IS NOT VALID', why_invalid, .true.)
  return
  endif

tao/version/tao_version_mod.f90

@@ -6,5 +6,5 @@
 !-
 
 module tao_version_mod
-character(*), parameter :: tao_version_date = "2023/09/28"
+character(*), parameter :: tao_version_date = "2023/09/29 01:32:20"
 end module