diff --git a/README.md b/README.md index d3a79ea..3b78136 100644 --- a/README.md +++ b/README.md @@ -31,6 +31,14 @@ With the virtual environment activated, run: pip install -e .[dev] ``` +## Running example landscape + +The example landscape can be generated using the command-line interface. Still in the virtual environment, run + +``` +pgrad --test --loglevel DEBUG +``` + ## Tests Tests can be run with `pytest` from the root directory of the repository. With the virtual environment activated, run: diff --git a/pyproject.toml b/pyproject.toml index 0f82ac2..a5be416 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,6 +7,7 @@ where = ["src"] [tool.setuptools.package-data] "pg_rad.data" = ["*.csv"] +"pg_rad.configs" = ["*.yml"] [project] name = "pg-rad" @@ -21,14 +22,19 @@ dependencies = [ "matplotlib>=3.9.2", "numpy>=2", "pandas>=2.3.1", - "pyyaml>=6.0.2" + "pyyaml>=6.0.2", + "scipy>=1.15.0" ] license = "MIT" license-files = ["LICEN[CS]E*"] +[project.scripts] +pgrad = "pg_rad.main:main" + [project.urls] Homepage = "https://github.com/pim-n/pg-rad" Issues = "https://github.com/pim-n/pg-rad/issues" [project.optional-dependencies] + dev = ["pytest", "mkinit", "notebook", "mkdocs-material", "mkdocstrings-python", "mkdocs-jupyter"] \ No newline at end of file diff --git a/src/pg_rad/configs/filepaths.py b/src/pg_rad/configs/filepaths.py new file mode 100644 index 0000000..e8103af --- /dev/null +++ b/src/pg_rad/configs/filepaths.py @@ -0,0 +1,3 @@ +ATTENUATION_TABLE = 'attenuation_table.csv' +TEST_EXP_DATA = 'test_path_coords.csv' +LOGGING_CONFIG = 'logging.yml' diff --git a/src/pg_rad/data/__init__.py b/src/pg_rad/data/__init__.py new file mode 100644 index 0000000..a9a2c5b --- /dev/null +++ b/src/pg_rad/data/__init__.py @@ -0,0 +1 @@ +__all__ = [] diff --git a/src/pg_rad/data/attenuation_table.csv b/src/pg_rad/data/attenuation_table.csv new file mode 100644 index 0000000..0c8f024 --- /dev/null +++ b/src/pg_rad/data/attenuation_table.csv @@ -0,0 +1,39 @@ +energy_mev,mu +1.00000E-03,3.606E+03 +1.50000E-03,1.191E+03 +2.00000E-03,5.279E+02 +3.00000E-03,1.625E+02 +3.20290E-03,1.340E+02 +3.20290E-03,1.485E+02 +4.00000E-03,7.788E+01 +5.00000E-03,4.027E+01 +6.00000E-03,2.341E+01 +8.00000E-03,9.921E+00 +1.00000E-02,5.120E+00 +1.50000E-02,1.614E+00 +2.00000E-02,7.779E-01 +3.00000E-02,3.538E-01 +4.00000E-02,2.485E-01 +5.00000E-02,2.080E-01 +6.00000E-02,1.875E-01 +8.00000E-02,1.662E-01 +1.00000E-01,1.541E-01 +1.50000E-01,1.356E-01 +2.00000E-01,1.233E-01 +3.00000E-01,1.067E-01 +4.00000E-01,9.549E-02 +5.00000E-01,8.712E-02 +6.00000E-01,8.055E-02 +8.00000E-01,7.074E-02 +1.00000E+00,6.358E-02 +1.25000E+00,5.687E-02 +1.50000E+00,5.175E-02 +2.00000E+00,4.447E-02 +3.00000E+00,3.581E-02 +4.00000E+00,3.079E-02 +5.00000E+00,2.751E-02 +6.00000E+00,2.522E-02 +8.00000E+00,2.225E-02 +1.00000E+01,2.045E-02 +1.50000E+01,1.810E-02 +2.00000E+01,1.705E-02 \ No newline at end of file diff --git a/tests/data/coordinates.csv b/src/pg_rad/data/test_path_coords.csv similarity index 100% rename from tests/data/coordinates.csv rename to src/pg_rad/data/test_path_coords.csv diff --git a/src/pg_rad/exceptions/__init__.py b/src/pg_rad/exceptions/__init__.py index 426b7e7..99b24a1 100644 --- a/src/pg_rad/exceptions/__init__.py +++ b/src/pg_rad/exceptions/__init__.py @@ -4,7 +4,7 @@ __ignore__ = ["logger"] from pg_rad.exceptions import exceptions from pg_rad.exceptions.exceptions import (ConvergenceError, DataLoadError, - InvalidCSVError,) + InvalidCSVError, OutOfBoundsError,) __all__ = ['ConvergenceError', 'DataLoadError', 'InvalidCSVError', - 'exceptions'] + 'OutOfBoundsError', 'exceptions'] diff --git a/src/pg_rad/exceptions/exceptions.py b/src/pg_rad/exceptions/exceptions.py index 60a1b34..dc8c3eb 100644 --- a/src/pg_rad/exceptions/exceptions.py +++ b/src/pg_rad/exceptions/exceptions.py @@ -8,3 +8,7 @@ class DataLoadError(Exception): class InvalidCSVError(DataLoadError): """Raised when a file is not a valid CSV.""" + + +class OutOfBoundsError(Exception): + """Raised when an object is attempted to be placed out of bounds.""" diff --git a/src/pg_rad/isotopes/__init__.py b/src/pg_rad/isotopes/__init__.py index 0c7976d..6f000c0 100644 --- a/src/pg_rad/isotopes/__init__.py +++ b/src/pg_rad/isotopes/__init__.py @@ -2,7 +2,9 @@ __ignore__ = ["logger"] from pg_rad.isotopes import isotope +from pg_rad.isotopes import presets from pg_rad.isotopes.isotope import (Isotope,) +from pg_rad.isotopes.presets import (CS137,) -__all__ = ['Isotope', 'isotope'] +__all__ = ['CS137', 'Isotope', 'isotope', 'presets'] diff --git a/src/pg_rad/isotopes/isotope.py b/src/pg_rad/isotopes/isotope.py index ea758b8..78ac54e 100644 --- a/src/pg_rad/isotopes/isotope.py +++ b/src/pg_rad/isotopes/isotope.py @@ -1,3 +1,6 @@ +from pg_rad.physics import get_mass_attenuation_coeff + + class Isotope: """Represents the essential information of an isotope. @@ -11,13 +14,14 @@ class Isotope: name: str, E: float, b: float - ): + ): if E <= 0: raise ValueError("primary_gamma must be a positive energy (keV).") if not (0 <= b <= 1): - raise ValueError("branching_ratio_pg must be a ratio (0 <= b <= 1)") + raise ValueError("branching_ratio_pg must be a ratio b in [0,1]") self.name = name self.E = E self.b = b + self.mu_mass_air = get_mass_attenuation_coeff(E / 1000) diff --git a/src/pg_rad/isotopes/presets.py b/src/pg_rad/isotopes/presets.py new file mode 100644 index 0000000..d8135d6 --- /dev/null +++ b/src/pg_rad/isotopes/presets.py @@ -0,0 +1,10 @@ +from .isotope import Isotope + + +class CS137(Isotope): + def __init__(self): + super().__init__( + name="Cs-137", + E=661.66, + b=0.851 + ) diff --git a/src/pg_rad/landscape/__init__.py b/src/pg_rad/landscape/__init__.py index 9a0c85c..9dda6f9 100644 --- a/src/pg_rad/landscape/__init__.py +++ b/src/pg_rad/landscape/__init__.py @@ -1,8 +1,11 @@ # do not expose internal logger when running mkinit __ignore__ = ["logger"] +from pg_rad.landscape import director from pg_rad.landscape import landscape -from pg_rad.landscape.landscape import (Landscape, create_landscape_from_path,) +from pg_rad.landscape.director import (LandscapeDirector,) +from pg_rad.landscape.landscape import (Landscape, LandscapeBuilder,) -__all__ = ['Landscape', 'create_landscape_from_path', 'landscape'] +__all__ = ['Landscape', 'LandscapeBuilder', 'LandscapeDirector', 'director', + 'landscape'] diff --git a/src/pg_rad/landscape/director.py b/src/pg_rad/landscape/director.py new file mode 100644 index 0000000..edddd7f --- /dev/null +++ b/src/pg_rad/landscape/director.py @@ -0,0 +1,23 @@ +from importlib.resources import files +import logging + +from pg_rad.configs.filepaths import TEST_EXP_DATA +from pg_rad.isotopes import CS137 +from pg_rad.landscape.landscape import LandscapeBuilder +from pg_rad.objects import PointSource + + +class LandscapeDirector: + def __init__(self): + self.logger = logging.getLogger(__name__) + self.logger.debug("LandscapeDirector initialized.") + + def build_test_landscape(self): + fp = files('pg_rad.data').joinpath(TEST_EXP_DATA) + source = PointSource(activity=100E9, isotope=CS137(), pos=(0, 0, 0)) + lb = LandscapeBuilder("Test landscape") + lb.set_air_density(1.243) + lb.set_path_from_experimental_data(fp, z=0) + lb.set_point_sources(source) + landscape = lb.build() + return landscape diff --git a/src/pg_rad/landscape/landscape.py b/src/pg_rad/landscape/landscape.py index 1883694..805081e 100644 --- a/src/pg_rad/landscape/landscape.py +++ b/src/pg_rad/landscape/landscape.py @@ -1,133 +1,161 @@ import logging +from typing import Self -from matplotlib import pyplot as plt -from matplotlib.patches import Circle -import numpy as np - -from pg_rad.path import Path +from pg_rad.dataloader import load_data +from pg_rad.exceptions import OutOfBoundsError from pg_rad.objects import PointSource +from pg_rad.path import Path, path_from_RT90 +from pg_rad.physics.fluence import phi_single_source logger = logging.getLogger(__name__) class Landscape: - """A generic Landscape that can contain a Path and sources. - - Args: - air_density (float, optional): Air density, kg/m^3. Defaults to 1.243. - size (int | tuple[int, int, int], optional): Size of the world. - Defaults to 500. - scale (str, optional): The scale of the size argument passed. - Defaults to 'meters'. + """ + A generic Landscape that can contain a Path and sources. """ def __init__( self, - air_density: float = 1.243, - size: int | tuple[int, int, int] = 500, - scale: str = 'meters' + name: str, + path: Path, + point_sources: list[PointSource] = [], + size: tuple[int, int, int] = [500, 500, 50], + air_density: float = 1.243 ): - if isinstance(size, int): - self.world = np.zeros((size, size, size)) - elif isinstance(size, tuple) and len(size) == 3: - self.world = np.zeros(size) - else: - raise TypeError("size must be integer or a tuple of 3 integers.") - - self.air_density = air_density - self.scale = scale - - self.path: Path = None - self.sources: list[PointSource] = [] - logger.debug("Landscape initialized.") - - def plot(self, z: float | int = 0): - """Plot a slice of the world at a height `z`. + """Initialize a landscape. Args: - z (int, optional): Height of slice. Defaults to 0. + path (Path): A Path object. + point_sources (list[PointSource], optional): List of point sources. + air_density (float, optional): Air density in kg/m^3. + Defaults to 1.243. + size (tuple[int, int, int], optional): (x,y,z) dimensions of world + in meters. Defaults to [500, 500, 50]. - Returns: - fig, ax: Matplotlib figure objects. + Raises: + TypeError: _description_ """ - x_lim, y_lim, _ = self.world.shape - fig, ax = plt.subplots() - ax.set_xlim(right=x_lim) - ax.set_ylim(top=y_lim) - ax.set_xlabel(f"X [{self.scale}]") - ax.set_ylabel(f"Y [{self.scale}]") + self.name = name + self.path = path + self.point_sources = point_sources + self.size = size + self.air_density = air_density - if self.path is not None: - ax.plot(self.path.x_list, self.path.y_list, 'bo-') + logger.debug(f"Landscape created: {self.name}") - for s in self.sources: - if np.isclose(s.z, z): - dot = Circle( - (s.x, s.y), - radius=5, - color=s.color, - zorder=5 + def calculate_fluence_at(self, pos: tuple): + total_phi = 0. + for source in self.point_sources: + r = source.distance_to(pos) + phi_source = phi_single_source( + r=r, + activity=source.activity, + branching_ratio=source.isotope.b, + mu_mass_air=source.isotope.mu_mass_air, + air_density=self.air_density + ) + total_phi += phi_source + return total_phi + + def calculate_fluence_along_path(self): + pass + + +class LandscapeBuilder: + def __init__(self, name: str = "Unnamed landscape"): + self.name = name + self._path = None + self._point_sources = [] + self._size = None + self._air_density = None + + logger.debug(f"LandscapeBuilder initialized: {self.name}") + + def set_air_density(self, air_density) -> Self: + """Set the air density of the world.""" + self._air_density = air_density + return self + + def set_landscape_size(self, size: tuple[int, int, int]) -> Self: + """Set the size of the landscape in meters (x,y,z).""" + if self._path and any(p > s for p, s in zip(self._path.size, size)): + raise OutOfBoundsError( + "Cannot set landscape size smaller than the path." + ) + + self._size = size + logger.debug("Size of the landscape has been updated.") + + return self + + def set_path_from_experimental_data( + self, + filename: str, + z: int, + east_col: str = "East", + north_col: str = "North" + ) -> Self: + df = load_data(filename) + self._path = path_from_RT90( + df=df, + east_col=east_col, + north_col=north_col, + z=z + ) + + # The size of the landscape will be updated if + # 1) _size is not set, or + # 2) _size is too small to contain the path. + needs_resize = ( + not self._size + or any(p > s for p, s in zip(self._path.size, self._size)) + ) + + if needs_resize: + if not self._size: + logger.debug("Because no Landscape size was set, " + "it will now set to path dimensions.") + else: + logger.warning( + "Path exceeds current landscape size. " + "Landscape size will be expanded to accommodate path." ) - ax.text( - s.x + 0.06, - s.y + 0.06, - s.name, - color=s.color, - fontsize=10, - ha="left", - va="bottom", - zorder=6 - ) + self.set_landscape_size(self._path.size) - ax.add_patch(dot) + return self - return fig, ax - - def add_sources(self, *sources: PointSource): + def set_point_sources(self, *sources): """Add one or more point sources to the world. Args: *sources (pg_rad.sources.PointSource): One or more sources, passed as Source1, Source2, ... Raises: - ValueError: If the source is outside the boundaries of the + OutOfBoundsError: If any source is outside the boundaries of the landscape. """ - max_x, max_y, max_z = self.world.shape[:3] if any( - not (0 <= source.x < max_x and - 0 <= source.y < max_y and - 0 <= source.z < max_z) + any(p < 0 or p >= s for p, s in zip(source.pos, self._size)) for source in sources ): - raise ValueError("One or more sources are outside the landscape!") + raise OutOfBoundsError( + "One or more sources attempted to " + "be placed outside the landscape." + ) - self.sources.extend(sources) + self._point_sources = sources - def set_path(self, path: Path): - """ - Set the path in the landscape. - """ - self.path = path + def build(self): + landscape = Landscape( + name=self.name, + path=self._path, + point_sources=self._point_sources, + size=self._size, + air_density=self._air_density + ) - -def create_landscape_from_path(path: Path, max_z: float | int = 50): - """Generate a landscape from a path, using its dimensions to determine - the size of the landscape. - - Args: - path (Path): A Path object describing the trajectory. - max_z (int, optional): Height of the world. Defaults to 50 meters. - - Returns: - landscape (pg_rad.landscape.Landscape): A landscape with dimensions - based on the provided Path. - """ - max_x = np.ceil(max(path.x_list)) - max_y = np.ceil(max(path.y_list)) - - landscape = Landscape(size=(max_x, max_y, max_z)) - landscape.path = path - return landscape + logger.info(f"Landscape built successfully: {landscape.name}") + return landscape diff --git a/src/pg_rad/logger/__init__.py b/src/pg_rad/logger/__init__.py new file mode 100644 index 0000000..57b4e3f --- /dev/null +++ b/src/pg_rad/logger/__init__.py @@ -0,0 +1,5 @@ +from pg_rad.logger import logger + +from pg_rad.logger.logger import (setup_logger,) + +__all__ = ['logger', 'setup_logger'] diff --git a/src/pg_rad/logging/logger.py b/src/pg_rad/logger/logger.py similarity index 67% rename from src/pg_rad/logging/logger.py rename to src/pg_rad/logger/logger.py index a460e20..370f35f 100644 --- a/src/pg_rad/logging/logger.py +++ b/src/pg_rad/logger/logger.py @@ -1,8 +1,10 @@ -import logging -import pathlib +import logging.config +from importlib.resources import files import yaml +from pg_rad.configs.filepaths import LOGGING_CONFIG + def setup_logger(log_level: str = "WARNING"): levels = ["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"] @@ -10,8 +12,7 @@ def setup_logger(log_level: str = "WARNING"): if log_level not in levels: raise ValueError(f"Log level must be one of {levels}.") - base_dir = pathlib.Path(__file__).resolve().parent - config_file = base_dir / "configs" / "logging.yml" + config_file = files('pg_rad.configs').joinpath(LOGGING_CONFIG) with open(config_file) as f: config = yaml.safe_load(f) diff --git a/src/pg_rad/logging/__init__.py b/src/pg_rad/logging/__init__.py deleted file mode 100644 index 666f5de..0000000 --- a/src/pg_rad/logging/__init__.py +++ /dev/null @@ -1,5 +0,0 @@ -from pg_rad.logging import logger - -from pg_rad.logging.logger import (setup_logger,) - -__all__ = ['logger', 'setup_logger'] diff --git a/src/pg_rad/main.py b/src/pg_rad/main.py new file mode 100644 index 0000000..19f4161 --- /dev/null +++ b/src/pg_rad/main.py @@ -0,0 +1,33 @@ +import argparse + +from pg_rad.logger import setup_logger +from pg_rad.landscape import LandscapeDirector + + +def main(): + parser = argparse.ArgumentParser( + prog="pg-rad", + description="Primary Gamma RADiation landscape tool" + ) + + parser.add_argument( + "--test", + action="store_true", + help="Load and run the test landscape" + ) + parser.add_argument( + "--loglevel", + default="INFO", + choices=["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"], + ) + + args = parser.parse_args() + setup_logger(args.loglevel) + + if args.test: + landscape = LandscapeDirector().build_test_landscape() + print(landscape.name) + + +if __name__ == "__main__": + main() diff --git a/src/pg_rad/objects/__init__.py b/src/pg_rad/objects/__init__.py index 2170a46..92edc21 100644 --- a/src/pg_rad/objects/__init__.py +++ b/src/pg_rad/objects/__init__.py @@ -1,13 +1,11 @@ # do not expose internal logger when running mkinit __ignore__ = ["logger"] -from pg_rad.objects import detectors from pg_rad.objects import objects from pg_rad.objects import sources -from pg_rad.objects.detectors import (Detector,) from pg_rad.objects.objects import (BaseObject,) from pg_rad.objects.sources import (PointSource,) -__all__ = ['BaseObject', 'Detector', 'PointSource', 'detectors', 'objects', +__all__ = ['BaseObject', 'PointSource', 'objects', 'sources'] diff --git a/src/pg_rad/objects/objects.py b/src/pg_rad/objects/objects.py index 16f03b9..d3ec285 100644 --- a/src/pg_rad/objects/objects.py +++ b/src/pg_rad/objects/objects.py @@ -1,36 +1,42 @@ -import math from typing import Self +import numpy as np + class BaseObject: def __init__( self, - x: float, - y: float, - z: float, + pos: tuple[float, float, float] = (0, 0, 0), name: str = "Unnamed object", color: str = 'grey'): """ A generic object. Args: - x (float): X coordinate. - y (float): Y coordinate. - z (float): Z coordinate. + pos (tuple[float, float, float]): Position vector (x,y,z). name (str, optional): Name for the object. Defaults to "Unnamed object". color (str, optional): Matplotlib compatible color string. Defaults to "red". """ - self.x = x - self.y = y - self.z = z + if len(pos) != 3: + raise ValueError("Position must be tuple of length 3 (x,y,z).") + self.pos = pos self.name = name self.color = color - def distance_to(self, other: Self) -> float: - return math.dist( - (self.x, self.y, self.z), - (other.x, other.y, other.z), - ) + def distance_to(self, other: Self | tuple) -> float: + if isinstance(other, tuple) and len(other) == 3: + r = np.linalg.norm( + np.subtract(self.pos, other) + ) + else: + try: + r = np.linalg.norm( + np.subtract(self.pos, other.pos) + ) + except AttributeError as e: + raise e("other must be an object in the world \ + or a position tuple (x,y,z).") + return r diff --git a/src/pg_rad/objects/sources.py b/src/pg_rad/objects/sources.py index 1bc3a2d..9361f1f 100644 --- a/src/pg_rad/objects/sources.py +++ b/src/pg_rad/objects/sources.py @@ -11,26 +11,23 @@ class PointSource(BaseObject): def __init__( self, - x: float, - y: float, - z: float, activity: int, isotope: Isotope, + pos: tuple[float, float, float] = (0, 0, 0), name: str | None = None, - color: str = "red"): + color: str = 'red' + ): """A point source. Args: - x (float): X coordinate. - y (float): Y coordinate. - z (float): Z coordinate. activity (int): Activity A in MBq. isotope (Isotope): The isotope. - name (str | None, optional): Can give the source a unique name. - Defaults to None, making the name sequential. - (Source-1, Source-2, etc.). - color (str, optional): Matplotlib compatible color string. - Defaults to "red". + pos (tuple[float, float, float], optional): + Position of the PointSource. + name (str, optional): Can give the source a unique name. + If not provided, point sources are sequentially + named: Source1, Source2, ... + color (str, optional): Matplotlib compatible color string """ self.id = PointSource._id_counter @@ -40,11 +37,10 @@ class PointSource(BaseObject): if name is None: name = f"Source {self.id}" - super().__init__(x, y, z, name, color) + super().__init__(pos, name, color) self.activity = activity self.isotope = isotope - self.color = color logger.debug(f"Source created: {self.name}") diff --git a/src/pg_rad/path/path.py b/src/pg_rad/path/path.py index e750ba2..5c9ade3 100644 --- a/src/pg_rad/path/path.py +++ b/src/pg_rad/path/path.py @@ -42,14 +42,18 @@ class Path: def __init__( self, coord_list: Sequence[tuple[float, float]], - z: float = 0 + z: float = 0., + z_box: float = 50. ): """Construct a path of sequences based on a list of coordinates. Args: coord_list (Sequence[tuple[float, float]]): List of x,y coordinates. - z (float, optional): Height of the path. Defaults to 0. + z (float, optional): position of the path in z-direction in meters. + Defaults to 0 meters. + z_box (float, optional): How much empty space to set + above the path in meters. Defaults to 50 meters. """ if len(coord_list) < 2: @@ -70,6 +74,11 @@ class Path: ] self.z = z + self.size = ( + np.ceil(max(self.x_list)), + np.ceil(max(self.y_list)), + z + z_box + ) logger.debug("Path created.") diff --git a/src/pg_rad/physics/__init__.py b/src/pg_rad/physics/__init__.py new file mode 100644 index 0000000..fe9b62d --- /dev/null +++ b/src/pg_rad/physics/__init__.py @@ -0,0 +1,10 @@ +# do not expose internal logger when running mkinit +__ignore__ = ["logger"] +from pg_rad.physics import attenuation +from pg_rad.physics import fluence + +from pg_rad.physics.attenuation import (get_mass_attenuation_coeff,) +from pg_rad.physics.fluence import (phi_single_source,) + +__all__ = ['attenuation', 'fluence', 'get_mass_attenuation_coeff', + 'phi_single_source'] diff --git a/src/pg_rad/physics/attenuation.py b/src/pg_rad/physics/attenuation.py new file mode 100644 index 0000000..9af70e9 --- /dev/null +++ b/src/pg_rad/physics/attenuation.py @@ -0,0 +1,17 @@ +from importlib.resources import files + +from pandas import read_csv +from scipy.interpolate import interp1d + +from pg_rad.configs.filepaths import ATTENUATION_TABLE + + +def get_mass_attenuation_coeff( + *args + ) -> float: + csv = files('pg_rad.data').joinpath(ATTENUATION_TABLE) + data = read_csv(csv) + x = data["energy_mev"].to_numpy() + y = data["mu"].to_numpy() + f = interp1d(x, y) + return f(*args) diff --git a/src/pg_rad/physics/fluence.py b/src/pg_rad/physics/fluence.py new file mode 100644 index 0000000..8bae5f3 --- /dev/null +++ b/src/pg_rad/physics/fluence.py @@ -0,0 +1,37 @@ +import numpy as np + + +def phi_single_source( + r: float, + activity: float | int, + branching_ratio: float, + mu_mass_air: float, + air_density: float, + ) -> float: + """Compute the contribution of a single point source to the + primary photon fluence rate phi at position (x,y,z). + + Args: + r (float): [m] Distance to the point source. + activity (float | int): [Bq] Activity of the point source. + branching_ratio (float): Branching ratio for the photon energy E_gamma. + mu_mass_air (float): [cm^2/g] Mass attenuation coefficient for air. + air_density (float): [kg/m^3] Air density. + + Returns: + phi (float): [s^-1 m^-2] Primary photon fluence rate at distance r from + a point source. + """ + + # Linear photon attenuation coefficient in m^-1. + mu_mass_air *= 0.1 + mu_air = mu_mass_air * air_density + + phi = ( + activity + * branching_ratio + * np.exp(-mu_air * r) + / (4 * np.pi * r**2) + ) + + return phi diff --git a/src/pg_rad/plotting/__init__.py b/src/pg_rad/plotting/__init__.py new file mode 100644 index 0000000..271d2c9 --- /dev/null +++ b/src/pg_rad/plotting/__init__.py @@ -0,0 +1,7 @@ +# do not expose internal logger when running mkinit +__ignore__ = ["logger"] +from pg_rad.plotting import landscape_plotter + +from pg_rad.plotting.landscape_plotter import (LandscapeSlicePlotter,) + +__all__ = ['LandscapeSlicePlotter', 'landscape_plotter'] diff --git a/src/pg_rad/plotting/landscape_plotter.py b/src/pg_rad/plotting/landscape_plotter.py new file mode 100644 index 0000000..9c5552b --- /dev/null +++ b/src/pg_rad/plotting/landscape_plotter.py @@ -0,0 +1,75 @@ +import logging + +from matplotlib import pyplot as plt +from matplotlib.patches import Circle + +from pg_rad.landscape import Landscape + + +logger = logging.getLogger(__name__) + + +class LandscapeSlicePlotter: + def plot(self, landscape: Landscape, z: int = 0): + """Plot a top-down slice of the landscape at a height z. + + Args: + landscape (Landscape): the landscape to plot + z (int, optional): Height at which to plot slice. Defaults to 0. + """ """ + + """ + self.z = z + fig, ax = plt.subplots() + + self._draw_base(ax, landscape) + self._draw_path(ax, landscape) + self._draw_point_sources(ax, landscape) + + ax.set_aspect("equal") + plt.show() + + def _draw_base(self, ax, landscape): + width, height = landscape.size[:2] + ax.set_xlim(right=width) + ax.set_ylim(top=height) + ax.set_xlabel("X [m]") + ax.set_ylabel("Y [m]") + ax.set_title(f"Landscape (top-down, z = {self.z})") + + def _draw_path(self, ax, landscape): + if landscape.path.z < self.z: + ax.plot(landscape.path.x_list, landscape.path.y_list, 'bo-') + else: + logger.warning( + "Path is above the slice height z." + "It will not show on the plot." + ) + + def _draw_point_sources(self, ax, landscape): + for s in landscape.point_sources: + if s.z <= self.z: + dot = Circle( + (s.x, s.y), + radius=5, + color=s.color, + zorder=5 + ) + + ax.text( + s.x + 0.06, + s.y + 0.06, + s.name, + color=s.color, + fontsize=10, + ha="left", + va="bottom", + zorder=6 + ) + + ax.add_patch(dot) + else: + logger.warning( + f"Source {s.name} is above slice height z." + "It will not show on the plot." + ) diff --git a/tests/test_attenuation_functions.py b/tests/test_attenuation_functions.py new file mode 100644 index 0000000..385ace5 --- /dev/null +++ b/tests/test_attenuation_functions.py @@ -0,0 +1,29 @@ +import pytest + +from pg_rad.physics import get_mass_attenuation_coeff + + +@pytest.mark.parametrize("energy,mu", [ + (1.00000E-03, 3.606E+03), + (1.00000E-02, 5.120E+00), + (1.00000E-01, 1.541E-01), + (1.00000E+00, 6.358E-02), + (1.00000E+01, 2.045E-02) +]) +def test_exact_attenuation_retrieval(energy, mu): + """ + Test if retrieval of values that are exactly in the table is correct. + """ + func_mu = get_mass_attenuation_coeff(energy) + assert pytest.approx(func_mu, rel=1E-6) == mu + + +@pytest.mark.parametrize("energy,mu", [ + (0.662, 0.0778) +]) +def test_attenuation_interpolation(energy, mu): + """ + Test retreival for Cs-137 mass attenuation coefficient. + """ + interp_mu = get_mass_attenuation_coeff(energy) + assert pytest.approx(interp_mu, rel=1E-2) == mu diff --git a/tests/test_fluence_rate.py b/tests/test_fluence_rate.py new file mode 100644 index 0000000..c2615f6 --- /dev/null +++ b/tests/test_fluence_rate.py @@ -0,0 +1,34 @@ +from math import dist, exp, pi + +import pytest + +from pg_rad.landscape import LandscapeDirector + + +@pytest.fixture +def phi_ref(): + A = 100 # MBq + b = 0.851 + mu_mass_air = 0.0778 # cm^2/g + air_density = 1.243 # kg/m^3 + r = dist((0, 0, 0), (10, 10, 0)) # m + + A *= 1E9 # Convert to Bq + mu_mass_air *= 0.1 # Convert to m^2/kg + + mu_air = mu_mass_air * air_density # [m^2/kg] x [kg/m^3] = [m^-1] + + # [s^-1] x exp([m^-1] x [m]) / [m^-2] = [s^-1 m^-2] + phi = A * b * exp(-mu_air * r) / (4 * pi * r**2) + return phi + + +@pytest.fixture +def test_landscape(): + landscape = LandscapeDirector().build_test_landscape() + return landscape + + +def test_single_source_fluence(phi_ref, test_landscape): + phi = test_landscape.calculate_fluence_at((10, 10, 0)) + assert pytest.approx(phi, rel=1E-3) == phi_ref diff --git a/tests/test_sources.py b/tests/test_sources.py index 156276a..39d3063 100644 --- a/tests/test_sources.py +++ b/tests/test_sources.py @@ -1,27 +1,36 @@ import numpy as np import pytest -from pg_rad.sources import PointSource +from pg_rad.objects import PointSource +from pg_rad.isotopes import CS137 + @pytest.fixture def test_sources(): + iso = CS137() pos_a = np.random.rand(3) pos_b = np.random.rand(3) - a = PointSource(*tuple(pos_a), strength = None) - b = PointSource(*tuple(pos_b), strength = None) + a = PointSource(pos=pos_a, activity=None, isotope=iso) + b = PointSource(pos=pos_b, activity=None, isotope=iso) return pos_a, pos_b, a, b + def test_if_distances_equal(test_sources): - """Verify whether from PointSource A to PointSource B is the same as B to A.""" + """ + Verify whether distance from PointSource A to B is the same as B to A. + """ _, _, a, b = test_sources assert a.distance_to(b) == b.distance_to(a) + def test_distance_calculation(test_sources): - """Verify whether distance between two PointSources is calculated correctly.""" + """ + Verify whether distance between PointSources is calculated correctly. + """ pos_a, pos_b, a, b = test_sources @@ -32,4 +41,4 @@ def test_distance_calculation(test_sources): assert np.isclose( a.distance_to(b), np.sqrt(dx**2 + dy**2 + dz**2), - rtol = 1e-12) \ No newline at end of file + rtol=1e-12)