Add background counts for 3x3 inch NaI. Add auto ROI lookup function from the FWHM. Add ROI test.

Feature detectors roadgen

.github/workflows/ci-tests.yml

@@ -0,0 +1,33 @@
+name: Tests
+on:
+  pull_request:
+    branches: ["main", "dev"]
+  workflow_dispatch:
+
+jobs:
+  tests:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ ubuntu-latest, windows-latest ]
+        python-version: [ '3.12' ]
+
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+          cache: pip
+      - name: Install Dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install -e .[dev]
+      - name: Run linting
+        run: |
+          flake8 ./src/
+
+      - name: Run tests
+        run: |
+          pytest     

.gitignore

@@ -176,7 +176,7 @@ cython_debug/
 #  be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
 #  and can be added to the global gitignore or merged into this file.  For a more nuclear
 #  option (not recommended) you can uncomment the following to ignore the entire idea folder.
-#.idea/
+.idea/
 
 # Abstra
 # Abstra is an AI-powered process automation framework.

README.md

@@ -1,3 +1,5 @@
+[![Python 3.12](https://img.shields.io/badge/python-3.12-blue.svg)](https://www.python.org/downloads/release/python-312/)
+[![Tests](https://github.com/pim-n/pg-rad/actions/workflows/ci-tests.yml/badge.svg)](https://github.com/pim-n/pg-rad/actions/workflows/ci-tests.yml)
 # pg-rad
 Primary Gamma RADiation landscape - Development
 
@@ -31,6 +33,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:

docs/API/landscape/create_landscape_from_path.md

@@ -1,4 +0,0 @@
----
-title: pg_rad.landscape.create_landscape_from_path
----
-::: pg_rad.landscape.create_landscape_from_path

docs/API/landscape/landscape.md

@@ -1,4 +0,0 @@
----
-title: pg_rad.landscape.Landscape
----
-::: pg_rad.landscape.Landscape

docs/API/objects/object.md

@@ -1,5 +0,0 @@
----
-title: pg_rad.objects.Object
----
-
-::: pg_rad.objects.Object

docs/API/path/path.md

@@ -1,4 +0,0 @@
----
-title: pg_rad.path.Path
----
-::: pg_rad.path.Path

docs/API/path/path_from_RT90.md

@@ -1,5 +0,0 @@
----
-title: pg_rad.path.path_from_RT90
----
-::: pg_rad.path.path_from_RT90
-

docs/API/path/simplify_path.md

@@ -1,4 +0,0 @@
----
-title: pg_rad.path.simplify_path
----
-::: pg_rad.path.simplify_path

docs/API/sources/point_source.md

@@ -1,5 +0,0 @@
----
-title: pg_rad.sources.PointSource
----
-
-::: pg_rad.sources.PointSource

docs/config-spec.md

@@ -0,0 +1,218 @@
+!!! note
+    To get started quickly, you may copy and modify one of the example configs found [here](quickstart.md#example-configs).
+
+
+The config file must be a [YAML](https://yaml.org/) file. YAML is a serialization language that works with key-value pairs, but in a syntax more readable than some other alternatives. In YAML, the indentation matters. I
+
+
+The remainder of this chapter will explain the different required and optionals keys, what they represent, and allowed values. 
+
+## Required keys
+
+### Simulation options
+
+The first step is to name the simulation, and define the speed of the vehicle (assumed constant) and acquisition time.
+
+#### Landscape name
+
+The name is a string, which may include spaces, numbers and special characters.
+
+Examples:
+
+```yaml
+name: test_landscape
+```
+```yaml
+name: Test Landscape 1
+```
+
+#### Acquisition time
+
+The acquisition time of the detector in seconds.
+
+Example:
+
+```yaml
+acquisition_time: 1
+```
+
+!!! note
+    All units in the config file must be specified in SI units, e.g. meters and seconds, unless the key contains a unit itself (e.g. `activity_MBq` means activity in MegaBequerels).
+
+#### Vehicle speed
+
+The speed of the vehicle in m/s. Currently, the vehicle speed must be assumed constant. An example could be
+
+```yaml
+speed: 13.89 # this is approximately 50 km/h
+```
+
+!!! note
+    The text after the `#` signifies a comment. PG-RAD will ignore this, but it can be helpful for yourself to write notes.
+
+
+### Path
+
+The `path` keyword is used to create a path for the detector to travel along. There are two ways to specify a path; from experimental data or by specifying a procedural path.
+
+#### Path - Experimental data
+
+Currently the only supported coordinate format is the RT90 (East, North) coordinate system. If you have experimental data in CSV format with columns for these coordinates, then you can load that path into PG-RAD as follows:
+
+```yaml
+path:
+  file: path/to/experimental_data.csv
+  east_col_name: East
+  north_col_name: North
+```
+
+#### Path - Procedural path
+
+Alternatively, you can let PG-RAD generate a path for you. A procedural path can be specified with at least two subkeys: `length` and `segments`. 
+
+Currently supported segments are: `straight`, `turn_left` and `turn_right`, and are provided in a list under the `segments` subkey as follows:
+
+```yaml
+path:
+  segments:
+    - straight
+    - turn_left
+    - straight
+```
+
+The length must also be specified, using the `length` subkey. `length` can be specified in two ways: a list with the same length as the `segments` list
+
+```yaml
+path:
+  segments:
+    - straight
+    - turn_left
+    - straight
+  length:
+    - 500
+    - 250
+    - 500
+```
+
+which will assign that length (meters) to each segment. Alternatively, a single number can be passed:
+
+```yaml
+path:
+  segments:
+    - straight
+    - turn_left
+    - straight
+  length: 1250
+```
+
+Setting the length for the total path will cause PG-RAD to *randomly assign* portions of the total length to each segment.
+
+Finally, there is also an option to specify the turn angle in degrees:
+
+```yaml
+path:
+  segments:
+    - straight
+    - turn_left: 90
+    - straight
+  length: 1250
+```
+
+Like with the lengths, if a turn segment has no angle specified, a random one (within pre-defined limits) will be taken.
+
+!!! warning
+    Letting PG-RAD randomly assign lengths and angles can cause (expected) issues. That is because of physics restrictions. If the combination of length, angle (radius) and velocity of the vehicle is such that the centrifugal force makes it impossible to take this turn, PG-RAD will raise an error. To fix it, you can 1) reduce the speed; 2) define a smaller angle for the turn; or 3) assign more length to the turn segment.
+
+!!! info
+    For more information about how procedural roads are generated, including the random sampling of lengths and angles, see X
+
+### Sources
+
+Currently, the only type of source supported is a point source. Point sources can be added under the `sources` key, where the **subkey is the name** of the source:
+
+```yaml
+sources:
+  my_source: ...
+```
+
+the source name should not contain spaces or special characters other than `_` or `-`. There are three required subkeys under `sources.my_source`, which are: `activity_MBq`, `isotope` and `position`.
+
+#### Source activity
+
+The source activity is in MegaBequerels and must be a strictly positive number:
+
+```yaml
+sources:
+  my_source:
+    activity_MBq: 100
+```
+
+#### Source isotope
+
+The isotope for the point source. This must be a string, following the naming convention of the symbol followed by the number of nucleons, e.g. `Cs137`:
+
+```yaml
+sources:
+  my_source:
+    activity_MBq: 100
+    isotope: Cs137
+```
+
+!!! info
+    Currently the following isotopes are supported: `Cs137`
+
+#### Source position
+
+There are two ways to specify the source position. Either with absolute (x,y,z) coordinates
+
+```yaml
+sources:
+  my_source:
+    activity_MBq: 100
+    isotope: Cs137
+    position: [0, 0, 0]
+```
+
+or relative to the path, using the subkeys `along_path`, `dist_from_path` and `side`
+
+```yaml
+sources:
+  my_source:
+    activity_MBq: 100
+    isotope: Cs137
+    position:
+      along_path: 100
+      dist_from_path: 50
+      side: left
+```
+
+Note that side is relative to the direction of travel. The path will by default start at (x,y) = (0,0) and initial heading is parallel to the x-axis.
+
+### Detector
+
+The final required key is the `detector`. Currently, only isotropic detectors are supported. Nonetheless, you must specify it with `name`, `is_isotropic` and `efficiency`:
+
+```yaml
+detector:
+  name: test
+  is_isotropic: True
+  efficiency: 0.02
+```
+
+Note there are some existing detectors available, where efficiency is not required and will be looked up by PG-RAD itself:
+
+```yaml
+detector:
+  name: NaIR
+  is_isotropic: True
+```
+
+## Optional keys
+
+The following subkeys are optional and should be put under the `options` key.
+
+```yaml
+options:
+  air_density_kg_per_m3: 1.243
+  seed: 1234
+```

docs/explainers/prefab_roads.ipynb

@@ -0,0 +1,118 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "1a063d05",
+   "metadata": {},
+   "source": [
+    "# Pseudo-ramdom procedural roads\n",
+    "\n",
+    "Suppose one wishes to describe a road between A and B in terms of segments\n",
+    "\n",
+    "$$\n",
+    "\\text{straight, turn left, straight, turn right, straight}\n",
+    "$$\n",
+    "\n",
+    "Let's see how we can create a random road of length $L$ from a pre-determined set of prefabs.\n",
+    "\n",
+    "#### Random apportionment of total length\n",
+    "\n",
+    "Suppose we want to build a road of length $L$ out of $K$ segments. The total number of waypoints $N$ depends on the step size $\\Delta s$:\n",
+    "\n",
+    "$$\n",
+    "N = \\frac{L}{\\Delta s}.\n",
+    "$$\n",
+    "\n",
+    "Let $\\left( p_1, p_2, \\dots, p_K \\right)$ represent the proportion of $N$ that each prefab will be assigned, where $\\sum p_i = 1$. One useful distribution here is the [Dirichlet distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution), which is parametrized by a vector $\\mathbf{\\alpha} = \\left(\\alpha_1, \\alpha_2, \\dots, \\alpha_K \\right)$. The special case where all $\\alpha_i$, the scalar parameter $\\alpha$ is called a *concentration parameter*. Setting the same $\\alpha$ across the entire parameter space makes the distribution symmetric, meaning no prior assumptions are made regarding the proportion of $N$ that will be assigned to each segment. $\\alpha = 1$ leads to what is known as a flat Dirichlet distribution, whereas higher values lead to more dense and evenly distributed $\\left( p_1, p_2, \\dots, p_K \\right)$. On the other hand, keeping $\\alpha \\leq 1$ gives a sparser distribution which can lead to larger variance in apportioned number of waypoints to $\\left( p_1, p_2, \\dots, p_K \\right)$.\n",
+    "\n",
+    "#### Expectation value and variance of Dirichlet distribution\n",
+    "\n",
+    "Suppose we draw our samples for proportion of length from the Dirichlet distribution\n",
+    "\n",
+    "$$\n",
+    "(p_1, p_2, \\ldots, p_K) \\sim \\text{Dirichlet}(\\alpha, \\alpha, \\ldots, \\alpha)\n",
+    "$$\n",
+    "\n",
+    "with $\\alpha _{0}=\\sum _{i=1}^{K}\\alpha _{i}$, the mean and variance are then\n",
+    "\n",
+    "$$\n",
+    "\\operatorname {E} [p_{i}]={\\frac {\\alpha _{i}}{\\alpha _{0}}}, \\; \\operatorname {Var} [p_{i}]={\\frac {\\alpha _{i}(\\alpha _{0}-\\alpha _{i})}{\\alpha _{0}^{2}(\\alpha _{0}+1)}}.\n",
+    "$$\n",
+    "\n",
+    "If $\\alpha$ is a scalar, then $\\alpha _{0}= K \\alpha$ and the above simplifies to\n",
+    "\n",
+    "$$\n",
+    "\\operatorname {E} [p_{i}]={\\frac {\\alpha}{K \\alpha}}={\\frac {1}{K}}, \\; \\operatorname {Var} [p_{i}]={\\frac {\\alpha(K \\alpha -\\alpha)}{(K \\alpha)^{2}(K \\alpha +1)}}.\n",
+    "$$\n",
+    "\n",
+    "We see that $\\operatorname {Var} [p_{i}] \\propto \\frac{1}{\\alpha}$ meaning that the variance reduces with increasing $\\alpha$. We can simply scale the proportions\n",
+    "\n",
+    "$$\n",
+    "(N \\cdot p_1, N \\cdot p_2, \\ldots, N \\cdot p_K)\n",
+    "$$\n",
+    "\n",
+    "to get the randomly assigned number of waypoints for each prefab. We now have a distribution which can give randomly assigned lengths to a given list of prefabs, with a parameter to control the degree of randomness. With a large concentration parameter $\\alpha$, the distribution of lengths will be more uniform, with each prefab getting $N \\cdot \\operatorname {E} [p_{i}]={\\frac {N}{K}}$ waypoints assigned to it. Likewise, keeping $\\alpha$ low increases variance and allows for a more random assignment of proportions of waypoints to each prefab segment.\n",
+    "\n",
+    "#### Random angles\n",
+    "\n",
+    "Suppose a turn of a pre-defined arc length $l$ made of $N/K$ waypoints. If one wants to create a random angle, one has to keep in mind that the minimum radius $R_{min}$ depends on the speed of the vehicle and the weather conditions:\n",
+    "\n",
+    "$$\n",
+    "R_{\\text{min,vehicle}} = \\frac{v^2}{g\\mu},\n",
+    "$$\n",
+    "\n",
+    "where\n",
+    "- $v$ is the velocity of the vehicle in $\\text{m/s}$,\n",
+    "- $g$ is the gravitational acceleration (about $9.8$ $\\text{m/s}^{2}$), and\n",
+    "- $\\mu$ is the friction coefficient (about $0.7$ for dry asphalt).\n",
+    "\n",
+    "A regular turn (not a U-turn or roundabout) should also have an lower and upper limit on the angle, say, 30 degrees to 90 degrees for a conservative estimate. In terms of radii, it becomes\n",
+    "\n",
+    "$$\n",
+    "R_{\\text{min}} = \\max\\left(R_{\\text{min,vehicle}}, \\frac{l}{\\pi/2}\\right)\n",
+    "$$\n",
+    "\n",
+    "and\n",
+    "\n",
+    "$$\n",
+    "R_{\\text{max}} = \\frac{l}{\\pi/6}.\n",
+    "$$\n",
+    "\n",
+    "We then sample\n",
+    "\n",
+    "$$\n",
+    "R \\sim \\text{Uniform}\\left(R_{\\text{min}}, R_{\\text{max\\_angle}}\\right)\n",
+    "$$\n",
+    "\n",
+    "and obtain a random radius for a turn of arc length $l$ with limits to ensure the radius is large enough given the velocity of the vehicle. Finally, the curvature profile is related to the radius by\n",
+    "\n",
+    "$$\n",
+    "\\kappa = \\frac{1}{R}\n",
+    "$$\n",
+    "\n",
+    "which means that the curvature profile of a turn is simply a vector $\\mathbf{\\kappa} = (1/R, \\dots, 1/R)$ with a length of $N/K$ waypoints."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": ".venv",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/index.md

@@ -2,29 +2,17 @@
 
 Primary Gamma RADiation Landscapes (PG-RAD) is a Python package for research in source localization. It can simulate mobile gamma spectrometry data acquired from vehicle-borne detectors along a predefined path (e.g. a road).
 
-## Requirements
+## About
 
-PG-RAD requires Python `3.12`. The guides below assume a unix-like system.
+This software has been developed as part of dissertation work for the degree of master of Computational Science and Physics at Lund University, Sweden. The work has been done at the department of Medical Radiation Physics (MSF), Faculty of Medicine. The radiological emergency preparedness research group of MSF is assigned by the Swedish Radiation Safety Authority (SSM) to aid in preparation for effective mitigation of radiological or nuclear disasters on Swedish soil.
 
-## Installation (CLI)
+## Value proposition
 
-<!--pipx seems like a possible option to install python package in a contained environment on unix-->
+PG-RAD is a toolbox that allows for simulation of detector response for a wide variety of source localization scenarios. The strength of the software lies in its simple and minimal configuration and user input, while its flexibility allows for reconstruction of specific scenarios with relative ease. PG-RAD is also general enough that novel methods such as UAV-borne detectors can be simulated and evaluated.
 
-Lorem ipsum
+User input takes the form of an input file (YAML), describing the path, detector and source(s), and optional parameters. The output of the program is visualizations of the world (the path and sources), as well as the detector count rate as a function of distance travelled along the path.
 
-## Installation (Python module)
-
-If you are interested in using PG-RAD in another Python project, create a virtual environment first:
-
-```
-python3 -m venv .venv
-```
-
-Then install PG-RAD in it:
-
-```
-source .venv/bin/activate
-(.venv) pip install git+https://github.com/pim-n/pg-rad
+Users can provide experimental / geographical coordinates representing real roads. Alternatively, users can let PG-RAD generate a procedural road, where the user can easily control what that road should look like. The user can specify a single point source, several point sources, as well as a field of radioactive material covering a large area.
 ```
 
 See how to get started with PG-RAD with your own Python code [here](pg-rad-in-python).

docs/installation.md

@@ -0,0 +1,47 @@
+## Requirements
+
+PG-RAD requires Python `>=3.12.4` and `<3.13`. It has been tested on `3.12.9`. The guides below assume a unix-like system. You can check the Python version you have installed as follows:
+
+```
+python --version
+```
+
+If you don't have the right version installed there are various ways to get a compatible version, such as [pyenv](https://github.com/pyenv/pyenv?tab=readme-ov-file#installation).
+
+## Installation (CLI)
+
+<!--pipx seems like a possible option to install python package in a contained environment on unix-->
+
+Lorem ipsum
+
+## Installation (Python module)
+
+If you are interested in using PG-RAD in another Python project, create a virtual environment first:
+
+```
+python -m venv .venv
+```
+
+Then install PG-RAD in it:
+
+```
+source .venv/bin/activate
+(.venv) pip install git+https://github.com/pim-n/pg-rad
+```
+
+See how to get started with PG-RAD with your own Python code [here](pg-rad-in-python).
+
+## For developers
+```
+git clone https://github.com/pim-n/pg-rad
+cd pg-rad
+git checkout dev
+```
+
+or
+
+```
+git@github.com:pim-n/pg-rad.git
+cd pg-rad
+git checkout dev
+```

docs/pg-rad-in-cli.md

@@ -1,4 +0,0 @@
----
-title: Using PG-RAD in CLI
----
-Lorem ipsum.

docs/quickstart.md

@@ -0,0 +1,187 @@
+## Installation
+
+See the [installation guide](installation.md).
+
+## Test your installation
+
+First, see if PG-RAD is available on your system by typing
+
+```
+pgrad --help
+```
+
+You should get output along the lines of
+
+```
+usage: pg-rad [-h] ...
+
+Primary Gamma RADiation landscape tool
+
+...
+```
+
+If you get something like `pgrad: command not found`, please consult the [installation guide](installation.md).
+
+You can run a quick test scenario as follows:
+
+```
+pgrad --test
+```
+
+This should produce a plot of a scenario containing a single point source and a path.
+
+## Running PG-RAD
+
+In order to use the CLI for your own simulations, you need to provide a *config file*. To run with your config, run
+
+```
+pgrad --config path/to/my_config.yml
+```
+
+where `path/to/my_config.yml` points to your config file.
+
+## Example configs
+
+The easiest way is to take one of these example configs, and adjust them as needed. Alternatively, there is a detailed guide on how to write your own config file [here](config-spec.md).
+
+=== "Example 1"
+
+    The position can be defined relative to the path. `along_path` means at what distance traveled along the path the source is found. If the path is 200 meters long and `along_path` is `100` then the source is halfway along the path. `dist_from_path` is the distance in meters from the path. `side` is the side of the path the source is located. This is relative to the direction the path is traveled.
+
+    ``` yaml
+    name: Example 1
+    speed: 13.89
+    acquisition_time: 1
+
+    path:
+      file: path/to/exp_coords.csv
+      east_col_name: East
+      north_col_name: North
+
+    sources:
+      source1:
+      activity_MBq: 1000
+      isotope: CS137
+      position:
+        along_path: 100
+        dist_from_path: 50
+        side: left
+    
+    detector:
+      name: dummy
+      is_isotropic: True
+    ```
+
+=== "Example 2"
+
+    The position can also just be defined with (x,y,z) coordinates.
+
+    ``` yaml
+    name: Example 2
+    speed: 13.89
+    acquisition_time: 1
+
+    path:
+      file: path/to/exp_coords.csv
+      east_col_name: East
+      north_col_name: North
+
+    sources:
+      source1:
+        activity_MBq: 1000
+        isotope: CS137
+        position: [104.3, 32.5, 0]
+      source2:
+        activity_MBq: 100
+        isotope: CS137
+        position: [0, 0, 0]
+        
+    detector:
+      name: dummy
+      is_isotropic: True
+    ```
+
+=== "Example 3"
+
+    This is an example of a procedural path with random apportionment of total length and random angles being assigned to turns. The parameter `alpha` is optional, and is related to randomness. A higher value leads to more uniform apportionment of lengths and a lower value to more random apportionment. More information about `alpha` can be found [here](pg-rad-config-spec.md).
+
+    ``` yaml
+    name: Example 3
+    speed: 8.33
+    acquisition_time: 1
+
+    path:
+      length: 1000
+      segments:
+        - straight
+        - turn_left
+        - straight
+      alpha: 100
+
+    sources:
+      source1:
+        activity_MBq: 1000
+        isotope: CS137
+        position: [0, 0, 0]
+        
+    detector:
+      name: dummy
+      is_isotropic: True
+    ```
+
+=== "Example 4"
+
+    This is an example of a procedural path that is partially specified. Note that turn_left now is a key for the corresponding angle of 45 degrees. The length is still divided randomly
+
+    ``` yaml
+    name: Example 4
+    speed: 8.33
+    acquisition_time: 1
+
+    path:
+      length: 1000
+      segments:
+        - straight
+        - turn_left: 45
+        - straight
+
+    sources:
+      source1:
+        activity_MBq: 1000
+        isotope: CS137
+        position: [0, 0, 0]
+        
+    detector:
+      name: dummy
+      is_isotropic: True
+    ```
+
+=== "Example 5"
+
+    This is an example of a procedural path that is fully specified. See how length is now a list matching the length of the segments.
+
+    ``` yaml
+    name: Example 5
+    speed: 8.33
+    acquisition_time: 1
+
+    path:
+      length:
+        - 400
+        - 200
+        - 400
+      segments:
+        - straight
+        - turn_left: 45
+        - straight
+
+    sources:
+      source1:
+        activity_MBq: 1000
+        isotope: CS137
+        position: [0, 0, 0]
+        
+    detector:
+      name: dummy
+      is_isotropic: True
+    ```

mkdocs.yml

@@ -30,6 +30,11 @@ markdown_extensions:
   - pymdownx.superfences
   - pymdownx.arithmatex:
       generic: true
+  - admonition
+  - pymdownx.details
+  - pymdownx.tabbed:
+      alternate_style: true
+      combine_header_slug: true
 
 extra_javascript:
   - javascripts/mathjax.js
@@ -47,3 +52,12 @@ plugins:
         options:
           show_source: false
           show_root_heading: false
+
+nav:
+  - Home: index.md
+  - Installation Guide: installation.md
+  - Quickstart Guide: quickstart.md
+  - 'Tutorial: Writing a Config File': config-spec.md
+  - Explainers:
+    - explainers/planar_curve.ipynb
+    - explainers/prefab_roads.ipynb

pyproject.toml

@@ -5,6 +5,10 @@ build-backend = "setuptools.build_meta"
 [tool.setuptools.packages.find]
 where = ["src"]
 
+[tool.setuptools.package-data]
+"pg_rad.data" = ["*.csv"]
+"pg_rad.configs" = ["*.yml"]
+
 [project]
 name = "pg-rad"
 version = "0.2.1"
@@ -18,15 +22,19 @@ dependencies = [
     "matplotlib>=3.9.2",
     "numpy>=2",
     "pandas>=2.3.1",
-    "piecewise_regression==1.5.0",
-    "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"]
+
+dev = ["pytest", "mkinit", "notebook", "mkdocs-material", "mkdocstrings-python", "mkdocs-jupyter", "flake8"]

src/pg_rad/background/background.py

@@ -0,0 +1,51 @@
+from importlib.resources import files
+from typing import Tuple
+import numpy as np
+from pandas import read_csv
+
+from pg_rad.configs.defaults import FWHM_PARAMS
+from pg_rad.detector.detector import Detector
+
+
+def generate_background(
+    cps_array: np.ndarray,
+) -> np.ndarray:
+    """
+    Generate synthetic background cps for a given detector and energy.
+    """
+    pass
+
+
+def fwhm(A: float, B: float, C: float, E: float) -> float:
+    return np.sqrt(A + B * E + C * E**2)
+
+
+def get_roi_from_fwhm(
+    detector: Detector,
+    energy_keV: float
+) -> Tuple[float, float]:
+    """
+    Get the region of interest for given primary gamma energy, using
+    the 3*FWHM rule.
+    """
+    A, B, C = FWHM_PARAMS.get(detector.type)
+    delta = 3*fwhm(A, B, C, energy_keV)
+    return (energy_keV-delta, energy_keV+delta)
+
+
+def get_cps_from_roi(
+    detector: Detector,
+    roi_lo: float,
+    roi_hi: float
+) -> float:
+
+    csv = files('pg_rad.data.backgrounds').joinpath(detector.name+'.csv')
+    data = read_csv(csv)
+
+    # get indices of nearest bins
+    idx_min = (data["Energy"] - roi_lo).abs().idxmin()
+    idx_max = (data["Energy"] - roi_hi).abs().idxmin()
+    idx_start, idx_end = sorted([idx_min, idx_max])
+
+    cps_sum = data.loc[idx_start:idx_end, "cps"].sum()
+    return cps_sum

src/pg_rad/configs/defaults.py

@@ -0,0 +1,32 @@
+# --- Physics defaults ---
+DEFAULT_AIR_DENSITY = 1.243  # kg/m^3, Skåne
+
+# --- Simulation defaults ---
+DEFAULT_SEED = None
+DEFAULT_ACQUISITION_TIME = 1.0
+
+# --- Source defaults ---
+DEFAULT_PATH_HEIGHT = 0.0
+DEFAULT_SOURCE_HEIGHT = 0.0
+
+# --- Segmented road defaults ---
+DEFAULT_MIN_TURN_ANGLE = 30.
+DEFAULT_MAX_TURN_ANGLE = 90.
+
+DEFAULT_FRICTION_COEFF = 0.7      # dry asphalt
+DEFAULT_GRAVITATIONAL_ACC = 9.81  # m/s^2
+DEFAULT_ALPHA = 100.
+
+# --- Detector efficiencies ---
+DETECTOR_EFFICIENCIES = {
+    "dummy": 1.0,
+    "NaIR": 0.0216,
+    "NaIF": 0.0254
+}
+
+# A, B, C parameters for different detector types
+# for formula FWHM(E) = sqrt(A + B*E + C*E^2)
+FWHM_PARAMS = {
+    "NaI": (-70.55, 2.678, 0.000602),
+    "HPGe": (1.778, 0.001843, 0.000001)
+}

src/pg_rad/configs/filepaths.py

@@ -0,0 +1,4 @@
+ATTENUATION_TABLE = 'attenuation_table.csv'
+ISOTOPE_TABLE = 'isotopes.csv'
+TEST_EXP_DATA = 'test_path_coords.csv'
+LOGGING_CONFIG = 'logging.yml'

src/pg_rad/configs/logging.yml

@@ -3,10 +3,13 @@ disable_existing_loggers: false
 formatters:
   simple:
     format: '%(asctime)s - %(levelname)s: %(message)s'
+  colored:
+    '()': pg_rad.logger.logger.ColorFormatter
+    format: '%(asctime)s - %(levelname)s: %(message)s'
 handlers:
   stdout:
     class: logging.StreamHandler
-    formatter: simple
+    formatter: colored
     stream: ext://sys.stdout
 loggers:
   root:

src/pg_rad/data/angular_efficiencies/LU_HPGe_90.csv

@@ -0,0 +1,37 @@
+angle,662,1173,1332
+0,0.015,0.030,0.033
+10,0.011,0.021,0.024
+20,0.086,0.127,0.146
+30,0.294,0.356,0.397
+40,0.661,0.700,0.734
+50,1.054,1.057,1.057
+60,1.154,1.140,1.137
+70,1.186,1.152,1.138
+80,1.151,1.114,1.097
+90,1.000,1.000,1.000
+100,1.020,1.040,1.047
+110,1.074,1.093,1.103
+120,1.113,1.092,1.102
+130,1.139,1.122,1.113
+140,1.146,1.152,1.140
+150,1.113,1.118,1.104
+160,1.113,1.096,1.099
+170,1.091,1.076,1.083
+180,1.076,1.066,1.078
+-170,1.102,1.091,1.093
+-160,1.122,1.100,1.102
+-150,1.128,1.105,1.093
+-140,1.144,1.112,1.123
+-130,1.140,1.117,1.095
+-120,1.146,1.127,1.098
+-110,1.068,1.068,1.045
+-100,1.013,1.025,1.016
+-90,1.004,1.018,1.021
+-80,1.150,1.137,1.132
+-70,1.184,1.167,1.164
+-60,1.158,1.140,1.138
+-50,1.090,1.068,1.064
+-40,0.595,0.620,0.631
+-30,0.332,0.430,0.430
+-20,0.055,0.081,0.096
+-10,0.009,0.018,0.019

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/data/backgrounds/LU_NaI_3inch.csv

@@ -0,0 +1,257 @@
+Energy,Data,cps
+0,0,0
+11.6866,0,0
+23.3732,1653,1.181558256
+35.0598,3827,2.735525375
+46.7464,5010,3.581129378
+58.433,7880,5.632594711
+70.1196,11068,7.911365261
+81.8062,13414,9.588277341
+93.4928,14536,10.39027877
+105.179,14580,10.42172981
+116.866,14138,10.10578985
+128.553,13559,9.691922802
+140.239,12725,9.095782702
+151.926,11597,8.289492495
+163.612,10696,7.645461044
+175.299,9847,7.038598999
+186.986,9063,6.478198713
+198.672,8192,5.855611151
+210.359,7627,5.451751251
+222.045,6941,4.961401001
+233.732,6450,4.610436026
+245.419,6046,4.321658327
+257.105,5560,3.974267334
+268.792,4845,3.463187991
+280.478,4267,3.05003574
+292.165,4066,2.906361687
+303.852,3711,2.652609006
+315.538,3413,2.439599714
+327.225,3080,2.201572552
+338.911,3035,2.169406719
+350.598,2789,1.993566833
+362.285,2766,1.977126519
+373.971,2432,1.73838456
+385.658,2182,1.55968549
+397.344,2010,1.436740529
+409.031,1911,1.365975697
+420.718,1871,1.337383846
+432.404,1678,1.199428163
+444.091,1601,1.144388849
+455.777,1538,1.099356683
+467.464,1473,1.052894925
+479.151,1390,0.993566833
+490.837,1369,0.978556112
+502.524,1358,0.970693352
+514.21,1356,0.96926376
+525.897,1323,0.945675482
+537.584,1256,0.897784132
+549.27,1171,0.837026447
+560.957,1072,0.766261615
+572.643,1094,0.781987134
+584.33,1157,0.827019299
+596.017,1073,0.766976412
+607.703,1083,0.774124375
+619.39,1042,0.744817727
+631.076,995,0.711222302
+642.763,938,0.670478914
+654.45,770,0.550393138
+666.136,787,0.562544675
+677.823,761,0.543959971
+689.509,703,0.502501787
+701.196,671,0.479628306
+712.883,611,0.436740529
+724.569,599,0.428162974
+736.256,573,0.40957827
+747.942,623,0.445318084
+759.629,636,0.454610436
+771.316,634,0.453180843
+783.002,602,0.430307362
+794.689,525,0.375268049
+806.375,544,0.388849178
+818.062,557,0.39814153
+829.749,527,0.376697641
+841.435,508,0.363116512
+853.122,509,0.363831308
+864.808,478,0.341672623
+876.495,498,0.355968549
+888.182,515,0.368120086
+899.868,521,0.372408863
+911.555,539,0.385275197
+923.241,546,0.390278771
+934.928,545,0.389563974
+946.615,529,0.378127234
+958.301,489,0.349535382
+969.988,490,0.350250179
+981.674,447,0.319513939
+993.361,437,0.312365976
+1005.05,406,0.290207291
+1016.73,380,0.271622588
+1028.42,357,0.255182273
+1040.11,346,0.247319514
+1051.79,348,0.248749107
+1063.48,337,0.240886347
+1075.17,376,0.268763402
+1086.85,376,0.268763402
+1098.54,363,0.259471051
+1110.23,388,0.277340958
+1121.91,330,0.235882773
+1133.6,364,0.260185847
+1145.29,323,0.230879199
+1156.97,359,0.256611866
+1168.66,330,0.235882773
+1180.35,284,0.203002144
+1192.03,256,0.182987848
+1203.72,252,0.180128663
+1215.41,243,0.173695497
+1227.09,264,0.188706219
+1238.78,226,0.16154396
+1250.47,214,0.152966405
+1262.15,184,0.131522516
+1273.84,179,0.127948535
+1285.53,168,0.120085776
+1297.21,185,0.132237312
+1308.9,177,0.126518942
+1320.59,175,0.12508935
+1332.27,171,0.122230164
+1343.96,187,0.133666905
+1355.65,206,0.147248034
+1367.33,255,0.182273052
+1379.02,331,0.23659757
+1390.71,395,0.282344532
+1402.39,540,0.385989993
+1414.08,608,0.43459614
+1425.77,612,0.437455325
+1437.45,575,0.411007863
+1449.14,516,0.368834882
+1460.82,467,0.333809864
+1472.51,369,0.263759828
+1484.2,254,0.181558256
+1495.88,220,0.157255182
+1507.57,125,0.089349535
+1519.26,136,0.097212294
+1530.94,92,0.065761258
+1542.63,74,0.052894925
+1554.32,67,0.047891351
+1566,68,0.048606147
+1577.69,70,0.05003574
+1589.38,58,0.041458184
+1601.06,66,0.047176555
+1612.75,59,0.042172981
+1624.44,47,0.033595425
+1636.12,60,0.042887777
+1647.81,74,0.052894925
+1659.5,74,0.052894925
+1671.18,84,0.060042888
+1682.87,88,0.062902073
+1694.56,88,0.062902073
+1706.24,65,0.046461758
+1717.93,78,0.05575411
+1729.62,84,0.060042888
+1741.3,62,0.04431737
+1752.99,71,0.050750536
+1764.68,54,0.038598999
+1776.36,53,0.037884203
+1788.05,53,0.037884203
+1799.74,34,0.024303074
+1811.42,36,0.025732666
+1823.11,46,0.032880629
+1834.8,45,0.032165833
+1846.48,41,0.029306648
+1858.17,33,0.023588277
+1869.86,33,0.023588277
+1881.54,29,0.020729092
+1893.23,39,0.027877055
+1904.92,37,0.026447462
+1916.6,28,0.020014296
+1928.29,41,0.029306648
+1939.98,39,0.027877055
+1951.66,46,0.032880629
+1963.35,46,0.032880629
+1975.04,36,0.025732666
+1986.72,39,0.027877055
+1998.41,47,0.033595425
+2010.1,59,0.042172981
+2021.78,56,0.040028592
+2033.47,44,0.031451036
+2045.15,57,0.040743388
+2056.84,50,0.035739814
+2068.53,68,0.048606147
+2080.21,56,0.040028592
+2091.9,53,0.037884203
+2103.59,51,0.03645461
+2115.27,30,0.021443888
+2126.96,40,0.028591851
+2138.65,47,0.033595425
+2150.33,42,0.030021444
+2162.02,49,0.035025018
+2173.71,31,0.022158685
+2185.39,29,0.020729092
+2197.08,49,0.035025018
+2208.77,26,0.018584703
+2220.45,32,0.022873481
+2232.14,23,0.016440315
+2243.83,19,0.013581129
+2255.51,31,0.022158685
+2267.2,16,0.011436741
+2278.89,21,0.015010722
+2290.57,19,0.013581129
+2302.26,20,0.014295926
+2313.95,14,0.010007148
+2325.63,18,0.012866333
+2337.32,18,0.012866333
+2349.01,22,0.015725518
+2360.69,27,0.0192995
+2372.38,21,0.015010722
+2384.07,21,0.015010722
+2395.75,35,0.02501787
+2407.44,49,0.035025018
+2419.13,55,0.039313796
+2430.81,47,0.033595425
+2442.5,68,0.048606147
+2454.19,56,0.040028592
+2465.87,57,0.040743388
+2477.56,63,0.045032166
+2489.25,52,0.037169407
+2500.93,61,0.043602573
+2512.62,37,0.026447462
+2524.31,34,0.024303074
+2535.99,34,0.024303074
+2547.68,26,0.018584703
+2559.37,15,0.010721944
+2571.05,16,0.011436741
+2582.74,9,0.006433167
+2594.43,12,0.008577555
+2606.11,8,0.00571837
+2617.8,11,0.007862759
+2629.48,2,0.001429593
+2641.17,7,0.005003574
+2652.86,5,0.003573981
+2664.54,2,0.001429593
+2676.23,2,0.001429593
+2687.92,0,0
+2699.6,6,0.004288778
+2711.29,2,0.001429593
+2722.98,6,0.004288778
+2734.66,5,0.003573981
+2746.35,5,0.003573981
+2758.04,2,0.001429593
+2769.72,5,0.003573981
+2781.41,6,0.004288778
+2793.1,7,0.005003574
+2804.78,2,0.001429593
+2816.47,2,0.001429593
+2828.16,6,0.004288778
+2839.84,5,0.003573981
+2851.53,4,0.002859185
+2863.22,1,0.000714796
+2874.9,4,0.002859185
+2886.59,4,0.002859185
+2898.28,5,0.003573981
+2909.96,2,0.001429593
+2921.65,3,0.002144389
+2933.34,4,0.002859185
+2945.02,5,0.003573981
+2956.71,3,0.002144389
+2968.4,0,0
+2980.08,0,0

src/pg_rad/data/detectors.csv

@@ -0,0 +1,4 @@
+name,type,is_isotropic
+dummy,NaI,true
+LU_NaI_3inch,NaI,true
+LU_HPGe_90,HPGe,false

src/pg_rad/data/field_efficiencies/LU_HPGe_90.csv

@@ -0,0 +1,17 @@
+energy_keV,field_efficiency_m2,uncertainty
+59.5,0.00140,0.00005
+81.0,0.00310,0.00010
+122.1,0.00420,0.00013
+136.5,0.00428,0.00017
+160.6,0.00426,0.00030
+223.2,0.00418,0.00024
+276.4,0.00383,0.00012
+302.9,0.00370,0.00012
+356.0,0.00338,0.00010
+383.8,0.00323,0.00010
+511.0,0.00276,0.00008
+661.7,0.00241,0.00007
+834.8,0.00214,0.00007
+1173.2,0.00179,0.00005
+1274.5,0.00168,0.00005
+1332.5,0.00166,0.00005

src/pg_rad/data/field_efficiencies/LU_NaI_3inch.csv

@@ -0,0 +1,64 @@
+energy_keV,field_efficiency_m2
+10,5.50129E-09
+11.22018454,2.88553E-07
+12.58925412,4.81878E-06
+14.12537545,3.55112E-05
+15.84893192,0.000146367
+17.7827941,0.000397029
+19.95262315,0.000803336
+22.38721139,0.00131657
+25.11886432,0.001862377
+28.18382931,0.002373449
+31.6227766,0.002811046
+33.164,0.00269554
+33.164,0.002698792
+35.48133892,0.002509993
+39.81071706,0.002801304
+44.66835922,0.003015877
+50.11872336,0.003227431
+56.23413252,0.00341077
+63.09573445,0.003562051
+70.79457844,0.00368852
+79.43282347,0.003788875
+89,0.003867423
+100,0.003925025
+112.2018454,0.003967222
+125.8925412,0.003991551
+141.2537545,0.004000729
+158.4893192,0.003993145
+177.827941,0.003969163
+199.5262315,0.003925289
+223.8721139,0.003856247
+251.1886432,0.00375596
+281.8382931,0.003619634
+316.227766,0.003446087
+354.8133892,0.003242691
+398.1071706,0.003021761
+446.6835922,0.002791816
+501.1872336,0.002568349
+562.3413252,0.002350052
+630.9573445,0.002147662
+707.9457844,0.001957893
+794.3282347,0.001785694
+891,0.001626634
+1000,0.001482571
+1122.018454,0.00135047
+1258.925412,0.001231358
+1412.537545,0.001116695
+1584.893192,0.001011833
+1778.27941,0.000917017
+1995.262315,0.000828435
+2238.721139,0.000746854
+2511.886432,0.000672573
+2818.382931,0.00060493
+3162.27766,0.000544458
+3548.133892,0.000488446
+3981.071706,0.000438438
+4466.835922,0.000392416
+5011.872336,0.00035092
+5623.413252,0.000313959
+6309.573445,0.000279409
+7079.457844,0.000247794
+7943.282347,0.000218768
+8913,0.000190209
+10000,0.000164309

src/pg_rad/data/field_efficiencies/dummy.csv

@@ -0,0 +1,3 @@
+energy_keV,field_efficiency_m2
+0,1.0
+10000,1.0

src/pg_rad/data/isotopes.csv

@@ -0,0 +1,6 @@
+isotope,gamma_energy_keV,branching_ratio
+Cs134,604.7,0.976
+Cs134,795.9,0.855
+Cs137,661.657,0.851
+Co60,1173.228,1.0
+Co60,1332.5,1.0

src/pg_rad/dataloader/__init__.py

@@ -1,8 +0,0 @@
-# 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

@@ -2,27 +2,21 @@ import logging
 
 import pandas as pd
 
-from pg_rad.exceptions import DataLoadError, InvalidCSVError
 
 logger = logging.getLogger(__name__)
 
+
 def load_data(filename: str) -> pd.DataFrame:
     logger.debug(f"Attempting to load file: {filename}")
 
     try:
         df = pd.read_csv(filename, delimiter=',')
-
     except FileNotFoundError as e:
-        logger.error(f"File not found: {filename}")
-        raise DataLoadError(f"File does not exist: {filename}") from e
-
+        logger.critical(e)
+        raise
     except pd.errors.ParserError as e:
-        logger.error(f"Invalid CSV format: {filename}")
-        raise InvalidCSVError(f"Invalid CSV file: {filename}") from e
-
-    except Exception as e:
-        logger.exception(f"Unexpected error while loading {filename}")
-        raise DataLoadError("Unexpected error while loading data") from e
+        logger.critical(e)
+        raise
 
     logger.debug(f"File loaded: {filename}")
     return df

src/pg_rad/detector/detector.py

@@ -0,0 +1,43 @@
+from importlib.resources import files
+
+from pandas import read_csv
+
+from pg_rad.utils.interpolators import (
+    get_field_efficiency, get_angular_efficiency
+)
+
+
+class Detector:
+    def __init__(
+        self,
+        name: str,
+        type: str,
+        is_isotropic: bool
+    ):
+        self.name = name
+        self.type = type
+        self.is_isotropic = is_isotropic
+
+    def get_efficiency(self, energy_keV, angle=None):
+        f_eff = get_field_efficiency(self.name, energy_keV)
+
+        if self.is_isotropic or angle is None:
+            return f_eff
+        else:
+            f_eff = get_field_efficiency(self.name, energy_keV)
+            a_eff = get_angular_efficiency(self.name, energy_keV, *angle)
+            return f_eff * a_eff
+
+
+def load_detector(name) -> Detector:
+    csv = files('pg_rad.data').joinpath('detectors.csv')
+    data = read_csv(csv)
+    dets = data['name'].values
+    if name in dets:
+        row = data[data['name'] == name].iloc[0]
+        return Detector(row['name'], row['type'], row['is_isotropic'])
+    else:
+        raise NotImplementedError(
+            f"Detector with name '{name}' not in detector library. Available:"
+            f"{', '.join(dets)}"
+        )

src/pg_rad/exceptions/__init__.py

@@ -1,10 +0,0 @@
-# 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,45 @@
 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."""
+
+
+class InvalidYAMLError(DataLoadError):
+    """Raised when a file is not a valid YAML."""
+
+
+class OutOfBoundsError(Exception):
+    """Raised when an object is attempted to be placed out of bounds."""
+
+
+class MissingConfigKeyError(KeyError):
+    """Raised when a (nested) config key is missing in the config."""
+    def __init__(self, key, subkey=None):
+        if subkey:
+            if isinstance(subkey, str):
+                pass
+            elif isinstance(subkey, set):
+                subkey = ', '.join(list(subkey))
+            self.message = f"Missing key in {key}: {subkey}"
+        else:
+            self.message = f"Missing key: {key}"
+
+        super().__init__(self.message)
+
+
+class DimensionError(ValueError):
+    """Raised if dimensions or coordinates do not match the system."""
+
+
+class InvalidIsotopeError(ValueError):
+    """Raised if attempting to load an isotope that is not valid."""
+
+
+class InvalidConfigValueError(ValueError):
+    """Raised if a config key has an incorrect type or value."""

src/pg_rad/inputparser/parser.py

@@ -0,0 +1,367 @@
+import logging
+import os
+from typing import Any, Dict, List, Union
+
+import yaml
+
+from pg_rad.exceptions.exceptions import (
+    MissingConfigKeyError,
+    DimensionError,
+    InvalidConfigValueError,
+    InvalidYAMLError
+)
+from pg_rad.configs import defaults
+from pg_rad.isotopes.isotope import get_isotope
+
+from .specs import (
+    MetadataSpec,
+    RuntimeSpec,
+    SimulationOptionsSpec,
+    PathSpec,
+    ProceduralPathSpec,
+    CSVPathSpec,
+    PointSourceSpec,
+    AbsolutePointSourceSpec,
+    RelativePointSourceSpec,
+    DetectorSpec,
+    SimulationSpec
+)
+
+
+logger = logging.getLogger(__name__)
+
+
+class ConfigParser:
+    _ALLOWED_ROOT_KEYS = {
+        "name",
+        "speed",
+        "acquisition_time",
+        "path",
+        "sources",
+        "detector",
+        "options"
+    }
+
+    def __init__(self, config_source: str):
+        self.config = self._load_yaml(config_source)
+
+    def parse(self) -> SimulationSpec:
+        self._warn_unknown_keys(
+            section="root",
+            provided=set(self.config.keys()),
+            allowed=self._ALLOWED_ROOT_KEYS,
+        )
+
+        metadata = self._parse_metadata()
+        runtime = self._parse_runtime()
+        options = self._parse_options()
+        path = self._parse_path()
+        sources = self._parse_point_sources()
+        detector = self._parse_detector()
+
+        return SimulationSpec(
+            metadata=metadata,
+            runtime=runtime,
+            options=options,
+            path=path,
+            point_sources=sources,
+            detector=detector
+        )
+
+    def _load_yaml(self, config_source: str) -> Dict[str, Any]:
+        if os.path.exists(config_source):
+            with open(config_source) as f:
+                data = yaml.safe_load(f)
+        else:
+            data = yaml.safe_load(config_source)
+
+        if not isinstance(data, dict):
+            raise InvalidYAMLError(
+                "Provided path or string is not a valid YAML representation."
+            )
+
+        return data
+
+    def _parse_metadata(self) -> MetadataSpec:
+        try:
+            return MetadataSpec(
+                name=self.config["name"]
+            )
+        except KeyError as e:
+            raise MissingConfigKeyError("root", {"name"}) from e
+
+    def _parse_runtime(self) -> RuntimeSpec:
+        required = {"speed", "acquisition_time"}
+        missing = required - self.config.keys()
+        if missing:
+            raise MissingConfigKeyError("root", missing)
+
+        return RuntimeSpec(
+            speed=float(self.config["speed"]),
+            acquisition_time=float(self.config["acquisition_time"]),
+        )
+
+    def _parse_options(self) -> SimulationOptionsSpec:
+        options = self.config.get("options", {})
+
+        allowed = {"air_density_kg_per_m3", "seed"}
+        self._warn_unknown_keys(
+            section="options",
+            provided=set(options.keys()),
+            allowed=allowed,
+        )
+
+        air_density = options.get(
+            "air_density_kg_per_m3",
+            defaults.DEFAULT_AIR_DENSITY
+        )
+        seed = options.get("seed")
+
+        if not isinstance(air_density, float) or air_density <= 0:
+            raise InvalidConfigValueError(
+                "options.air_density_kg_per_m3 must be a positive float "
+                "in kg/m^3."
+            )
+        if (
+            seed is not None or
+            (isinstance(seed, int) and seed <= 0)
+        ):
+            raise InvalidConfigValueError(
+                "Seed must be a positive integer value."
+            )
+
+        return SimulationOptionsSpec(
+            air_density=air_density,
+            seed=seed,
+        )
+
+    def _parse_path(self) -> PathSpec:
+        allowed_csv = {"file", "east_col_name", "north_col_name", "z"}
+        allowed_proc = {"segments", "length", "z", "alpha", "direction"}
+
+        path = self.config.get("path")
+
+        direction = path.get("direction", 'positive')
+
+        if direction == 'positive':
+            opposite_direction = False
+        elif direction == 'negative':
+            opposite_direction = True
+        else:
+            raise InvalidConfigValueError(
+                "Direction must be positive or negative."
+            )
+
+        if path is None:
+            raise MissingConfigKeyError("path")
+
+        if not isinstance(path, dict):
+            raise InvalidConfigValueError("Path must be a dictionary.")
+
+        if "file" in path:
+            self._warn_unknown_keys(
+                section="path (csv)",
+                provided=set(path.keys()),
+                allowed=allowed_csv,
+            )
+
+            return CSVPathSpec(
+                file=path["file"],
+                east_col_name=path["east_col_name"],
+                north_col_name=path["north_col_name"],
+                z=path.get("z", 0),
+                opposite_direction=opposite_direction
+            )
+
+        elif "segments" in path:
+            self._warn_unknown_keys(
+                section="path (procedural)",
+                provided=set(path.keys()),
+                allowed=allowed_proc,
+            )
+            return self._parse_procedural_path(path, opposite_direction)
+
+        else:
+            raise InvalidConfigValueError("Invalid path configuration.")
+
+    def _parse_procedural_path(
+        self,
+        path: Dict[str, Any],
+        opposite_direction: bool
+    ) -> ProceduralPathSpec:
+        raw_segments = path.get("segments")
+
+        if not isinstance(raw_segments, List):
+            raise InvalidConfigValueError(
+                "path.segments must be a list of segments."
+            )
+
+        raw_length = path.get("length")
+
+        if raw_length is None:
+            raise MissingConfigKeyError("path", {"length"})
+
+        if isinstance(raw_length, int | float):
+            raw_length = [float(raw_length)]
+
+        segments, angles = self._process_segment_angles(raw_segments)
+        lengths = self._process_segment_lengths(raw_length, len(segments))
+
+        return ProceduralPathSpec(
+            segments=segments,
+            angles=angles,
+            lengths=lengths,
+            z=path.get("z", defaults.DEFAULT_PATH_HEIGHT),
+            alpha=path.get("alpha", defaults.DEFAULT_ALPHA),
+            opposite_direction=opposite_direction
+        )
+
+    def _process_segment_angles(
+        self,
+        raw_segments: List[Union[str, dict]]
+    ) -> List[Dict[str, Any]]:
+
+        segments, angles = [], []
+
+        for segment in raw_segments:
+
+            if isinstance(segment, str):
+                segments.append(segment)
+                angles.append(None)
+
+            elif isinstance(segment, dict):
+                if len(segment) != 1:
+                    raise InvalidConfigValueError(
+                        "Invalid segment definition."
+                    )
+                seg_type, angle = list(segment.items())[0]
+                segments.append(seg_type)
+                angles.append(angle)
+
+            else:
+                raise InvalidConfigValueError(
+                    "Invalid segment entry format."
+                )
+
+        return segments, angles
+
+    def _process_segment_lengths(
+        self,
+        raw_length_list: List[Union[int, float]],
+        num_segments: int
+    ) -> List[float]:
+        num_lengths = len(raw_length_list)
+
+        if num_lengths == num_segments or num_lengths == 1:
+            return raw_length_list
+        else:
+            raise ValueError(
+                "Path length must either be a single number or a list with "
+                "number of elements equal to the number of segments."
+            )
+
+    @staticmethod
+    def _is_turn(segment_type: str) -> bool:
+        return segment_type in {"turn_left", "turn_right"}
+
+    def _parse_point_sources(self) -> List[PointSourceSpec]:
+        source_dict = self.config.get("sources")
+        if source_dict is None:
+            raise MissingConfigKeyError("sources")
+
+        if not isinstance(source_dict, dict):
+            raise InvalidConfigValueError(
+                "sources must have subkeys representing point source names."
+            )
+
+        specs: List[PointSourceSpec] = []
+
+        for name, params in source_dict.items():
+            required = {
+                "activity_MBq", "isotope", "position", "gamma_energy_keV"
+            }
+            if not isinstance(params, dict):
+                raise InvalidConfigValueError(
+                    f"sources.{name} is not defined correctly."
+                    f" Must have subkeys {required}"
+                )
+
+            missing = required - params.keys()
+            if missing:
+                raise MissingConfigKeyError(name, missing)
+
+            activity = params.get("activity_MBq")
+            isotope_name = params.get("isotope")
+            gamma_energy_keV = params.get("gamma_energy_keV")
+
+            isotope = get_isotope(isotope_name, gamma_energy_keV)
+
+            if not isinstance(activity, int | float) or activity <= 0:
+                raise InvalidConfigValueError(
+                    f"sources.{name}.activity_MBq must be positive value "
+                    "in MegaBequerels."
+                )
+
+            position = params.get("position")
+
+            if isinstance(position, list):
+                if len(position) != 3:
+                    raise DimensionError(
+                        "Absolute position must be [x, y, z]."
+                    )
+
+                specs.append(
+                    AbsolutePointSourceSpec(
+                        name=name,
+                        activity_MBq=float(activity),
+                        isotope=isotope,
+                        x=float(position[0]),
+                        y=float(position[1]),
+                        z=float(position[2]),
+                    )
+                )
+
+            elif isinstance(position, dict):
+                alignment = position.get("acquisition_alignment")
+                if alignment not in {'best', 'worst', None}:
+                    raise InvalidConfigValueError(
+                        f"sources.{name}.acquisition_alignment must be "
+                        "'best' or 'worst', with 'best' aligning source "
+                        f"{name} in the middle of the two nearest acquisition "
+                        "points, and 'worst' aligning exactly perpendicular "
+                        "to the nearest acquisition point."
+                    )
+
+                specs.append(
+                    RelativePointSourceSpec(
+                        name=name,
+                        activity_MBq=float(activity),
+                        isotope=isotope,
+                        along_path=float(position["along_path"]),
+                        dist_from_path=float(position["dist_from_path"]),
+                        side=position["side"],
+                        z=position.get("z", defaults.DEFAULT_SOURCE_HEIGHT),
+                        alignment=alignment
+                    )
+                )
+
+            else:
+                raise InvalidConfigValueError(
+                    f"Invalid position format for source '{name}'."
+                    )
+
+        return specs
+
+    def _parse_detector(self) -> DetectorSpec:
+        det_name = self.config.get("detector")
+        if not det_name:
+            raise MissingConfigKeyError("detector")
+
+        return DetectorSpec(name=det_name)
+
+    def _warn_unknown_keys(self, section: str, provided: set, allowed: set):
+        unknown = provided - allowed
+        if unknown:
+            logger.warning(
+                f"Unknown keys in '{section}' section: {unknown}"
+            )

src/pg_rad/inputparser/specs.py

@@ -0,0 +1,79 @@
+from abc import ABC
+from dataclasses import dataclass
+from typing import Literal
+
+
+@dataclass
+class MetadataSpec:
+    name: str
+
+
+@dataclass
+class RuntimeSpec:
+    speed: float
+    acquisition_time: float
+
+
+@dataclass
+class SimulationOptionsSpec:
+    air_density: float
+    seed: int | None = None
+
+
+@dataclass
+class PathSpec(ABC):
+    z: int | float
+    opposite_direction: bool
+
+
+@dataclass
+class ProceduralPathSpec(PathSpec):
+    segments: list[str]
+    angles: list[float]
+    lengths: list[int | None]
+    alpha: float
+
+
+@dataclass
+class CSVPathSpec(PathSpec):
+    file: str
+    east_col_name: str
+    north_col_name: str
+
+
+@dataclass
+class PointSourceSpec(ABC):
+    activity_MBq: float
+    isotope: str
+    name: str
+
+
+@dataclass
+class AbsolutePointSourceSpec(PointSourceSpec):
+    x: float
+    y: float
+    z: float
+
+
+@dataclass
+class RelativePointSourceSpec(PointSourceSpec):
+    along_path: float
+    dist_from_path: float
+    side: str
+    z: float
+    alignment: Literal["best", "worst"] | None
+
+
+@dataclass
+class DetectorSpec:
+    name: str
+
+
+@dataclass
+class SimulationSpec:
+    metadata: MetadataSpec
+    runtime: RuntimeSpec
+    options: SimulationOptionsSpec
+    path: PathSpec
+    point_sources: list[PointSourceSpec]
+    detector: DetectorSpec

src/pg_rad/isotopes/__init__.py

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

src/pg_rad/isotopes/isotope.py

@@ -1,3 +1,12 @@
+from importlib.resources import files
+
+from pandas import read_csv
+
+from pg_rad.configs.filepaths import ISOTOPE_TABLE
+from pg_rad.exceptions.exceptions import InvalidIsotopeError
+from pg_rad.utils.interpolators import get_mass_attenuation_coeff
+
+
 class Isotope:
     """Represents the essential information of an isotope.
 
@@ -16,8 +25,43 @@ class Isotope:
             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)
+
+
+def get_isotope(isotope: str, energy_gamma_keV: float) -> Isotope:
+    """Lazy factory function to create isotope objects."""
+    path = files('pg_rad.data').joinpath(ISOTOPE_TABLE)
+    df = read_csv(path)
+
+    isotope_df = df[df['isotope'] == isotope]
+
+    if isotope_df.empty:
+        raise InvalidIsotopeError(f"No data available for isotope {isotope}.")
+
+    tol = 0.01 * energy_gamma_keV
+    closest_energy = isotope_df[
+        (isotope_df['gamma_energy_keV'] >= energy_gamma_keV - tol) &
+        (isotope_df['gamma_energy_keV'] <= energy_gamma_keV + tol)
+    ]
+
+    if closest_energy.empty:
+        available_gammas = ', '.join(
+            str(x)+' keV' for x in isotope_df['gamma_energy_keV'].to_list()
+        )
+        raise InvalidIsotopeError(
+            f"No gamma of {energy_gamma_keV}±{tol} keV "
+            f"found for isotope {isotope}. "
+            f"Available gammas are: {available_gammas}"
+        )
+
+    matched_row = closest_energy.iloc[0]
+    return Isotope(
+        name=isotope,
+        E=matched_row['gamma_energy_keV'],
+        b=matched_row['branching_ratio']
+    )

src/pg_rad/landscape/__init__.py

@@ -1,8 +0,0 @@
-# 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/builder.py

@@ -0,0 +1,243 @@
+import logging
+from typing import Literal, Self
+
+import numpy as np
+
+from .landscape import Landscape
+from pg_rad.dataloader.dataloader import load_data
+from pg_rad.objects.sources import PointSource
+from pg_rad.exceptions.exceptions import OutOfBoundsError
+from pg_rad.utils.projection import rel_to_abs_source_position
+from pg_rad.inputparser.specs import (
+    SimulationSpec,
+    CSVPathSpec,
+    AbsolutePointSourceSpec,
+    RelativePointSourceSpec,
+    DetectorSpec
+)
+from pg_rad.detector.detector import load_detector
+from pg_rad.path.path import Path, path_from_RT90
+
+from road_gen.generators.segmented_road_generator import SegmentedRoadGenerator
+
+
+logger = logging.getLogger(__name__)
+
+
+class LandscapeBuilder:
+    def __init__(self, name: str = "Unnamed landscape"):
+        self.name = name
+        self._path = None
+        self._point_sources = []
+        self._size = None
+        self._air_density = 1.243
+        self._detector = None
+
+        logger.debug(f"LandscapeBuilder initialized: {self.name}")
+
+    def set_air_density(self, air_density: float | None) -> Self:
+        """Set the air density of the world."""
+        if air_density and isinstance(air_density, float):
+            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 get_path(self):
+        return self._path
+
+    def set_path_from_segments(
+        self,
+        sim_spec: SimulationSpec
+    ):
+        path = sim_spec.path
+        segments = path.segments
+        lengths = path.lengths
+        angles = path.angles
+        alpha = path.alpha
+        z = path.z
+
+        sg = SegmentedRoadGenerator(
+            ds=sim_spec.runtime.speed * sim_spec.runtime.acquisition_time,
+            velocity=sim_spec.runtime.speed,
+            seed=sim_spec.options.seed
+        )
+
+        x, y = sg.generate(
+            segments=segments,
+            lengths=lengths,
+            angles=angles,
+            alpha=alpha
+        )
+
+        self._path = Path(
+            list(zip(x, y)),
+            z=z,
+            opposite_direction=path.opposite_direction
+        )
+        self._fit_landscape_to_path()
+        return self
+
+    def set_path_from_experimental_data(
+        self,
+        spec: CSVPathSpec
+    ) -> Self:
+        df = load_data(spec.file)
+        self._path = path_from_RT90(
+            df=df,
+            east_col=spec.east_col_name,
+            north_col=spec.north_col_name,
+            z=spec.z,
+            opposite_direction=spec.opposite_direction
+        )
+
+        self._fit_landscape_to_path()
+
+        return self
+
+    def set_point_sources(
+        self,
+        *sources: AbsolutePointSourceSpec | RelativePointSourceSpec
+    ):
+        """Add one or more point sources to the world.
+
+        Args:
+            *sources (AbsolutePointSourceSpec | RelativePointSourceSpec):
+                One or more sources, passed as SourceSpecs tuple
+        Raises:
+            OutOfBoundsError: If any source is outside the boundaries of the
+            landscape.
+        """
+
+        for s in sources:
+            if isinstance(s, AbsolutePointSourceSpec):
+                pos = (s.x, s.y, s.z)
+            elif isinstance(s, RelativePointSourceSpec):
+                path = self.get_path()
+
+                if s.alignment:
+                    along_path = self._align_relative_source(
+                        s.along_path,
+                        path,
+                        s.alignment
+                    )
+                    logger.info(
+                        f"Because source {s.name} was set to align with path "
+                        f"({s.alignment} alignment), it was moved to be at "
+                        f"{along_path} m along the path from {s.along_path} m."
+                    )
+                else:
+                    along_path = s.along_path
+                pos = rel_to_abs_source_position(
+                    x_list=path.x_list,
+                    y_list=path.y_list,
+                    path_z=path.z,
+                    along_path=along_path,
+                    side=s.side,
+                    dist_from_path=s.dist_from_path)
+            if any(
+                p < 0 or p >= s for p, s in zip(pos, self._size)
+            ):
+                raise OutOfBoundsError(
+                    "One or more sources attempted to "
+                    "be placed outside the landscape."
+                    )
+
+            self._point_sources.append(PointSource(
+                activity_MBq=s.activity_MBq,
+                isotope=s.isotope,
+                position=pos,
+                name=s.name
+            ))
+
+    def set_detector(self, spec: DetectorSpec) -> Self:
+        self._detector = load_detector(spec.name)
+        return self
+
+    def _fit_landscape_to_path(self) -> None:
+        """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."
+                )
+
+            max_size = max(self._path.size)
+            self.set_landscape_size((max_size, max_size))
+
+    def _align_relative_source(
+        self,
+        along_path: float,
+        path: "Path",
+        mode: Literal["best", "worst"],
+    ) -> tuple[float, float, float]:
+        """Given the arc length at which the point source is placed,
+            align the source relative to the waypoints of the path. Here,
+            'best' means the point source is moved such that it is
+            perpendicular to the midpoint between two acuisition points.
+            'worst' means the point source is moved such that it is
+            perpendicular to the nearest acquisition point.
+
+            The distance to the path is not affected by this algorithm.
+
+            For more details on alignment, see
+            Fig. 4 and page 24 in Bukartas (2021).
+
+        Args:
+            along_path (float): Current arc length position of the source.
+            path (Path): The path to align to.
+            mode (Literal["best", "worst"]): Alignment mode.
+
+        Returns:
+            along_new (float): The updated arc length position.
+        """
+        ds = np.hypot(
+            path.x_list[1] - path.x_list[0],
+            path.y_list[1] - path.y_list[0],
+        )
+
+        if mode == "worst":
+            along_new = round(along_path / ds) * ds
+
+        elif mode == "best":
+            along_new = (round(along_path / ds - 0.5) + 0.5) * ds
+
+        else:
+            raise ValueError(f"Unknown alignment mode: {mode}")
+
+        return along_new
+
+    def build(self):
+        landscape = Landscape(
+            name=self.name,
+            path=self._path,
+            point_sources=self._point_sources,
+            detector=self._detector,
+            size=self._size,
+            air_density=self._air_density
+        )
+
+        logger.info(f"Landscape built successfully: {landscape.name}")
+        return landscape

src/pg_rad/landscape/director.py

@@ -0,0 +1,50 @@
+from importlib.resources import files
+import logging
+
+from pg_rad.configs.filepaths import TEST_EXP_DATA
+from .builder import LandscapeBuilder
+from pg_rad.inputparser.specs import (
+    SimulationSpec,
+    CSVPathSpec,
+    ProceduralPathSpec
+)
+from pg_rad.objects.sources import PointSource
+
+
+logger = logging.getLogger(__name__)
+
+
+class LandscapeDirector:
+    def __init__(self):
+        logger.debug("LandscapeDirector initialized.")
+
+    @staticmethod
+    def build_test_landscape():
+        fp = files('pg_rad.data').joinpath(TEST_EXP_DATA)
+        source = PointSource(
+            activity_MBq=100E6,
+            isotope="CS137",
+            position=(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
+
+    @staticmethod
+    def build_from_config(config: SimulationSpec):
+        lb = LandscapeBuilder(config.metadata.name)
+        lb.set_air_density(config.options.air_density)
+
+        if isinstance(config.path, CSVPathSpec):
+            lb.set_path_from_experimental_data(spec=config.path)
+        elif isinstance(config.path, ProceduralPathSpec):
+            lb.set_path_from_segments(
+                sim_spec=config,
+            )
+        lb.set_point_sources(*config.point_sources)
+        lb.set_detector(config.detector)
+        landscape = lb.build()
+        return landscape

src/pg_rad/landscape/landscape.py

@@ -1,131 +1,42 @@
 import logging
 
-from matplotlib import pyplot as plt
-from matplotlib.patches import Circle
-import numpy as np
+from pg_rad.path.path import Path
+from pg_rad.detector.detector import Detector
+from pg_rad.objects.sources import PointSource
 
-from pg_rad.path import Path
-from pg_rad.objects import PointSource
 
 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'.
+class Landscape:
+    """
+    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 = 'meters',
+            name: str,
+            path: Path,
+            air_density: float,
+            point_sources: list[PointSource],
+            detector: Detector,
+            size: tuple[int, int, int]
             ):
-        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.")
-
-        self.air_density = air_density
-        self.scale = scale
-
-        self.path: Path = None
-        self.sources: list[PointSource] = []
-        logger.debug("Landscape initialized.")
-    
-    def plot(self, z = 0):
-        """Plot a slice of the world at a height `z`.
+        """Initialize a landscape.
 
         Args:
-            z (int, optional): Height of slice. Defaults to 0.
-
-        Returns:
-            fig, ax: Matplotlib figure objects.
-        """        
-        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}]")
-
-        if not self.path == 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(
-                        (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)
-                
-        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, ...
-        Raises:
-            ValueError: If the source is outside the boundaries of the landscape.
+            name (str): Name of the landscape.
+            path (Path): The path of the detector.
+            air_density (float): Air density in kg/m^3.
+            point_sources (list[PointSource]): List of point sources.
+            detector (Detector): The detector object.
+            size (tuple[int, int, int]): Size of the world.
         """
 
-        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)
-            for source in sources
-        ):
-            raise ValueError("One or more sources are outside the landscape!")
-
-        self.sources.extend(sources)
-
-    def set_path(self, path: Path):
-        """
-        Set the path in the landscape.
-        """
+        self.name = name
         self.path = path
