import math
import numpy as np
import numbers
from typing import Optional, List
from . import ray, material, source, shield, detector
import importlib
pyvista_spec = importlib.util.find_spec("pyvista")
pyvista_found = pyvista_spec is not None
if pyvista_found:
import pyvista
''' '''
'''
ZapMeNot - a point kernel photon shielding library
Copyright (C) 2019-2025 C. Alan Ford
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <https://www.gnu.org/licenses/>.
'''
[docs]
class Model:
"""Performs point-kernel shielding analysis.
The Model class combines various shielding elements to perform
the point-kernel photon shielding analysis. These elements include
sources, shields, and detectors.
"""
# used to calculate exposure (R/sec) from flux (photon/cm2 sec),
# photon energy (MeV),
# and linear energy absorption coeff (cm2/g)
# aka, "flux to exposure conversion factor"
# for more information, see "Radiation Shielding", J. K. Shultis
# and R.E. Faw, 2000, page 141.
# This value is based on a value of energy deposition
# per ion in air of 33.85 [ICRU Report 39, 1979].
FLUX_TO_EXPOSURE_CONVERSION_FACTOR = 1.835E-8
def __init__(self) -> None:
self.source: Optional[source.Source] = None
self.shield_list: List[shield.Shield] = []
self.detector: Optional[detector.Detector] = None
self.filler_material: Optional[material.Material] = None
self.buildup_factor_material: Optional[material.Material] = None
[docs]
def set_filler_material(self, filler_material: str,
density: Optional[float] = None) -> None:
r"""Set the filler material used by the model
Parameters
----------
filler_material
The material to be used.
density
The density of the material in g/cm\ :sup:`3`.
"""
if not isinstance(filler_material, str):
raise ValueError("Invalid filler material")
self.filler_material = material.Material(filler_material)
if density is not None:
if not isinstance(density, numbers.Number):
raise ValueError(f"Invalid density: {density}")
self.filler_material.density = density
[docs]
def add_source(self, new_source: source.Source) -> None:
"""Set the source used by the model. Only one source instance
is allowed per model. Source can be replaced by calling this method.
Parameters
----------
new_source
The source to be used.
"""
if not isinstance(new_source, source.Source):
raise ValueError("Invalid source")
if self.source is not None and isinstance(self.source, shield.Shield):
# removing the old source from the shield list
if self.source in self.shield_list:
self.shield_list.remove(self.source)
self.source = new_source
# don't forget that sources are shields too!
if isinstance(new_source, shield.Shield):
self.shield_list.append(new_source)
[docs]
def add_shield(self, new_shield: shield.Shield) -> None:
"""Add a shield to the collection of shields used by the model.
Parameters
----------
new_shield
The shield to be added.
"""
if not isinstance(new_shield, shield.Shield):
raise ValueError("Invalid shield")
self.shield_list.append(new_shield)
[docs]
def add_detector(self, new_detector: detector.Detector) -> None:
"""Set the detector used by the model. Only one detector instance
is allowed per model. Detector can be replaced by calling this method.
Parameters
----------
new_detector
The detector to be used in the model.
"""
if not isinstance(new_detector, detector.Detector):
raise ValueError("Invalid detector")
self.detector = new_detector
[docs]
def set_buildup_factor_material(self, new_material: material.Material) \
-> None:
"""Set the material used to calculation exposure buildup factors.
Parameters
----------
new_material
The material to be used in buildup factor calculations.
"""
if not isinstance(new_material, material.Material):
raise ValueError("Invalid buildup factor material")
self.buildup_factor_material = new_material
[docs]
def calculate_exposure(self) -> float:
"""Calculates the exposure at the detector location.
Note: Significant use of Numpy arrays to speed up evaluating the
dose from each source point. A "for loop" is used to loop
through photon energies, but many of the iterations through
all source points is performed using matrix math.
Returns
-------
The exposure in units of mR/hr.
"""
results_by_photon_energy = self.generate_summary()
if len(results_by_photon_energy) == 0:
return 0 # may occur if source has no photons
elif len(results_by_photon_energy) == 1:
return results_by_photon_energy[0][4] # mR/hr
else:
# sum exposure over all photons
an_array = np.array(results_by_photon_energy)
integral_results = np.sum(an_array[:, 4])
return integral_results # mR/hr
[docs]
def generate_summary(self) -> List[List[float]]:
"""Calculates the energy flux and exposure at the detector location.
Note: Significant use of Numpy arrays to speed up evaluating the
dose from each source point. A "for loop" is used to loop
through photon energies, but many of the iterations through
all source points is performed using matrix math.
Returns
-------
List, by photon energy, of photon energy, photon emission rate,
uncollided energy flux, uncollided exposure, and total exposure
"""
# build an array of shield crossing lengths.
# The first index is the source point.
# The second index is the shield (including the source body).
# The total transit distance in the "filler" material (if any)
# is determined by subtracting the sum of the shield crossing
# lengths from the total ray length.
if self.source is None:
raise ValueError("Model is missing a source")
if self.detector is None:
raise ValueError("Model is missing a detector")
source_points = self.source._get_source_points()
source_point_weights = self.source._get_source_point_weights()
crossing_distances = np.zeros((len(source_points),
len(self.shield_list)))
total_distance = np.zeros((len(source_points)))
for index, nextPoint in enumerate(source_points):
vector = ray.FiniteLengthRay(nextPoint,
self.detector.location)
total_distance[index] = vector._length
# check to see if source point and detector are coincident
if total_distance[index] == 0.0:
raise ValueError("detector and source are coincident")
for index2, currentShield in enumerate(self.shield_list):
crossing_distances[index, index2] = \
currentShield._get_crossing_length(vector)
gaps = total_distance - np.sum(crossing_distances, axis=1)
if np.amin(gaps) < 0:
raise ValueError("Looks like shields and/or sources overlap")
results_by_photon_energy = []
# get a list of photons (energy & intensity) from the source
spectrum = self.source.get_photon_source_list()
air = material.Material('air')
# iterate through the photon list
for photon in spectrum:
photon_energy = photon[0]
# photon source strength
photon_yield = photon[1]
dose_coeff = air.get_mass_energy_abs_coeff(photon_energy)
# determine the xsecs
xsecs = np.zeros((len(self.shield_list)))
for index, currentShield in enumerate(self.shield_list):
xsecs[index] = currentShield.material.density * \
currentShield.material.get_mass_atten_coeff(photon_energy)
# determine an array of mean free paths, one per source point
total_mfp = crossing_distances * xsecs
total_mfp = np.sum(total_mfp, axis=1)
# add the gaps if required
if self.filler_material is not None:
gap_xsec = self.filler_material.density * \
self.filler_material.get_mass_atten_coeff(photon_energy)
total_mfp = total_mfp + (gaps * gap_xsec)
uncollided_flux_factor = np.exp(-total_mfp)
if (self.buildup_factor_material is not None):
buildup_factor = \
self.buildup_factor_material.get_buildup_factor(
photon_energy, total_mfp)
else:
buildup_factor = 1.0
# Notes for the following code:
# uncollided_point_energy_flux - an ARRAY of uncollided energy
# flux for a at the detector from a range of quadrature
# locations and a specific photon energy
# total_uncollided_energy_flux - an INTEGRAL of uncollided energy
# flux for a at the detector and a specific photon energy
#
uncollided_point_energy_flux = photon_yield * \
np.asarray(source_point_weights) \
* uncollided_flux_factor * photon_energy * \
(1/(4*math.pi*np.power(total_distance, 2)))
total_uncollided_energy_flux = np.sum(uncollided_point_energy_flux)
uncollided_point_exposure = uncollided_point_energy_flux * \
Model.FLUX_TO_EXPOSURE_CONVERSION_FACTOR * dose_coeff * \
1000 * 3600 # mR/hr
total_uncollided_exposure = np.sum(uncollided_point_exposure)
collided_point_exposure = uncollided_point_exposure * \
buildup_factor
total_collided_exposure = np.sum(collided_point_exposure)
results_by_photon_energy.append(
[photon_energy, photon_yield, total_uncollided_energy_flux,
total_uncollided_exposure, total_collided_exposure])
return results_by_photon_energy
[docs]
def display(self) -> None:
"""
Produces a graphic display of the model.
"""
if pyvista_found:
# find the bounding box for all objects
bounds = self._findBoundingBox()
pl: pyvista.Plotter = pyvista.Plotter()
self._trimBlocks(pl, bounds)
self._addPoints(pl)
pl.show_bounds(grid='front', location='outer', all_edges=True)
if self.source is not None or self.detector is not None:
pl.add_legend(face=None, size=(0.1, 0.1))
pl.show()
def _trimBlocks(self, pl: pyvista.Plotter, bounds: List[float]) -> None:
"""
Adds shields to a Plotter instance after trimming any
infinite shields to a predefined bounding box.
"""
# display color for shield geometries
SHIELD_COLOR = 'blue'
# display color for source geometry
SOURCE_COLOR = 'red'
for currentShield in self.shield_list:
opacity = 1.0
if currentShield.is_hollow():
opacity = 0.5
# first handle infinite shields
if isinstance(currentShield, shield.SemiInfiniteShield):
clipped = currentShield.draw()
if isinstance(clipped, pyvista.PolyData):
clipped = clipped.clip_closed_surface(
normal='x', origin=[bounds[0], 0, 0])
clipped = clipped.clip_closed_surface( # type: ignore
normal='y', origin=[0, bounds[2], 0])
clipped = clipped.clip_closed_surface( # type: ignore
normal='z', origin=[0, 0, bounds[4]])
clipped = clipped.clip_closed_surface( # type: ignore
normal='-x', origin=[bounds[1], 0, 0])
clipped = clipped.clip_closed_surface( # type: ignore
normal='-y', origin=[0, bounds[3], 0])
clipped = clipped.clip_closed_surface( # type: ignore
normal='-z', origin=[0, 0, bounds[5]])
pl.add_mesh(clipped, color=SHIELD_COLOR, opacity=opacity)
# now handle the sources and non-infinite shields
else:
if isinstance(currentShield, source.Source):
if isinstance(self.source, source.PointSource):
# point sources are handled later
pass
else:
pl.add_mesh(currentShield.draw(),
SOURCE_COLOR, label='source', line_width=3)
else:
pl.add_mesh(currentShield.draw(), SHIELD_COLOR,
opacity=opacity)
# now add the "bounds" as a transparent block to for a display size
mesh = pyvista.Box(bounds)
pl.add_mesh(mesh, opacity=0)
def _findBoundingBox(self) -> List[float]:
"""Calculates a bounding box is X, Y, Z geometry that
includes the volumes of all shields, the source, and the detector
"""
blocks = pyvista.MultiBlock()
for currentShield in self.shield_list:
if isinstance(currentShield, shield.SemiInfiniteShield):
# for infinite shield bodies,
# project the detector location onto the infinite surface
# to get points to add to the geometry
if isinstance(self.detector, detector.Detector):
points = currentShield._projection(self.detector.x,
self.detector.y,
self.detector.z)
for point in points:
# we are appending a degenerate line as a
# representation of a point
blocks.append(pyvista.Line(point, point))
else:
# add finite shields to the MultiBlock composite
blocks.append(currentShield.draw())
# include the detector geometry in the MultiBlock composite
if self.detector is not None:
blocks.append(self.detector.draw())
# check for a zero width bounding box in any direction
xmin, xmax, ymin, ymax, zmin, zmax = blocks.bounds
x_width = abs(xmax - xmin)
y_width = abs(ymax - ymin)
z_width = abs(zmax - zmin)
max_width = max(x_width, y_width, z_width)
# define a minimum dimension as 20% of the maximum dimension
min_width = max_width * 0.20
# check for dimensions smaller than the defined minimum
if x_width < min_width:
xmin = xmin - min_width/2
xmax = xmax + min_width/2
if y_width < min_width:
ymin = ymin - min_width/2
ymax = ymax + min_width/2
if z_width < min_width:
zmin = zmin - min_width/2
zmax = zmax + min_width/2
# increase the display bounds by a smidge to avoid
# inadvertent clipping
bounds = [xmin, xmax, ymin, ymax, zmin, zmax]
boundingBox = [x * 1.01 for x in bounds]
return boundingBox
def _addPoints(self, pl: pyvista.Plotter) -> None:
"""
the goal here is to add 'points' to the display, but they
must be represented as spheres to have some physical
volume to display. Points will be displayed with a radius
of 5% of the smallest dimension of the bounding box.
A problem can occur if the bounding box has a width of 0 in one
or more of three dimensions. An exception is thrown if bounds
in all three directions are of zero width. Otherwise the zero
is ignored and the next largest dimension is used to size the
point representation.
"""
# ratio of the radius of point objects to the minimum width
# of other objects to be displayed
POINT_RATIO = 0.025
# display color for source geometry
SOURCE_COLOR = 'red'
# display color for detector geometry
DETECTOR_COLOR = 'green'
widths = [abs(pl.bounds[1] - pl.bounds[0]),
abs(pl.bounds[3] - pl.bounds[2]),
abs(pl.bounds[5] - pl.bounds[4])]
good_widths = [width for width in widths if width > 0]
if len(good_widths) == 0:
raise ValueError("detector and source are coincident or missing")
point_radius = min(good_widths) * POINT_RATIO
# check if the source is a point source
if isinstance(self.source, source.PointSource):
body = pyvista.Sphere(center=(self.source._x,
self.source._y,
self.source._z),
radius=point_radius)
pl.add_mesh(
body, line_width=5, color=SOURCE_COLOR,
label='source')
if self.detector is not None:
body = pyvista.Sphere(center=(self.detector.x,
self.detector.y,
self.detector.z),
radius=point_radius)
pl.add_mesh(
body, line_width=5, color=DETECTOR_COLOR,
label='detector')
# pl.set_background(color='white')