Proton transport

mcdc/constant.py

@@ -94,6 +94,7 @@
 # Miscellanies
 EVENT_TIME_CENSUS = 1 << 5
 EVENT_TIME_BOUNDARY = 1 << 6
+EVENT_CSDA_EDEP = 1 << 7
 
 # Materials
 MATERIAL = 0
@@ -115,6 +116,10 @@
 ELECTRON_REACTION_IONIZATION = 104
 ELECTRON_REACTION_BREMSSTRAHLUNG = 105
 ELECTRON_REACTION_EXCITATION = 106
+PROTON_REACTION_TOTAL = 200
+PROTON_REACTION_ELASTIC_SCATTERING = 201
+PROTON_REACTION_CAPTURE = 202
+PROTON_REACTION_INELASTIC_SCATTERING = 203
 
 # Particle types
 PARTICLE_NEUTRON = 0
@@ -138,12 +143,12 @@
 DISTRIBUTION_TABULATED_ENERGY_ANGLE = 8
 DISTRIBUTION_N_BODY = 9
 
-# Anguler distribution type
+# Angular distribution type
 ANGLE_ISOTROPIC = 0
 ANGLE_DISTRIBUTED = 1
 ANGLE_ENERGY_CORRELATED = 2
 
-# Referance frame
+# Reference frame
 REFERENCE_FRAME_LAB = 0
 REFERENCE_FRAME_COM = 1
 
@@ -193,8 +198,10 @@
 LIGHT_SPEED = 2.99792458e10  # cm/s
 NEUTRON_MASS = 939.565413e6  # eV/c^2
 ELECTRON_MASS = 510.99895069e3  # eV/c^2
+PROTON_MASS = 938.27208943e6  # eV/c^2
 BOLTZMANN_K = 8.61733326e-5  # eV/K
 ELECTRON_CUTOFF_ENERGY = 100  # eV
+PROTON_CUTOFF_ENERGY = 250000  # eV
 MU_CUTOFF = 0.999999
 THERMAL_THRESHOLD_FACTOR = 400
 

mcdc/main.py

@@ -168,6 +168,10 @@ def preparation():
             if isinstance(material, Material):
                 update_fissionable_from_nuclides(material)
 
+    if settings.proton_transport:
+        for nuclide in simulationPy.nuclides:
+            nuclide.set_proton_data()
+
     if settings.electron_transport:
         for element in simulationPy.elements:
             element.set_electron_data()

mcdc/mcdc_get/__init__.py

@@ -82,10 +82,18 @@
 
 import mcdc.mcdc_get.neutron_inelastic_scattering_reaction as neutron_inelastic_scattering_reaction
 
+import mcdc.mcdc_get.proton_capture_reaction as proton_capture_reaction
+
+import mcdc.mcdc_get.proton_elastic_scattering_reaction as proton_elastic_scattering_reaction
+
+import mcdc.mcdc_get.proton_inelastic_scattering_reaction as proton_inelastic_scattering_reaction
+
 import mcdc.mcdc_get.collision_data as collision_data
 
 import mcdc.mcdc_get.particle_bank as particle_bank
 
+import mcdc.mcdc_get.proton_reaction as proton_reaction
+
 import mcdc.mcdc_get.settings as settings
 
 import mcdc.mcdc_get.global_weight_roulette as global_weight_roulette

mcdc/mcdc_get/native_material.py

@@ -117,3 +117,61 @@ def element_densities_chunk(start, length, native_material, data):
     start += native_material["element_densities_offset"]
     end = start + length
     return data[start:end]
+
+
+@njit
+def stopping_power(index, native_material, data):
+    offset = native_material["stopping_power_offset"]
+    return data[offset + index]
+
+
+@njit
+def stopping_power_all(native_material, data):
+    start = native_material["stopping_power_offset"]
+    size = native_material["stopping_power_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def stopping_power_last(native_material, data):
+    start = native_material["stopping_power_offset"]
+    size = native_material["stopping_power_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def stopping_power_chunk(start, length, native_material, data):
+    start += native_material["stopping_power_offset"]
+    end = start + length
+    return data[start:end]
+
+
+@njit
+def stopping_power_energy_grid(index, native_material, data):
+    offset = native_material["stopping_power_energy_grid_offset"]
+    return data[offset + index]
+
+
+@njit
+def stopping_power_energy_grid_all(native_material, data):
+    start = native_material["stopping_power_energy_grid_offset"]
+    size = native_material["stopping_power_energy_grid_length"]
+    end = start + size
+    return data[start:end]
+
+
+@njit
+def stopping_power_energy_grid_last(native_material, data):
+    start = native_material["stopping_power_energy_grid_offset"]
+    size = native_material["stopping_power_energy_grid_length"]
+    end = start + size
+    return data[end - 1]
+
+
+@njit
+def stopping_power_energy_grid_chunk(start, length, native_material, data):
+    start += native_material["stopping_power_energy_grid_offset"]
+    end = start + length
+    return data[start:end]