+        self.point_sources = point_sources
+        self.size = size
+        self.air_density = air_density
+        self.detector = detector
 
-def create_landscape_from_path(path: Path, max_z = 500):
-    """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.
-
-    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.debug(f"Landscape created: {self.name}")

src/pg_rad/logger/logger.py

@@ -0,0 +1,40 @@
+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"]
+
+    if log_level not in levels:
+        raise ValueError(f"Log level must be one of {levels}.")
+
+    config_file = files('pg_rad.configs').joinpath(LOGGING_CONFIG)
+
+    with open(config_file) as f:
+        config = yaml.safe_load(f)
+
+    config["loggers"]["root"]["level"] = log_level
+
+    logging.config.dictConfig(config)
+
+
+class ColorFormatter(logging.Formatter):
+    # ANSI escape codes
+    COLORS = {
+        logging.DEBUG: "\033[36m",     # Cyan
+        logging.INFO: "\033[32m",      # Green
+        logging.WARNING: "\033[33m",   # Yellow
+        logging.ERROR: "\033[31m",     # Red
+        logging.CRITICAL: "\033[41m",  # Red background
+    }
+    RESET = "\033[0m"
+
+    def format(self, record):
+        color = self.COLORS.get(record.levelno, self.RESET)
+        record.levelname = f"{color}{record.levelname}{self.RESET}"
+        record.msg = f"{record.msg}"
+        return super().format(record)

src/pg_rad/logging/__init__.py

@@ -1,5 +0,0 @@
-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,20 +0,0 @@
-import logging
-import pathlib
-
-import yaml
-
-def setup_logger(log_level: str = "WARNING"):
-    levels = ["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"]
-    
-    if not log_level 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)

src/pg_rad/main.py

@@ -0,0 +1,133 @@
+import argparse
+import logging
+import sys
+
+from pandas.errors import ParserError
+
+from pg_rad.exceptions.exceptions import (
+    MissingConfigKeyError,
+    OutOfBoundsError,
+    DimensionError,
+    InvalidConfigValueError,
+    InvalidIsotopeError,
+    InvalidYAMLError
+)
+from pg_rad.logger.logger import setup_logger
+from pg_rad.inputparser.parser import ConfigParser
+from pg_rad.landscape.director import LandscapeDirector
+from pg_rad.plotting.result_plotter import ResultPlotter
+from pg_rad.simulator.engine import SimulationEngine
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        prog="pg-rad",
+        description="Primary Gamma RADiation landscape tool"
+    )
+    parser.add_argument(
+        "--config",
+        help="Build from a config file."
+    )
+    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"],
+    )
+    parser.add_argument(
+        "--saveplot",
+        action="store_true",
+        help="Save the plot or not."
+    )
+
+    args = parser.parse_args()
+    setup_logger(args.loglevel)
+    logger = logging.getLogger(__name__)
+
+    if args.test:
+        test_yaml = """
+        name: Test landscape
+        speed: 8.33
+        acquisition_time: 1
+
+        path:
+            length:
+                - 500
+                - 500
+            segments:
+                - straight
+                - turn_left: 45
+            direction: negative
+
+        sources:
+            test_source:
+                activity_MBq: 100
+                position: [250, 100, 0]
+                isotope: Cs137
+                gamma_energy_keV: 661
+
+        detector: LU_NaI_3inch
+        """
+
+        cp = ConfigParser(test_yaml).parse()
+        landscape = LandscapeDirector.build_from_config(cp)
+        output = SimulationEngine(
+            landscape=landscape,
+            runtime_spec=cp.runtime,
+            sim_spec=cp.options,
+        ).simulate()
+
+        plotter = ResultPlotter(landscape, output)
+        plotter.plot()
+
+    elif args.config:
+        try:
+            cp = ConfigParser(args.config).parse()
+            landscape = LandscapeDirector.build_from_config(cp)
+            output = SimulationEngine(
+                landscape=landscape,
+                runtime_spec=cp.runtime,
+                sim_spec=cp.options
+            ).simulate()
+
+            plotter = ResultPlotter(landscape, output)
+            plotter.plot()
+        except (
+            MissingConfigKeyError,
+            KeyError
+        ) as e:
+            logger.critical(e)
+            logger.critical(
+                "The config file is missing required keys or may be an "
+                "invalid YAML file. Check the log above. Consult the "
+                "documentation for examples of how to write a config file."
+                )
+            sys.exit(1)
+        except (
+            OutOfBoundsError,
+            DimensionError,
+            InvalidIsotopeError,
+            InvalidConfigValueError
+        ) as e:
+            logger.critical(e)
+            logger.critical(
+                "One or more items in config are not specified correctly. "
+                "Please consult this log and fix the problem."
+                )
+            sys.exit(1)
+
+        except (
+            FileNotFoundError,
+            ParserError,
+            InvalidYAMLError
+        ) as e:
+            logger.critical(e)
+            sys.exit(1)
+
+
+if __name__ == "__main__":
+    main()

src/pg_rad/objects/__init__.py

@@ -1,13 +0,0 @@
-# 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 (Object,)
-from pg_rad.objects.sources import (PointSource,)
-
-__all__ = ['Detector', 'Object', 'PointSource', 'detectors', 'objects',
-           'sources']

src/pg_rad/objects/objects.py

@@ -1,33 +1,44 @@
-import math
 from typing import Self
 
-class Object:
+import numpy as np
+
+from pg_rad.exceptions.exceptions import DimensionError
+
+
+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.
-            name (str, optional): Name for the object. Defaults to "Unnamed object".
-            color (str, optional): Matplotlib compatible color string. Defaults to "red".
+            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 DimensionError("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

src/pg_rad/objects/sources.py

@@ -1,33 +1,33 @@
 import logging
 
-from .objects import Object
-from pg_rad.isotopes import Isotope
+from .objects import BaseObject
+from pg_rad.isotopes.isotope import Isotope
 
 logger = logging.getLogger(__name__)
 
-class PointSource(Object):
+
+class PointSource(BaseObject):
     _id_counter = 1
+
     def __init__(
             self,
-            x: float,
-            y: float,
-            z: float,
-            activity: int,
+            activity_MBq: int,
             isotope: Isotope,
+            position: 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.
+            activity_MBq (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
@@ -37,13 +37,18 @@ class PointSource(Object):
         if name is None:
             name = f"Source {self.id}"
 
-        super().__init__(x, y, z, name, color)
+        super().__init__(position, name, color)
 
-        self.activity = activity
+        self.activity = activity_MBq
         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)"
+        x, y, z = self.position
+        repr_str = (f"PointSource(name={self.name}, "
+                    + f"pos={(x, y, z)}, "
+                    + f"A={self.activity} MBq), "
+                    + f"isotope={self.isotope.name}.")
+
+        return repr_str

src/pg_rad/path/__init__.py

@@ -1,9 +0,0 @@
-# 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

@@ -5,12 +5,11 @@ import math
 from matplotlib import pyplot as plt
 import numpy as np
 import pandas as pd
-import piecewise_regression
 
-from pg_rad.exceptions import ConvergenceError
 
 logger = logging.getLogger(__name__)
 
+
 class PathSegment:
     def __init__(self, a: tuple[float, float], b: tuple[float, float]):
         """A straight Segment of a Path, from (x_a, y_a) to (x_b, y_b).
