add fluence function to landscape and update to new position system

Improve PEP8 adherance

.github/workflows/ci-docs.yml

@@ -17,7 +17,7 @@ jobs:
           git config user.email 41898282+github-actions[bot]@users.noreply.github.com
       - uses: actions/setup-python@v5
         with:
-          python-version: 3.x
+          python-version: 3.12.9
       - run: echo "cache_id=$(date --utc '+%V')" >> $GITHUB_ENV
       - uses: actions/cache@v4
         with:

.gitignore

@@ -1,3 +1,6 @@
+# Custom
+dev-tools/
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[codz]

README.md

@@ -23,13 +23,27 @@ With Python verion `>=3.12.4` and `<3.13`, create a virtual environment and inst
 ```
 python3 -m venv .venv
 source .venv/bin/activate
-(venv) pip install -e .[dev]
+```
+
+With the virtual environment activated, run:
+
+```
+pip install -e .[dev]
 ```
 
 ## Tests
 
-Tests can be run with `pytest` from the root directory of the repository.
+Tests can be run with `pytest` from the root directory of the repository. With the virtual environment activated, run:
 
 ```
-(venv) pytest
-```
+pytest
+```
+
+## Local viewing of documentation
+
+PG-RAD uses [Material for MkDocs](https://squidfunk.github.io/mkdocs-material/) for generating documentation. It can be locally viewed by (in the venv) running:
+```
+mkdocs serve
+```
+
+where you can add the `--livereload` flag to automatically update the documentation as you write to the Markdown files.

docs/javascripts/mathjax.js

@@ -0,0 +1,18 @@
+window.MathJax = {
+  tex: {
+    inlineMath: [['$', '$'], ["\\(", "\\)"]],
+    displayMath: [['$$', '$$'], ["\\[", "\\]"]],
+    processEscapes: true,
+    processEnvironments: true
+  },
+  options: {
+    processHtmlClass: "arithmatex"
+  }
+};
+
+document$.subscribe(() => { 
+  MathJax.startup.output.clearCache()
+  MathJax.typesetClear()
+  MathJax.texReset()
+  MathJax.typesetPromise()
+})

docs/pg-rad-in-python.md

@@ -1,5 +0,0 @@
----
-title: Using PG-RAD as a module
----
-
-Consult the API documentation in the side bar.

mkdocs.yml

@@ -28,8 +28,16 @@ markdown_extensions:
   - pymdownx.inlinehilite
   - pymdownx.snippets
   - pymdownx.superfences
+  - pymdownx.arithmatex:
+      generic: true
+
+extra_javascript:
+  - javascripts/mathjax.js
+  - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js
 
 plugins:
+- mkdocs-jupyter:
+    execute: false
 - mkdocstrings:
     enabled: !ENV [ENABLE_MKDOCSTRINGS, true]
     default_handler: python

pyproject.toml

@@ -7,7 +7,7 @@ where = ["src"]
 
 [project]
 name = "pg-rad"
-version = "0.2.0"
+version = "0.2.1"
 authors = [
   { name="Pim Nelissen", email="pi0274ne-s@student.lu.se" },
 ]
@@ -29,4 +29,4 @@ Homepage = "https://github.com/pim-n/pg-rad"
 Issues = "https://github.com/pim-n/pg-rad/issues"
 
 [project.optional-dependencies]
-dev = ["pytest", "notebook", "mkdocs-material", "mkdocstrings-python"]
+dev = ["pytest", "mkinit", "notebook", "mkdocs-material", "mkdocstrings-python", "mkdocs-jupyter"]

src/pg_rad/data/__init__.py

@@ -0,0 +1 @@
+__all__ = []

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

src/pg_rad/dataloader/__init__.py

@@ -0,0 +1,8 @@
+# do not expose internal logger when running mkinit
+__ignore__ = ["logger"]
+
+from pg_rad.dataloader import dataloader
+
+from pg_rad.dataloader.dataloader import (load_data,)
+
+__all__ = ['dataloader', 'load_data']

src/pg_rad/dataloader/dataloader.py

@@ -1,12 +1,14 @@
+import logging
+
 import pandas as pd
 
-from pg_rad.logger import setup_logger
 from pg_rad.exceptions import DataLoadError, InvalidCSVError
 
-logger = setup_logger(__name__)
+logger = logging.getLogger(__name__)
+
 
 def load_data(filename: str) -> pd.DataFrame:
-    logger.debug(f"Attempting to load data from {filename}")
+    logger.debug(f"Attempting to load file: {filename}")
 
     try:
         df = pd.read_csv(filename, delimiter=',')
@@ -23,4 +25,5 @@ def load_data(filename: str) -> pd.DataFrame:
         logger.exception(f"Unexpected error while loading {filename}")
         raise DataLoadError("Unexpected error while loading data") from e
 
-    return df
+    logger.debug(f"File loaded: {filename}")
+    return df

src/pg_rad/exceptions/__init__.py

@@ -0,0 +1,10 @@
+# do not expose internal logger when running mkinit
+__ignore__ = ["logger"]
+
+from pg_rad.exceptions import exceptions
+
+from pg_rad.exceptions.exceptions import (ConvergenceError, DataLoadError,
+                                          InvalidCSVError,)
+
+__all__ = ['ConvergenceError', 'DataLoadError', 'InvalidCSVError',
+           'exceptions']

src/pg_rad/exceptions/exceptions.py

@@ -1,8 +1,10 @@
 class ConvergenceError(Exception):
     """Raised when an algorithm fails to converge."""
 
+
 class DataLoadError(Exception):
     """Base class for data loading errors."""
 
+
 class InvalidCSVError(DataLoadError):
-    """Raised when a file is not a valid CSV."""
+    """Raised when a file is not a valid CSV."""

src/pg_rad/isotopes/__init__.py

@@ -0,0 +1,10 @@
+# do not expose internal logger when running mkinit
+__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__ = ['CS137', 'Isotope', 'isotope', 'presets']

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.
 
@@ -5,19 +8,20 @@ class Isotope:
         name (str): Full name (e.g. Caesium-137).
         E (float): Energy of the primary gamma in keV.
         b (float): Branching ratio for the gamma at energy E.
-    """  
+    """
     def __init__(
             self,
             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.b = b
+        self.mu_mass_air = get_mass_attenuation_coeff(E / 1000)

src/pg_rad/isotopes/presets.py

@@ -0,0 +1,4 @@
+from .isotope import Isotope
+
+
+CS137 = Isotope("Cs-137", E=661.66, b=0.851)

src/pg_rad/landscape/__init__.py

@@ -0,0 +1,8 @@
+# do not expose internal logger when running mkinit
+__ignore__ = ["logger"]
+
+from pg_rad.landscape import landscape
+
+from pg_rad.landscape.landscape import (Landscape, create_landscape_from_path,)
+
+__all__ = ['Landscape', 'create_landscape_from_path', 'landscape']

src/pg_rad/landscape/landscape.py

@@ -1,38 +1,47 @@
+import logging
+
 from matplotlib import pyplot as plt
 from matplotlib.patches import Circle
 import numpy as np
+from numpy.typing import ArrayLike
 
 from pg_rad.path import Path
-from pg_rad.sources import PointSource
+from pg_rad.objects import PointSource
+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 in 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'.
-    """    
+        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'.
+    """
     def __init__(
             self,
             air_density: float = 1.243,
             size: int | tuple[int, int, int] = 500,
-            scale = 'meters',
-            ):    
+            scale: str = 'meters'
+            ):
         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 an integer or a tuple of 3 integers.")
+            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] = []
-    
-    def plot(self, z = 0):
+        logger.debug("Landscape initialized.")
+
+    def plot(self, z: float | int = 0):
         """Plot a slice of the world at a height `z`.
 
         Args:
@@ -40,7 +49,7 @@ class Landscape:
 
         Returns:
             fig, ax: Matplotlib figure objects.
-        """        
+        """
         x_lim, y_lim, _ = self.world.shape
 
         fig, ax = plt.subplots()
@@ -49,9 +58,9 @@ class Landscape:
         ax.set_xlabel(f"X [{self.scale}]")
         ax.set_ylabel(f"Y [{self.scale}]")
 
-        if not self.path == None:
+        if self.path is not None:
             ax.plot(self.path.x_list, self.path.y_list, 'bo-')
-        
+
         for s in self.sources:
             if np.isclose(s.z, z):
                 dot = Circle(
@@ -73,25 +82,23 @@ class Landscape:
                 )
 
                 ax.add_patch(dot)
-                
+
         return fig, ax
 
     def add_sources(self, *sources: PointSource):
         """Add one or more point sources to the world.
 
         Args:
-            *sources (pg_rad.sources.PointSource): One or more sources, passed as
-            Source1, Source2, ...
+            *sources (pg_rad.sources.PointSource): One or more sources,
+            passed as Source1, Source2, ...
         Raises:
-            ValueError: If the 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)
+            ValueError: If the source is outside the boundaries of the
+            landscape.
+        """
+        if not any(
+            (0 <= source.pos[0] <= self.world.shape[0] or
+             0 <= source.pos[1] <= self.world.shape[1] or
+             0 <= source.pos[2] <= self.world.shape[2])
             for source in sources
         ):
             raise ValueError("One or more sources are outside the landscape!")
@@ -102,25 +109,44 @@ class Landscape:
         """
         Set the path in the landscape.
         """
+        if not isinstance(path, Path):
+            raise TypeError("path must be of type Path.")
         self.path = path
-    
-def create_landscape_from_path(path: Path, max_z = 500):
+
+    def calculate_fluence_at(self, pos: tuple):
+        total_phi = 0.
+        for source in self.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):
+        if self.path is None:
+            raise ValueError("Path is not set!")
+
+
+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 500 meters.
+        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.
-    """    
+        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 = Landscape(size=(max_x, max_y, max_z))
     landscape.path = path
-    return landscape
+    return landscape

src/pg_rad/logging/__init__.py

@@ -0,0 +1,5 @@
+from pg_rad.logging import logger
+
+from pg_rad.logging.logger import (setup_logger,)
+
+__all__ = ['logger', 'setup_logger']

src/pg_rad/logging/logger.py

@@ -1,17 +1,21 @@
 import logging
-import logging.config
 import pathlib
 
 import yaml
 
-def setup_logger(name):
-    logger = logging.getLogger(name)
+
+def setup_logger(log_level: str = "WARNING"):
+    levels = ["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"]
+
+    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"
 
     with open(config_file) as f:
         config = yaml.safe_load(f)
-    
+
+    config["loggers"]["root"]["level"] = log_level
+
     logging.config.dictConfig(config)
-    return logger

src/pg_rad/objects.py

@@ -1,33 +0,0 @@
-import math
-from typing import Self
-
-class Object:
-    def __init__(
-            self,
-            x: float,
-            y: float,
-            z: float,
-            name: str = "Unnamed object",
-            color: str = 'grey'):
-        """
-        A generic object.
-
-        Args:
-            x (float): X coordinate.
-            y (float): Y coordinate.
-            z (float): Z coordinate.
-            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
-        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),
-        )

src/pg_rad/objects/__init__.py

@@ -0,0 +1,13 @@
+# 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',
+           'sources']

src/pg_rad/objects/objects.py

@@ -0,0 +1,42 @@
+from typing import Self
+
+import numpy as np
+
+
+class BaseObject:
+    def __init__(
+            self,
+            pos: tuple[float, float, float],
+            name: str = "Unnamed object",
+            color: str = 'grey'):
+        """
+        A generic object.
+
+        Args:
+            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".
+        """
+
+        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 | 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

src/pg_rad/objects/sources.py

@@ -1,13 +1,17 @@
-from pg_rad.objects import Object
-from pg_rad.isotope import Isotope
+import logging
 
-class PointSource(Object):
+from .objects import BaseObject
+from pg_rad.isotopes import Isotope
+
+logger = logging.getLogger(__name__)
+
+
+class PointSource(BaseObject):
     _id_counter = 1
+
     def __init__(
             self,
-            x: float,
-            y: float,
-            z: float,
+            pos: tuple,
             activity: int,
             isotope: Isotope,
             name: str | None = None,
@@ -15,17 +19,16 @@ class PointSource(Object):
         """A point source.
 
         Args:
-            x (float): X coordinate.
-            y (float): Y coordinate.
-            z (float): Z coordinate.
+            pos (tuple): a position vector of length 3 (x,y,z).
             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".
+            color (str, optional): Matplotlib compatible color string.
+            Defaults to "red".
         """
-        
+
         self.id = PointSource._id_counter
         PointSource._id_counter += 1
 
@@ -33,11 +36,17 @@ class PointSource(Object):
         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}")
 
     def __repr__(self):
-        return f"PointSource(name={self.name}, pos={(self.x, self.y, self.z)}, isotope={self.isotope.name}, A={self.activity} MBq)"
+        repr_str = (f"PointSource(name={self.name}, "
+                    + f"pos={(self.x, self.y, self.z)}, "
+                    + f"A={self.activity} MBq), "
+                    + f"isotope={self.isotope.name}.")
+
+        return repr_str

src/pg_rad/path/__init__.py

@@ -0,0 +1,9 @@
+# do not expose internal logger when running mkinit
+__ignore__ = ["logger"]
+
+from pg_rad.path import path
+
+from pg_rad.path.path import (Path, PathSegment, path_from_RT90,
+                              simplify_path,)
+
+__all__ = ['Path', 'PathSegment', 'path', 'path_from_RT90', 'simplify_path']

src/pg_rad/path/path.py

@@ -1,4 +1,5 @@
 from collections.abc import Sequence
+import logging
 import math
 
 from matplotlib import pyplot as plt
@@ -7,9 +8,9 @@ import pandas as pd
 import piecewise_regression
 
 from pg_rad.exceptions import ConvergenceError
-from pg_rad.logger import setup_logger
 
-logger = setup_logger(__name__)
+logger = logging.getLogger(__name__)
+
 
 class PathSegment:
     def __init__(self, a: tuple[float, float], b: tuple[float, float]):
@@ -18,18 +19,18 @@ class PathSegment:
         Args:
             a (tuple[float, float]): The starting point (x_a, y_a).
             b (tuple[float, float]): The final point (x_b, y_b).
-        """        
+        """
         self.a = a
         self.b = b
 
     def get_length(self) -> float:
         return math.dist(self.a, self.b)
-    
+
     length = property(get_length)
 
     def __str__(self) -> str:
         return str(f"({self.a}, {self.b})")
-    
+
     def __getitem__(self, index) -> float:
         if index == 0:
             return self.a
@@ -38,26 +39,30 @@ class PathSegment:
         else:
             raise IndexError
 
+
 class Path:
     def __init__(
             self,
             coord_list: Sequence[tuple[float, float]],
             z: float = 0,
-            path_simplify = False
+            path_simplify: bool = False
                  ):
         """Construct a path of sequences based on a list of coordinates.
 
         Args:
-            coord_list (Sequence[tuple[float, float]]): List of x,y coordinates.
+            coord_list (Sequence[tuple[float, float]]): List of x,y
+            coordinates.
             z (float, optional): Height of the path. Defaults to 0.
-            path_simplify (bool, optional): Whether to pg_rad.path.simplify_path(). Defaults to False.
-        """        
-        
+            path_simplify (bool, optional): Whether to
+            pg_rad.path.simplify_path(). Defaults to False.
+        """
+
         if len(coord_list) < 2:
-            raise ValueError("Must provide at least two coordinates as a list of tuples, e.g. [(x1, y1), (x2, y2)]")
+            raise ValueError("Must provide at least two coordinates as a \
+                of tuples, e.g. [(x1, y1), (x2, y2)]")
 
         x, y = tuple(zip(*coord_list))
-    
+
         if path_simplify:
             try:
                 x, y = simplify_path(list(x), list(y))
@@ -69,13 +74,19 @@ class Path:
 
         coord_list = list(zip(x, y))
 
-        self.segments = [PathSegment(i, ip1) for i, ip1 in zip(coord_list, coord_list[1:])]
+        self.segments = [
+            PathSegment(i, ip1)
+            for i, ip1 in
+            zip(coord_list, coord_list[1:])
+            ]
 
         self.z = z
 
+        logger.debug("Path created.")
+
     def get_length(self) -> float:
         return sum([s.length for s in self.segments])
-    
+
     length = property(get_length)
 
     def __getitem__(self, index) -> PathSegment:
@@ -83,40 +94,43 @@ class Path:
 
     def __str__(self) -> str:
         return str([str(s) for s in self.segments])
-    
+
     def plot(self, **kwargs):
         """
         Plot the path using matplotlib.
-        """  
+        """
         plt.plot(self.x_list, self.y_list, **kwargs)
 
+
 def simplify_path(
         x: Sequence[float],
         y: Sequence[float],
         keep_endpoints_equal: bool = False,
         n_breakpoints: int = 3
         ):
-    """From full resolution x and y arrays, return a piecewise linearly approximated/simplified pair of x and y arrays.
+    """From full resolution x and y arrays, return a piecewise linearly
+    approximated/simplified pair of x and y arrays.
 
     This function uses the `piecewise_regression` package. From a full set of
     coordinate pairs, the function fits linear sections, automatically finding
     the number of breakpoints and their positions.
 
-    On why the default value of n_breakpoints is 3, from the `piecewise_regression`
-    docs:
+    On why the default value of n_breakpoints is 3, from the
+    `piecewise_regression` docs:
     "If you do not have (or do not want to use) initial guesses for the number
     of breakpoints, you can set it to n_breakpoints=3, and the algorithm will
     randomly generate start_values. With a 50% chance, the bootstrap restarting
     algorithm will either use the best currently converged breakpoints or
-    randomly generate new start_values, escaping the local optima in two ways in
-    order to find better global optima."
+    randomly generate new start_values, escaping the local optima in two ways
+    in order to find better global optima."
 
     Args:
         x (Sequence[float]): Full list of x coordinates.
         y (Sequence[float]): Full list of y coordinates.
         keep_endpoints_equal (bool, optional): Whether or not to force start
-        and end to be exactly equal to the original. This will worsen the linear
-        approximation at the beginning and end of path. Defaults to False.
+        and end to be exactly equal to the original. This will worsen the
+        linear approximation at the beginning and end of path. Defaults to
+        False.
         n_breakpoints (int, optional): Number of breakpoints. Defaults to 3.
 
     Returns:
@@ -127,25 +141,27 @@ def simplify_path(
         ConvergenceError: If the fitting algorithm failed to simplify the path.
 
     Reference:
-        Pilgrim, C., (2021). piecewise-regression (aka segmented regression) in Python. Journal of Open Source Software, 6(68), 3859, https://doi.org/10.21105/joss.03859.
+        Pilgrim, C., (2021). piecewise-regression (aka segmented regression)
+        in Python. Journal of Open Source Software,
+        6(68), 3859, https://doi.org/10.21105/joss.03859.
     """
-    
-    logger.debug(f"Attempting piecewise regression on path.")
+
+    logger.debug("Attempting piecewise regression on path.")
 
     pw_fit = piecewise_regression.Fit(x, y, n_breakpoints=n_breakpoints)
     pw_res = pw_fit.get_results()
 
-    if pw_res == None:
-        logger.error("Piecewise regression failed to converge.")
-        raise ConvergenceError("Piecewise regression failed to converge.") 
-    
+    if pw_res is None:
+        logger.warning("Piecewise regression failed to converge.")
+        raise ConvergenceError("Piecewise regression failed to converge.")
+
     est = pw_res['estimates']
 
     # extract and sort breakpoints
     breakpoints_x = sorted(
         v['estimate'] for k, v in est.items() if k.startswith('breakpoint')
     )
- 
+
     x_points = [x[0]] + breakpoints_x + [x[-1]]
 
     y_points = pw_fit.predict(x_points)
@@ -156,32 +172,38 @@ def simplify_path(
         y_points[-1] = y[-1]
 
     logger.info(
-        f"Piecewise regression reduced path from {len(x)-1} to {len(x_points)-1} segments."
+        f"Piecewise regression reduced path from \
+            {len(x)-1} to {len(x_points)-1} segments."
     )
 
     return x_points, y_points
 
+
 def path_from_RT90(
         df: pd.DataFrame,
         east_col: str = "East",
         north_col: str = "North",
         **kwargs
         ) -> Path:
-    """Construct a path from East and North formatted coordinates (RT90) in a Pandas DataFrame.
+    """Construct a path from East and North formatted coordinates (RT90)
+    in a Pandas DataFrame.
 
     Args:
-        df (pandas.DataFrame): DataFrame containing at least the two columns noted in the cols argument.
+        df (pandas.DataFrame): DataFrame containing at least the two columns
+        noted in the cols argument.
         east_col (str): The column name for the East coordinates.
         north_col (str): The column name for the North coordinates.
 
     Returns:
-        Path: A Path object built from the aquisition coordinates in the DataFrame.
-    """    
-    
+        Path: A Path object built from the aquisition coordinates in the
+        DataFrame.
+    """
+
     east_arr = np.array(df[east_col]) - min(df[east_col])
     north_arr = np.array(df[north_col]) - min(df[north_col])
 
     coord_pairs = list(zip(east_arr, north_arr))
 
     path = Path(coord_pairs, **kwargs)
-    return path
+    logger.debug("Loaded path from provided RT90 coordinates.")
+    return path

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']

src/pg_rad/physics/attenuation.py

@@ -0,0 +1,15 @@
+from importlib.resources import files
+
+from pandas import read_csv
+from scipy.interpolate import interp1d
+
+
+def get_mass_attenuation_coeff(
+        *args
+        ) -> float:
+    csv = files('pg_rad.data').joinpath('attenuation_table.csv')
+    data = read_csv(csv)
+    x = data["energy_mev"].to_numpy()
+    y = data["mu"].to_numpy()
+    f = interp1d(x, y)
+    return f(*args)

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

tests/test_attenuation_functions.py

@@ -0,0 +1,30 @@
+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),
+    (1.25, 0.06)
+])
+def test_attenuation_interpolation(energy, mu):
+    """
+    Test Cs-137 and Co-60 mass attenuation coefficients.
+    """
+    interp_mu = get_mass_attenuation_coeff(energy)
+    assert pytest.approx(interp_mu, rel=1E-2) == mu

tests/test_fluence_rate.py

@@ -0,0 +1,41 @@
+from math import dist, exp, pi
+
+import pytest
+
+from pg_rad.isotopes import CS137
+from pg_rad.landscape import Landscape
+from pg_rad.objects import PointSource
+
+
+@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 = Landscape()
+    source = PointSource(
+        pos=(0, 0, 0),
+        activity=100E9,
+        isotope=CS137)
+    landscape.add_sources(source)
+    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

tests/test_sources.py

@@ -1,28 +1,36 @@
 import numpy as np
 import pytest
 
-from pg_rad.objects import Source
+from pg_rad.objects import PointSource
+from pg_rad.isotopes import Isotope
+
 
 @pytest.fixture
 def test_sources():
+    iso = Isotope("test", E=662, b=0)
     pos_a = np.random.rand(3)
     pos_b = np.random.rand(3)
 
-    a = Source(*tuple(pos_a), strength = None)
-    b = Source(*tuple(pos_b), strength = None)
+    a = PointSource(pos_a, activity=None, isotope=iso)
+    b = PointSource(pos_b, activity=None, isotope=iso)
 
     return pos_a, pos_b, a, b
 
+
 def test_if_distances_equal(test_sources):
-    """_Verify whether from object A to object 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 static objects (e.g. sources) 
-    is calculated correctly._"""
+    """
+    Verify whether distance between PointSources is calculated correctly.
+    """
 
     pos_a, pos_b, a, b = test_sources
 
@@ -33,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)
+        rtol=1e-12)