@@ -38,40 +37,51 @@ class PathSegment:
         else:
             raise IndexError
 
+
 class Path:
     def __init__(
             self,
             coord_list: Sequence[tuple[float, float]],
-            z: float = 0,
-            path_simplify = False
+            opposite_direction: bool,
+            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.
-            path_simplify (bool, optional): Whether to pg_rad.path.simplify_path(). Defaults to False.
+            coord_list (Sequence[tuple[float, float]]): List of x,y
+            coordinates.
+            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:
-            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))
-            except ConvergenceError:
-                logger.warning("Continuing without simplifying path.")
-
         self.x_list = list(x)
         self.y_list = list(y)
 
         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
+        self.opposite_direction = opposite_direction
+
+        self.size = (
+            np.ceil(max(self.x_list)),
+            np.ceil(max(self.y_list)),
+            z + z_box
+        )
 
         logger.debug("Path created.")
 
@@ -92,76 +102,6 @@ class Path:
         """
         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.
-
-    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:
-    "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."
-
-    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.
-        n_breakpoints (int, optional): Number of breakpoints. Defaults to 3.
-
-    Returns:
-        x (list[float]): Reduced list of x coordinates.
-        y (list[float]): Reduced list of y coordinates.
-
-    Raises:
-        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.
-    """
-    
-    logger.debug(f"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.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)
-
-    if keep_endpoints_equal:
-        logger.debug("Forcing endpoint equality.")
-        y_points[0] = y[0]
-        y_points[-1] = y[-1]
-
-    logger.info(
-        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,
@@ -169,15 +109,18 @@ def path_from_RT90(
         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])

src/pg_rad/physics/fluence.py

@@ -0,0 +1,157 @@
+from typing import Tuple, TYPE_CHECKING
+
+import numpy as np
+
+if TYPE_CHECKING:
+    from pg_rad.landscape.landscape import Landscape
+    from pg_rad.detector.detector import Detector
+
+
+def phi(
+        r: float,
+        activity: float | int,
+        branching_ratio: float,
+        mu_mass_air: float,
+        air_density: float,
+        eff: 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_air = 0.1 * mu_mass_air * air_density
+
+    phi_r = (
+        activity
+        * eff
+        * branching_ratio
+        * np.exp(-mu_air * r)
+    ) / (4 * np.pi * r**2)
+
+    return phi_r
+
+
+def calculate_count_rate_per_second(
+    landscape: "Landscape",
+    pos: np.ndarray,
+    detector: "Detector",
+    scaling=1E6
+):
+    """Compute count rate in s^-1 m^-2 at a position in the landscape.
+
+    Args:
+        landscape (Landscape): The landscape to compute.
+        pos (np.ndarray): (N, 3) array of positions.
+        detector (Detector):
+            Detector object, needed to compute correct efficiency.
+
+    Returns:
+        total_phi (np.ndarray): (N,) array of count rates per second.
+    """
+    pos = np.atleast_2d(pos)
+    total_phi = np.zeros(pos.shape[0])
+
+    for source in landscape.point_sources:
+        # See Bukartas (2021) page 25 for incidence angle math
+        source_to_detector = pos - np.array(source.pos)
+        r = np.linalg.norm(source_to_detector, axis=1)
+        r = np.maximum(r, 1E-3)  # enforce minimum distance of 1cm
+
+        if not detector.is_isotropic:
+            v = np.zeros_like(pos)
+            v[1:] = pos[1:] - pos[:-1]
+            v[0] = v[1]  # handle first point
+            vx, vy = v[:, 0], v[:, 1]
+
+            r_vec = pos - np.array(source.pos)
+            rx, ry = r_vec[:, 0], r_vec[:, 1]
+
+            theta = np.arctan2(vy, vx) - np.arctan2(ry, rx)
+
+            # normalise to [-pi, pi] and convert to degrees
+            theta = (theta + np.pi) % (2 * np.pi) - np.pi
+            theta_deg = np.degrees(theta)
+
+            eff = detector.get_efficiency(source.isotope.E, theta_deg)
+        else:
+            eff = detector.get_efficiency(source.isotope.E)
+
+        phi_source = phi(
+            r=r,
+            activity=source.activity * scaling,
+            branching_ratio=source.isotope.b,
+            mu_mass_air=source.isotope.mu_mass_air,
+            air_density=landscape.air_density,
+            eff=eff
+        )
+
+        total_phi += phi_source
+
+    return total_phi
+
+
+def calculate_counts_along_path(
+    landscape: "Landscape",
+    detector: "Detector",
+    velocity: float,
+    points_per_segment: int = 10,
+) -> Tuple[np.ndarray, np.ndarray]:
+    """Compute the counts recorded in each acquisition period in the landscape.
+
+    Args:
+        landscape (Landscape): _description_
+        detector (Detector): _description_
+        points_per_segment (int, optional): _description_. Defaults to 100.
+
+    Returns:
+        Tuple[np.ndarray, np.ndarray]: Array of acquisition points and
+            integrated count rates.
+    """
+
+    path = landscape.path
+    num_points = len(path.x_list)
+    num_segments = len(path.segments)
+
+    segment_lengths = np.array([seg.length for seg in path.segments])
+    original_distances = np.zeros(num_points)
+    original_distances[1:] = np.cumsum(segment_lengths)
+
+    # arc lengths at which to evaluate the path
+    total_subpoints = num_segments * points_per_segment
+    s = np.linspace(0, original_distances[-1], total_subpoints)
+
+    # Interpolate x and y as functions of arc length
+    xnew = np.interp(s, original_distances, path.x_list)
+    ynew = np.interp(s, original_distances, path.y_list)
+    z = np.full_like(xnew, path.z)
+    full_positions = np.c_[xnew, ynew, z]
+
+    if path.opposite_direction:
+        full_positions = np.flip(full_positions, axis=0)
+
+    # [counts/s]
+    cps = calculate_count_rate_per_second(
+        landscape, full_positions, detector
+    )
+
+    # reshape so each segment is on a row
+    cps_per_seg = cps.reshape(num_segments, points_per_segment)
+
+    du = s[1] - s[0]
+    integrated_counts = np.trapezoid(cps_per_seg, dx=du, axis=1) / velocity
+    int_counts_result = np.zeros(num_points)
+    int_counts_result[1:] = integrated_counts
+
+    return original_distances, s, cps, int_counts_result

src/pg_rad/plotting/landscape_plotter.py

@@ -0,0 +1,150 @@
+import logging
+
+from matplotlib import pyplot as plt
+from matplotlib.axes import Axes
+from matplotlib.patches import Circle
+
+from numpy import median
+
+from pg_rad.landscape.landscape import Landscape
+
+logger = logging.getLogger(__name__)
+
+
+class LandscapeSlicePlotter:
+    def plot(
+        self,
+        landscape: Landscape,
+        z: int = 0,
+        show: bool = True,
+        save: bool = False,
+        ax: Axes | None = None
+    ):
+        """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.
+            show (bool, optional): Show the plot. Defaults to True.
+            save (bool, optional): Save the plot. Defaults to False.
+        """        """
+
+        """
+        self.z = z
+
+        if not ax:
+            fig, ax = plt.subplots()
+
+        self._draw_base(ax, landscape)
+        self._draw_path(ax, landscape)
+        self._draw_point_sources(ax, landscape)
+
+        ax.set_aspect("equal")
+
+        if save and not ax:
+            landscape_name = landscape.name.lower().replace(' ', '_')
+            filename = f"{landscape_name}_z{self.z}.png"
+            plt.savefig(filename)
+            logger.info("Plot saved to file: "+filename)
+
+        if show and not ax:
+            plt.show()
+
+        return ax
+
+    def _draw_base(self, ax, landscape: Landscape):
+        width, height = landscape.size[:2]
+
+        ax.set_xlim(right=max(width, .5*height))
+
+        # if the road is very flat, we center it vertically (looks better)
+        if median(landscape.path.y_list) == 0:
+            h = max(height, .5*width)
+            ax.set_ylim(bottom=-h//2,
+                        top=h//2)
+        else:
+            ax.set_ylim(top=max(height, .5*width))
+
+        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: Landscape):
+        if landscape.path.z <= self.z:
+            ax.plot(
+                landscape.path.x_list,
+                landscape.path.y_list,
+                linestyle='-',
+                marker='|',
+                markersize=3,
+                linewidth=1
+            )
+            if len(landscape.path.x_list) >= 2:
+                ax = self._draw_path_direction_arrow(ax, landscape.path)
+
+        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:
+            x, y, z = s.pos
+            if z <= self.z:
+                dot = Circle(
+                        (x, y),
+                        radius=5,
+                        color=s.color,
+                        zorder=5
+                )
+
+                ax.text(
+                    x + 0.06,
+                    y + 0.06,
+                    s.name+", z="+str(z),
+                    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."
+                    )
+
+    def _draw_path_direction_arrow(self, ax, path) -> Axes:
+        inset_ax = ax.inset_axes([0.8, 0.1, 0.15, 0.15])
+
+        x_start, y_start = path.x_list[0], path.y_list[0]
+        x_end, y_end = path.x_list[1], path.y_list[1]
+
+        dx = x_end - x_start
+        dy = y_end - y_start
+
+        if path.opposite_direction:
+            dx = -dx
+            dy = -dy
+
+        length = 10
+        dx_norm = dx / (dx**2 + dy**2)**0.5 * length
+        dy_norm = dy / (dx**2 + dy**2)**0.5 * length
+
+        inset_ax.arrow(
+            0, 0,
+            dx_norm, dy_norm,
+            head_width=5, head_length=5,
+            fc='red', ec='red',
+            zorder=4, linewidth=1
+        )
+
+        inset_ax.set_xlim(-2*length, 2*length)
+        inset_ax.set_ylim(-2*length, 2*length)
+        inset_ax.set_title("Direction", fontsize=8)
+        inset_ax.set_xticks([])
+        inset_ax.set_yticks([])
+        return ax

src/pg_rad/plotting/result_plotter.py

@@ -0,0 +1,206 @@
+from importlib.resources import files
+
+import numpy as np
+import pandas as pd
+from matplotlib import pyplot as plt
+from matplotlib.gridspec import GridSpec
+
+from .landscape_plotter import LandscapeSlicePlotter
+from pg_rad.simulator.outputs import SimulationOutput
+from pg_rad.landscape.landscape import Landscape
+
+
+class ResultPlotter:
+    def __init__(self, landscape: Landscape, output: SimulationOutput):
+        self.landscape = landscape
+        self.count_rate_res = output.count_rate
+        self.source_res = output.sources
+
+    def plot(self, landscape_z: float = 0):
+        self._plot_main(landscape_z)
+        self._plot_detector()
+        self._plot_metadata()
+
+        plt.show()
+
+    def _plot_main(self, landscape_z):
+        fig = plt.figure(figsize=(12, 8))
+        fig.suptitle(self.landscape.name)
+        fig.canvas.manager.set_window_title("Main results")
+
+        gs = GridSpec(
+            2, 2, figure=fig,
+            height_ratios=[1, 2], width_ratios=[1, 1],
+            hspace=0.3, wspace=0.3)
+
+        ax_cps = fig.add_subplot(gs[0, 0])
+        self._draw_cps(ax_cps)
+
+        ax_counts = fig.add_subplot(gs[0, 1])
+        self._draw_count_rate(ax_counts)
+
+        ax_landscape = fig.add_subplot(gs[1, :])
+        self._plot_landscape(ax_landscape, landscape_z)
+
+    def _plot_detector(self):
+        det = self.landscape.detector
+
+        fig = plt.figure(figsize=(10, 4))
+        fig.canvas.manager.set_window_title("Detector")
+
+        gs = GridSpec(1, 2, figure=fig, width_ratios=[0.5, 0.5])
+
+        ax_table = fig.add_subplot(gs[0, 0])
+        self._draw_detector_table(ax_table)
+
+        if not det.is_isotropic:
+            ax_polar = fig.add_subplot(gs[0, 1], projection='polar')
+
+            energies = [
+                source.primary_gamma for source in self.source_res
+            ]
+
+            self._draw_angular_efficiency_polar(ax_polar, det, energies[0])
+
+    def _plot_metadata(self):
+        fig, axs = plt.subplots(2, 1, figsize=(10, 6))
+        fig.canvas.manager.set_window_title("Simulation Metadata")
+
+        self._draw_table(axs[0])
+        self._draw_source_table(axs[1])
+
+    def _plot_landscape(self, ax, z):
+        lp = LandscapeSlicePlotter()
+        ax = lp.plot(landscape=self.landscape, z=z, ax=ax, show=False)
+        return ax
+
+    def _draw_cps(self, ax):
+        x = self.count_rate_res.sub_points
+        y = self.count_rate_res.cps
+        ax.plot(x, y, color='b')
+        ax.set_title('Counts per second (CPS)')
+        ax.set_xlabel('Arc length s [m]')
+        ax.set_ylabel('CPS [s$^{-1}$]')
+
+    def _draw_count_rate(self, ax):
+        x = self.count_rate_res.acquisition_points
+        y = self.count_rate_res.integrated_counts
+        ax.plot(x, y, color='r', linestyle='--', alpha=0.2)
+        ax.scatter(x, y, color='r', marker='x')
+        ax.set_title('Integrated counts')
+        ax.set_xlabel('Arc length s [m]')
+        ax.set_ylabel('N')
+
+    def _draw_table(self, ax):
+        ax.set_axis_off()
+        ax.set_title('Simulation parameters')
+        cols = ('Parameter', 'Value')
+        data = [
+            ["Air density (kg/m^3)", round(self.landscape.air_density, 3)],
+            ["Total path length (m)", round(self.landscape.path.length, 3)],
+            ["Readout points", len(self.count_rate_res.integrated_counts)],
+        ]
+
+        ax.table(
+            cellText=data,
+            colLabels=cols,
+            loc='center'
+        )
+
+        return ax
+
+    def _draw_detector_table(self, ax):
+        det = self.landscape.detector
+
+        if det.is_isotropic:
+            det_type = "Isotropic"
+        else:
+            det_type = "Angular"
+
+        source_energies = [
+            source.primary_gamma for source in self.source_res
+        ]
+
+        # list field efficiencies for each primary gamma in the landscape
+        effs = {e: det.get_efficiency(e) for e in source_energies}
+        formatted_effs = ", ".join(
+            f"{value:.3f} @ {key:.1f} keV"
+            for key, value in effs.items()
+        )
+        ax.set_axis_off()
+        ax.set_title('Detector')
+        cols = ('Parameter', 'Value')
+        data = [
+            ["Detector", f"{det.name} ({det_type})"],
+            ["Field efficiency", formatted_effs],
+        ]
+
+        ax.table(
+            cellText=data,
+            colLabels=cols,
+            loc='center'
+        )
+
+        return ax
+
+    def _draw_source_table(self, ax):
+        ax.set_axis_off()
+        ax.set_title('Point sources')
+
+        cols = (
+            'Name',
+            'Isotope',
+            'Activity (MBq)',
+            'Position (m)',
+            'Dist. to path (m)'
+        )
+
+        # this formats position to tuple
+        data = [
+            [
+                s.name,
+                s.isotope+f" ({s.primary_gamma} keV)",
+                s.activity,
+                "("+", ".join(f"{val:.2f}" for val in s.position)+")",
+                round(s.dist_from_path, 2)
+            ]
+            for s in self.source_res
+        ]
+        ax.table(
+            cellText=data,
+            colLabels=cols,
+            loc='center'
+        )
+
+        return ax
+
+    def _draw_angular_efficiency_polar(self, ax, detector, energy_keV):
+        # find the energies available for this detector.
+        csv = files('pg_rad.data.angular_efficiencies').joinpath(
+            detector.name + '.csv'
+        )
+        data = pd.read_csv(csv)
+        energy_cols = [col for col in data.columns if col != "angle"]
+        energies = np.array([float(col) for col in energy_cols])
+
+        # take energy column that is within 1% tolerance of energy_keV
+        rel_diff = np.abs(energies - energy_keV) / energies
+        match_idx = np.where(rel_diff <= 0.01)[0]
+        best_idx = match_idx[np.argmin(rel_diff[match_idx])]
+        col = energy_cols[best_idx]
+
+        theta_deg = data["angle"].to_numpy()
+        eff = data[col].to_numpy()
+
+        idx = np.argsort(theta_deg)
+        theta_deg = theta_deg[idx]
+        eff = eff[idx]
+
+        theta_deg = np.append(theta_deg, theta_deg[0])
+        eff = np.append(eff, eff[0])
+
+        theta_rad = np.radians(theta_deg)
+
+        print(theta_rad)
+        ax.plot(theta_rad, eff)
+        ax.set_title(f"Rel. angular efficiency @ {energy_keV:.1f} keV")

src/pg_rad/simulator/engine.py

@@ -0,0 +1,71 @@
+from typing import List
+
+from pg_rad.landscape.landscape import Landscape
+from pg_rad.simulator.outputs import (
+    CountRateOutput,
+    SimulationOutput,
+    SourceOutput
+)
+
+from pg_rad.physics.fluence import calculate_counts_along_path
+from pg_rad.utils.projection import minimal_distance_to_path
+from pg_rad.inputparser.specs import RuntimeSpec, SimulationOptionsSpec
+
+
+class SimulationEngine:
+    """Takes a fully built landscape and produces results."""
+    def __init__(
+        self,
+        landscape: Landscape,
+        runtime_spec: RuntimeSpec,
+        sim_spec: SimulationOptionsSpec,
+    ):
+
+        self.landscape = landscape
+        self.detector = self.landscape.detector
+        self.runtime_spec = runtime_spec
+        self.sim_spec = sim_spec
+
+    def simulate(self) -> SimulationOutput:
+        """Compute everything and return structured output."""
+
+        count_rate_results = self._calculate_count_rate_along_path()
+        source_results = self._calculate_point_source_distance_to_path()
+
+        return SimulationOutput(
+            name=self.landscape.name,
+            count_rate=count_rate_results,
+            sources=source_results
+        )
+
+    def _calculate_count_rate_along_path(self) -> CountRateOutput:
+        acq_points, sub_points, cps, int_counts = calculate_counts_along_path(
+            self.landscape,
+            self.detector,
+            velocity=self.runtime_spec.speed
+        )
+
+        return CountRateOutput(acq_points, sub_points, cps, int_counts)
+
+    def _calculate_point_source_distance_to_path(self) -> List[SourceOutput]:
+
+        path = self.landscape.path
+        source_output = []
+        for s in self.landscape.point_sources:
+            dist_to_path = minimal_distance_to_path(
+                path.x_list,
+                path.y_list,
+                path.z,
+                s.pos)
+
+            source_output.append(
+                SourceOutput(
+                    s.name,
+                    s.isotope.name,
+                    s.isotope.E,
+                    s.activity,
+                    s.pos,
+                    dist_to_path)
+            )
+
+        return source_output

src/pg_rad/simulator/outputs.py

@@ -0,0 +1,28 @@
+from typing import List, Tuple
+
+from dataclasses import dataclass
+
+
+@dataclass
+class CountRateOutput:
+    acquisition_points: List[float]
+    sub_points: List[float]
+    cps: List[float]
+    integrated_counts: List[float]
+
+
+@dataclass
+class SourceOutput:
+    name: str
+    isotope: str
+    primary_gamma: float
+    activity: float
+    position: Tuple[float, float, float]
+    dist_from_path: float
+
+
+@dataclass
+class SimulationOutput:
+    name: str
+    count_rate: CountRateOutput
+    sources: List[SourceOutput]

src/pg_rad/utils/interpolators.py

@@ -0,0 +1,56 @@
+from importlib.resources import files
+
+import numpy as np
+from pandas import read_csv
+from scipy.interpolate import interp1d, CubicSpline
+
+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)
+
+
+def get_field_efficiency(name: str, energy_keV: float) -> float:
+    csv = files('pg_rad.data.field_efficiencies').joinpath(name+'.csv')
+    data = read_csv(csv)
+    data = data.groupby("energy_keV", as_index=False).mean()
+    x = data["energy_keV"].to_numpy()
+    y = data["field_efficiency_m2"].to_numpy()
+    f = CubicSpline(x, y)
+    return f(energy_keV)
+
+
+def get_angular_efficiency(name: str, energy_keV: float, *angle: float):
+    csv = files('pg_rad.data.angular_efficiencies').joinpath(name+'.csv')
+    data = read_csv(csv)
+
+    # check all energies at which angular eff. is available for this detector.
+    # this is done within 1% tolerance
+    energy_cols = [col for col in data.columns if col != "angle"]
+    energies = np.array([float(col) for col in energy_cols])
+    rel_diff = np.abs(energies - energy_keV) / energies
+    match_idx = np.where(rel_diff <= 0.01)[0]
+
+    if len(match_idx) == 0:
+        raise NotImplementedError(
+            f"No angular efficiency defined for {energy_keV} keV "
+            f"in detector '{name}'. Available: {energies}"
+        )
+
+    best_idx = match_idx[np.argmin(rel_diff[match_idx])]
+    selected_energy_col = energy_cols[best_idx]
+
+    x = data["angle"].to_numpy()
+    y = data[selected_energy_col].to_numpy()
+    idx = np.argsort(x)
+    x = x[idx]
+    y = y[idx]
+    f = interp1d(x, y)
+
+    return f(angle)

src/pg_rad/utils/projection.py

@@ -0,0 +1,168 @@
+from typing import List, Tuple
+
+import numpy as np
+
+from pg_rad.exceptions.exceptions import OutOfBoundsError
+
+
+def rel_to_abs_source_position(
+    x_list: List,
+    y_list: List,
+    path_z: int | float,
+    along_path: int | float,
+    side: str,
+    dist_from_path: int | float,
+    z_rel: int | float = 0
+):
+    path = np.column_stack((x_list, y_list))
+    segments = np.diff(path, axis=0)
+    segment_lengths = np.hypot(segments[:, 0], segments[:, 1])
+
+    arc_length = np.cumsum(segment_lengths)
+    arc_length = np.insert(arc_length, 0, 0)
+
+    # find index of the segment where the source should be placed
+    s_i = np.searchsorted(arc_length, along_path, side='right') - 1
+    if not 0 <= s_i < len(segments):
+        raise OutOfBoundsError(
+            "along_path must be less than the length of the path!"
+            )
+
+    # fraction of segment length where to place the point source
+    segment_fraction = (
+        (along_path - arc_length[s_i])
+        / segment_lengths[s_i]
+    )
+
+    # tangential vector at point where point source to be placed
+    point_on_path = path[s_i] + segment_fraction * segments[s_i]
+    tangent_vec = segments[s_i] / segment_lengths[s_i]
+
+    if side not in {"left", "right"}:
+        raise ValueError("side must be 'left' or 'right'")
+
+    orientation = 1 if side == "left" else -1
+    perp_vec = orientation * np.array([-tangent_vec[1], tangent_vec[0]])
+
+    x_abs, y_abs = point_on_path + dist_from_path * perp_vec
+    z_abs = path_z + z_rel
+
+    return x_abs, y_abs, z_abs
+
+
+def minimal_distance_to_path(
+    x_list: List[float],
+    y_list: List[float],
+    path_z: float,
+    pos: Tuple[float, float, float]
+) -> float:
+    """
+    Compute minimal Euclidean distance from pos (x,y,z)
+    to interpolated polyline path at height path_z.
+    """
+    x, y, z = pos
+
+    path = np.column_stack((x_list, y_list))
+    point_xy = np.array([x, y])
+
+    segments = np.diff(path, axis=0)
+    segment_lengths = np.hypot(segments[:, 0], segments[:, 1])
+
+    min_dist = np.inf
+
+    for p0, seg, seg_len in zip(path[:-1], segments, segment_lengths):
+        if seg_len == 0:
+            continue
+
+        t = np.dot(point_xy - p0, seg) / (seg_len ** 2)
+        t = np.clip(t, 0.0, 1.0)
+
+        projection_xy = p0 + t * seg
+
+        dx, dy = point_xy - projection_xy
+        dz = z - path_z
+
+        dist = np.sqrt(dx**2 + dy**2 + dz**2)
+
+        if dist < min_dist:
+            min_dist = dist
+
+    return min_dist
+
+
+def abs_to_rel_source_position(
+    x_list: List,
+    y_list: List,
+    path_z: int | float,
+    pos: Tuple[float, float, float],
+) -> Tuple[float, float, str, float]:
+    """
+    Convert absolute (x,y,z) position into:
+        along_path,
+        dist_from_path,
+        side ("left" or "right"),
+        z_rel
+    """
+
+    x, y, z = pos
+
+    path = np.column_stack((x_list, y_list))
+    point = np.array([x, y])
+
+    segments = np.diff(path, axis=0)
+    segment_lengths = np.hypot(segments[:, 0], segments[:, 1])
+
+    arc_length = np.cumsum(segment_lengths)
+    arc_length = np.insert(arc_length, 0, 0)
+
+    min_dist = np.inf
+    best_projection = None
+    best_s_i = None
+    best_fraction = None
+    best_signed_dist = None
+
+    for (i, (p0, seg, seg_len)) in enumerate(
+        zip(path[:-1], segments, segment_lengths)
+    ):
+        if seg_len == 0:
+            continue
+
+        # project point onto infinite line
+        t = np.dot(point - p0, seg) / (seg_len ** 2)
+
+        # clamp to segment
+        t_clamped = np.clip(t, 0.0, 1.0)
+
+        projection = p0 + t_clamped * seg
+        diff_vec = point - projection
+        dist = np.linalg.norm(diff_vec)
+
+        if dist < min_dist:
+            min_dist = dist
+            best_projection = projection
+            best_s_i = i
+            best_fraction = t_clamped
+
+            # tangent and perpendicular
+            tangent_vec = seg / seg_len
+            perp_vec = np.array([-tangent_vec[1], tangent_vec[0]])
+
+            signed_dist = np.dot(diff_vec, perp_vec)
+            best_signed_dist = signed_dist
+
+    if best_projection is None:
+        raise ValueError("Could not project point onto path.")
+
+    # Compute along_path
+    along_path = (
+        arc_length[best_s_i] + best_fraction * segment_lengths[best_s_i]
+    )
+
+    # Side and distance
+    side = "left" if best_signed_dist > 0 else "right"
+    dist_from_path = abs(best_signed_dist)
+
+    # z relative
+    z_rel = z - path_z
+
+    return along_path, dist_from_path, side, z_rel

src/road_gen/generators/base_road_generator.py

@@ -0,0 +1,49 @@
+import secrets
+
+import numpy as np
+
+
+class BaseRoadGenerator:
+    """A base generator object for generating a road of a specified length."""
+    def __init__(
+            self,
+            ds: int | float,
+            velocity: int | float,
+            mu: float = 0.7,
+            g: float = 9.81,
+            seed: int | None = None
+            ):
+        """Initialize a BaseGenerator with a given or random seed.
+
+        Args:
+            ds (int | float): The step size in meters.
+            velocity (int | float): Velocity in meters per second.
+            mu (float): Coefficient of friction. Defaults to 0.7 (dry asphalt).
+            g (float): Acceleration due to gravity (m/s^2). Defaults to 9.81.
+            seed (int | None, optional): Set a seed for the generator.
+                Defaults to a random seed.
+        """
+        if seed is None:
+            seed = secrets.randbits(32)
+
+        if not isinstance(seed, int):
+            raise TypeError("seed must be an integer or None.")
+
+        if not isinstance(ds, int | float):
+            raise TypeError("Step size must be integer or float in meters.")
+
+        if not isinstance(velocity, int | float):
+            raise TypeError(
+                "Velocity must be integer or float in meters per second."
+            )
+
+        self.ds = ds
+
+        self.velocity = velocity
+        self.min_radius = (velocity ** 2) / (g * mu)
+
+        self.seed = seed
+        self._rng = np.random.default_rng(seed)
+
+    def generate(self):
+        pass

src/road_gen/generators/segmented_road_generator.py

@@ -0,0 +1,144 @@
+import logging
+from typing import Tuple
+
+import numpy as np
+
+from .base_road_generator import BaseRoadGenerator
+from road_gen.prefabs import prefabs
+from road_gen.integrator.integrator import integrate_road
+
+from pg_rad.configs import defaults
+
+
+logger = logging.getLogger(__name__)
+
+
+class SegmentedRoadGenerator(BaseRoadGenerator):
+    def __init__(
+        self,
+        ds: int | float,
+        velocity: int | float,
+        mu: float = 0.7,
+        g: float = 9.81,
+        seed: int | None = None
+    ):
+        """Initialize a SegmentedRoadGenerator with given or random seed.
+
+        Args:
+            ds (int | float): The step size in meters.
+            velocity (int | float): Velocity in meters per second.
+            mu (float): Coefficient of friction. Defaults to 0.7 (dry asphalt).
+            g (float): Acceleration due to gravity (m/s^2). Defaults to 9.81.
+            seed (int | None, optional): Set a seed for the generator.
+                Defaults to random seed.
+        """
+        super().__init__(ds, velocity, mu, g, seed)
+
+    def generate(
+            self,
+            segments: list[str],
+            lengths: list[int | float],
+            angles: list[float | None],
+            alpha: float = defaults.DEFAULT_ALPHA,
+            min_turn_angle: float = defaults.DEFAULT_MIN_TURN_ANGLE,
+            max_turn_angle: float = defaults.DEFAULT_MAX_TURN_ANGLE
+    ) -> Tuple[np.ndarray, np.ndarray]:
+        """Generate a curvature profile from a list of segments.
+
+        Args:
+            segments (list[str]): List of segments.
+            lengths (list[int | float]): List of segment lengths.
+            angles (list[float | None]): List of angles.
+            alpha (float, optional): Dirichlet concentration parameter.
+                A higher value leads to more uniform apportionment of the
+                length amongst the segments, while a lower value allows more
+                random apportionment. Defaults to 1.0.
+            min_turn_angle (float, optional): Minimum turn angle in degrees for
+                random sampling of turn radius. Does nothing if `angle_list` is
+                provided or no `turn_*` segement is specified in the `segments`
+                list.
+            min_turn_angle (float, optional): Maximum turn angle in degrees for
+                random sampling of turn radius. Does nothing if `angle_list` is
+                provided or no `turn_*` segement is specified in the `segments`
+                list.
+        Raises:
+            ValueError: Raised when a turn
+                is too tight given its segment length and the velocity.
+                To fix this, you can try to reduce the amount of segments or
+                increase length. Increasing alpha
+                (Dirichlet concentration parameter) can also help because this
+                reduces the odds of very small lengths being assigned to
+                turn segments.
+
+        Returns:
+            Tuple[np.ndarray, np.ndarray]: x and y coordinates of the
+                waypoints describing the random road.
+        """
+        existing_prefabs = prefabs.PREFABS.keys()
+        if not all(segment in existing_prefabs for segment in segments):
+            raise ValueError(
+                "Invalid segment type provided. Available choices"
+                f"{existing_prefabs}"
+            )
+
+        self.segments = segments
+        self.alpha = alpha
+
+        total_length = sum(lengths)
+        num_points = np.ceil(total_length / self.ds).astype(int)
+
+        # divide num_points into len(segments) randomly sized parts.
+        if len(lengths) == len(segments):
+            parts = np.array([seg_len / total_length for seg_len in lengths])
+        else:
+            parts = self._rng.dirichlet(
+                np.full(len(segments), alpha),
+                size=1)[0]
+
+        parts = parts * num_points
+        parts = np.round(parts).astype(int)
+
+        # correct round off so the sum of parts is still total length L.
+        if sum(parts) != num_points:
+            parts[0] += num_points - sum(parts)
+
+        curvature = np.zeros(num_points)
+        current_index = 0
+
+        for seg_name, seg_length, seg_angle in zip(segments, parts, angles):
+            seg_function = prefabs.PREFABS[seg_name]
+
+            if seg_name == 'straight':
+                curvature_s = seg_function(seg_length)
+            else:
+                R_min_angle = seg_length / np.deg2rad(max_turn_angle)
+                R_max_angle = seg_length / np.deg2rad(min_turn_angle)
+
+                # physics limit
+                R_min = max(self.min_radius, R_min_angle)
+
+                if R_min > R_max_angle:
+                    raise ValueError(
+                        f"{seg_name} with length {seg_length} does not have "
+                        "a possible radius. The minimum for the provided "
+                        "velocity and friction coefficient is "
+                        f"{self.min_radius}, but the possible range is "
+                        f"({R_min}, {R_max_angle})"
+                    )
+
+                if seg_angle:
+                    radius = seg_length / np.deg2rad(seg_angle)
+                else:
+                    radius = self._rng.uniform(R_min, R_max_angle)
+
+                if seg_name.startswith("u_turn"):
+                    curvature_s = seg_function(radius)
+                else:
+                    curvature_s = seg_function(seg_length, radius)
+
+            curvature[current_index:(current_index + seg_length)] = curvature_s
+            current_index += seg_length
+
+            x, y = integrate_road(curvature)
+
+        return x * self.ds, y * self.ds

src/road_gen/integrator/integrator.py

@@ -0,0 +1,21 @@
+import numpy as np
+
+
+def integrate_road(curvature: np.ndarray, ds: float = 1.0):
+    """Integrate along the curvature field to obtain X and Y coordinates.
+
+    Args:
+        curvature (np.ndarray): The curvature field.
+        ds (float, optional): _description_. Defaults to 1.0.
+
+    Returns:
+        x, y (np.ndarray): X and Y waypoints.
+    """
+
+    theta = np.zeros(len(curvature))
+    theta[1:] = np.cumsum(curvature[:-1] * ds)
+
+    x = np.cumsum(np.cos(theta) * ds)
+    y = np.cumsum(np.sin(theta) * ds)
+
+    return x, y

src/road_gen/prefabs/prefabs.py

@@ -0,0 +1,20 @@
+import numpy as np
+
+
+def straight(length: int) -> np.ndarray:
+    return np.zeros(length)
+
+
+def turn_left(length: int, radius: float) -> np.ndarray:
+    return np.full(length, 1.0 / radius)
+
+
+def turn_right(length: int, radius: float) -> np.ndarray:
+    return -turn_left(length, radius)
+
+
+PREFABS = {
+    'straight': straight,
+    'turn_left': turn_left,
+    'turn_right': turn_right,
+}

tests/test_attenuation_functions.py

@@ -0,0 +1,29 @@
+import pytest
+
+from pg_rad.utils.interpolators 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

tests/test_fluence_rate.py

@@ -0,0 +1,70 @@
+import numpy as np
+import pytest
+
+from pg_rad.inputparser.parser import ConfigParser
+from pg_rad.landscape.director import LandscapeDirector
+from pg_rad.physics.fluence import phi
+
+
+@pytest.fixture
+def isotropic_detector():
+    from pg_rad.detector.detector import load_detector
+    return load_detector('dummy')
+
+
+@pytest.fixture
+def phi_ref(test_landscape, isotropic_detector):
+    source = test_landscape.point_sources[0]
+
+    r = np.linalg.norm(np.array([0, 0, 0]) - np.array(source.pos))
+
+    A = source.activity * 1E6
+    b = source.isotope.b
+    mu_air = source.isotope.mu_mass_air * test_landscape.air_density
+    mu_air *= 0.1
+
+    eff = isotropic_detector.get_efficiency(source.isotope.E)
+    return A * eff * b * np.exp(-mu_air * r) / (4 * np.pi * r**2)
+
+
+@pytest.fixture
+def test_landscape():
+
+    test_yaml = """
+    name: Test landscape
+    speed: 8.33
+    acquisition_time: 1
+
+    path:
+        length: 1000
+        segments:
+            - straight
+
+    sources:
+        test_source:
+            activity_MBq: 100
+            position: [0, 100, 0]
+            isotope: Cs137
+            gamma_energy_keV: 661
+
+    detector: dummy
+    """
+
+    cp = ConfigParser(test_yaml).parse()
+    landscape = LandscapeDirector.build_from_config(cp)
+    return landscape
+
+
+def test_single_source_fluence(phi_ref, test_landscape, isotropic_detector):
+    s = test_landscape.point_sources[0]
+    r = np.linalg.norm(np.array([0, 0, 0]) - np.array(s.pos))
+    phi_calc = phi(
+        r,
+        s.activity*1E6,
+        s.isotope.b,
+        s.isotope.mu_mass_air,
+        test_landscape.air_density,
+        isotropic_detector.get_efficiency(s.isotope.E)
+        )
+
+    assert pytest.approx(phi_calc, rel=1E-6) == phi_ref

tests/test_path_functionality.py

@@ -1,47 +0,0 @@
-import pathlib
-
-import numpy as np
-import pytest
-
-from pg_rad.dataloader import load_data
-from pg_rad.path import Path, path_from_RT90
-
-@pytest.fixture
-def test_df():
-    csv_path = pathlib.Path(__file__).parent / "data/coordinates.csv"
-    return load_data(csv_path)
-
-def test_piecewise_regression(test_df):
-    """_Verify whether the intermediate points deviate less than 0.1 SD._"""
-        
-    p_full = path_from_RT90(
-        test_df,
-        east_col="East",
-        north_col="North",
-        simplify_path=False
-    )
- 
-    p_simpl = path_from_RT90(
-        test_df,
-        east_col="East",
-        north_col="North",
-        simplify_path=True
-    )
-
-    x_f = np.array(p_full.x_list)
-    y_f = np.array(p_full.y_list)
-
-    x_s = np.array(p_simpl.x_list)
-    y_s = np.array(p_simpl.y_list)
-
-    sd = np.std(y_f)
-
-    for xb, yb in zip(x_s[1:-1], y_s[1:-1]):
-        # find nearest original x index
-        idx = np.argmin(np.abs(x_f - xb))
-        deviation = abs(yb - y_f[idx])
-
-        assert deviation < 0.1 * sd, (
-            f"Breakpoint deviation too large: {deviation:.4f} "
-            f"(threshold {0.1 * sd:.4f}) at x={xb:.2f}"
-            )

tests/test_roi.py

@@ -0,0 +1,42 @@
+import pytest
+
+from pg_rad.background.background import get_roi_from_fwhm, get_cps_from_roi
+from pg_rad.detector.detector import load_detector
+
+E_gamma = 661.7  # keV
+
+
+@pytest.fixture
+def detector():
+    return load_detector("LU_NaI_3inch")
+
+
+@pytest.fixture
+def roi_Cs137_NaI(detector):
+    return get_roi_from_fwhm(detector, E_gamma)
+
+
+@pytest.mark.parametrize("ref_roi_lo, ref_roi_hi", [
+    (537.58, 792.69)
+])
+def test_roi_acquisition_Cs137_NaI(ref_roi_lo, ref_roi_hi, roi_Cs137_NaI):
+    """Test against ROI used in GammaVision"""
+    est_roi_lo, est_roi_hi = roi_Cs137_NaI
+
+    assert pytest.approx(est_roi_lo, rel=0.05) == ref_roi_lo
+    assert pytest.approx(est_roi_hi, rel=0.05) == ref_roi_hi
+
+
+@pytest.mark.parametrize("ref_cps", [
+    13.61
+])
+def test_cps_acquisition_Cs137_NaI(
+    ref_cps,
+    roi_Cs137_NaI,
+    detector
+):
+    """Test against cps from GammaVision"""
+    est_roi_lo, est_roi_hi = roi_Cs137_NaI
+    est_cps = get_cps_from_roi(detector, est_roi_lo, est_roi_hi)
+
+    assert pytest.approx(est_cps, rel=0.1) == ref_cps

tests/test_sources.py

@@ -1,27 +1,35 @@
 import numpy as np
 import pytest
 
-from pg_rad.sources import PointSource
+from pg_rad.objects.sources import PointSource
+
 
 @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(position=pos_a, activity_MBq=None, isotope=iso)
+    b = PointSource(position=pos_b, activity_MBq=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