From 431ede2ebe7532c9706bef4e75c7607b56ea224d Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Wed, 14 Jun 2023 16:23:42 -0300
Subject: [PATCH 01/17] Allowing customization of the center frequency and of
 material layers

---
 .../plot_customized_center_frequency.py       | 27 +++++++++++
 docs/examples/plot_full_scenario.py           | 39 +++++++--------
 src/neurotechdevkit/__init__.py               |  3 +-
 src/neurotechdevkit/materials.py              | 38 +++++++++++++++
 src/neurotechdevkit/scenarios/__init__.py     |  2 +-
 src/neurotechdevkit/scenarios/_base.py        | 16 ++-----
 src/neurotechdevkit/scenarios/_scenario_0.py  | 18 +++----
 src/neurotechdevkit/scenarios/_scenario_1.py  | 19 ++++----
 src/neurotechdevkit/scenarios/_scenario_2.py  | 16 +++----
 src/neurotechdevkit/scenarios/materials.py    | 48 -------------------
 tests/neurotechdevkit/scenarios/test_base.py  |  6 +--
 11 files changed, 111 insertions(+), 121 deletions(-)
 create mode 100644 docs/examples/plot_customized_center_frequency.py
 create mode 100644 src/neurotechdevkit/materials.py
 delete mode 100644 src/neurotechdevkit/scenarios/materials.py

diff --git a/docs/examples/plot_customized_center_frequency.py b/docs/examples/plot_customized_center_frequency.py
new file mode 100644
index 00000000..4b16e4ba
--- /dev/null
+++ b/docs/examples/plot_customized_center_frequency.py
@@ -0,0 +1,27 @@
+# -*- coding: utf-8 -*-
+"""
+Plot customized center frequency
+====================================
+
+This example demonstrates how to use a customized center frequency
+using ndk
+"""
+# %%
+import neurotechdevkit as ndk
+
+CENTER_FREQUENCY = 6e5
+
+scenario = ndk.make('scenario-0-v0')
+
+# Customizing material layers with random values
+scenario._material_layers = [
+  ("water", ndk.materials.Material(1600.0, 1100.0, 0.0, "#2E86AB")),
+  ("skull", ndk.materials.Material(3000.0, 1850.0, 4.0, "#FAF0CA")),
+  ("brain", ndk.materials.Material(1660.0, 1140.0, 0.3, "#DB504A")),
+  ("tumor", ndk.materials.Material(1850.0, 1250.0, 0.8, "#94332F")),
+]
+
+result = scenario.simulate_steady_state(center_frequency=CENTER_FREQUENCY)
+result.render_steady_state_amplitudes(show_material_outlines=False)
+
+# %%
diff --git a/docs/examples/plot_full_scenario.py b/docs/examples/plot_full_scenario.py
index 7c026a79..2774bd1c 100644
--- a/docs/examples/plot_full_scenario.py
+++ b/docs/examples/plot_full_scenario.py
@@ -41,6 +41,24 @@ class FullScenario(Scenario2D):
         "target_1": Target("target_1", np.array([0.064, 0.0]), 0.004, ""),
     }
 
+    """
+    The order of returned materials defines the layering of the scenario.
+
+    Each material has the following properties:
+    - vp: the speed of sound (in m/s).
+    - rho: the mass density (in kg/m³).
+    - alpha: the absorption (in dB/cm).
+    - render_color: the color used when rendering this material in the scenario layout
+        plot.
+    """
+    _material_layers = [
+        ("water", materials.water),
+        ("skin", materials.skin),
+        ("cortical bone", materials.cortical_bone),
+        ("trabecular bone", materials.trabecular_bone),
+        ("brain", materials.brain),
+    ]
+
     def __init__(self, complexity="fast"):
         """
         Instantiate a new scenario.
@@ -53,27 +71,6 @@ def __init__(self, complexity="fast"):
             complexity=complexity,
         )
 
-    @property
-    def _material_layers(self) -> list[tuple[str, Struct]]:
-        """"
-        This function defines the material layers for the scenario.
-        The order of returned materials defines the layering of the scenario.
-
-        Each material has the following properties:
-        - vp: the speed of sound (in m/s).
-        - rho: the mass density (in kg/m³).
-        - alpha: the absorption (in dB/cm).
-        - render_color: the color used when rendering this material in the scenario layout
-            plot.
-        """
-        return [
-            ("water", materials.water),
-            ("skin", materials.skin),
-            ("cortical bone", materials.cortical_bone),
-            ("trabecular bone", materials.trabecular_bone),
-            ("brain", materials.brain),
-        ]
-
     @property
     def _material_outline_upsample_factor(self) -> int:
         """
diff --git a/src/neurotechdevkit/__init__.py b/src/neurotechdevkit/__init__.py
index 3cefd11c..14b97161 100644
--- a/src/neurotechdevkit/__init__.py
+++ b/src/neurotechdevkit/__init__.py
@@ -3,12 +3,13 @@
 
 import os
 
-from . import scenarios, sources
+from . import materials, scenarios, sources
 from .results import load_result_from_disk
 
 __all__ = [
     "results",
     "scenarios",
+    "materials",
     "sources",
     "make",
     "ScenarioNotFoundError",
diff --git a/src/neurotechdevkit/materials.py b/src/neurotechdevkit/materials.py
new file mode 100644
index 00000000..b7db740c
--- /dev/null
+++ b/src/neurotechdevkit/materials.py
@@ -0,0 +1,38 @@
+"""Materials for the neurotechdevkit scenarios."""
+from dataclasses import dataclass
+
+from mosaic.types import Struct
+
+
+@dataclass
+class Material:
+    """A material with properties used in the neurotechdevkit scenarios."""
+
+    vp: float
+    rho: float
+    alpha: float
+    render_color: str
+
+    def to_struct(self) -> Struct:
+        """Return a Struct representation of the material.
+
+        Returns:
+            Struct: a Struct representation of the material.
+        """
+        struct = Struct()
+        struct.vp = self.vp
+        struct.rho = self.rho
+        struct.alpha = self.alpha
+        struct.render_color = self.render_color
+        return struct
+
+
+water = Material(1500.0, 1000.0, 0.0, "#2E86AB")
+skin = Material(1610.0, 1090.0, 0.2, "#FA8B53")
+cortical_bone = Material(2800.0, 1850.0, 4.0, "#FAF0CA")
+trabecular_bone = Material(2300.0, 1700.0, 8.0, "#EBD378")
+brain = Material(1560.0, 1040.0, 0.3, "#DB504A")
+
+# these numbers are completely made up
+# TODO: research reasonable values
+tumor = Material(1650.0, 1150.0, 0.8, "#94332F")
diff --git a/src/neurotechdevkit/scenarios/__init__.py b/src/neurotechdevkit/scenarios/__init__.py
index dc4d0181..d10110e5 100644
--- a/src/neurotechdevkit/scenarios/__init__.py
+++ b/src/neurotechdevkit/scenarios/__init__.py
@@ -1,5 +1,5 @@
 """Scenarios module."""
-from . import materials
+from .. import materials
 from ._base import Scenario, Scenario2D, Scenario3D, Target
 from ._scenario_0 import Scenario0
 from ._scenario_1 import Scenario1_2D, Scenario1_3D
diff --git a/src/neurotechdevkit/scenarios/_base.py b/src/neurotechdevkit/scenarios/_base.py
index 2e0c28e0..bb02ceac 100644
--- a/src/neurotechdevkit/scenarios/_base.py
+++ b/src/neurotechdevkit/scenarios/_base.py
@@ -16,6 +16,7 @@
 from stride.problem import StructuredData
 
 from .. import rendering, results
+from ..materials import Material
 from ..sources import Source
 from ._resources import budget_time_and_memory_resources
 from ._shots import create_shot
@@ -59,6 +60,7 @@ class Scenario(abc.ABC):
 
     _SCENARIO_ID: str
     _TARGET_OPTIONS: dict[str, Target]
+    _material_layers: list[tuple[str, Material]]
 
     def __init__(
         self,
@@ -224,7 +226,7 @@ def materials(self) -> Mapping[str, Struct]:
         - render_color: the color used when rendering this material in the
         scenario layout plot.
         """
-        return {name: material for name, material in self._material_layers}
+        return {name: material.to_struct() for name, material in self._material_layers}
 
     @property
     def layer_ids(self) -> Mapping[str, int]:
@@ -236,11 +238,6 @@ def ordered_layers(self) -> list[str]:
         """A list of material names in order of their layer id."""
         return [name for name, _ in self._material_layers]
 
-    @property
-    @abc.abstractmethod
-    def _material_layers(self) -> list[tuple[str, Struct]]:
-        pass
-
     @property
     @abc.abstractmethod
     def _material_outline_upsample_factor(self) -> int:
@@ -406,13 +403,6 @@ def simulate_steady_state(
         Returns:
             An object containing the result of the steady-state simulation.
         """
-        if center_frequency != 5.0e5:
-            raise NotImplementedError(
-                "500kHz is the only currently supported center frequency. Support for"
-                " other frequencies will be implemented once material properties as a"
-                " function of frequency has been implemented."
-            )
-
         problem = self.problem
         sim_time = select_simulation_time_for_steady_state(
             grid=problem.grid,
diff --git a/src/neurotechdevkit/scenarios/_scenario_0.py b/src/neurotechdevkit/scenarios/_scenario_0.py
index 2c631ca4..29072edb 100644
--- a/src/neurotechdevkit/scenarios/_scenario_0.py
+++ b/src/neurotechdevkit/scenarios/_scenario_0.py
@@ -1,9 +1,8 @@
 import numpy as np
 import stride
-from mosaic.types import Struct
 
+from .. import materials
 from ..sources import FocusedSource2D
-from . import materials
 from ._base import Scenario2D, Target
 from ._utils import (
     add_material_fields_to_problem,
@@ -25,6 +24,12 @@ class Scenario0(Scenario2D):
             description="Represents a simulated tumor.",
         ),
     }
+    _material_layers = [
+        ("water", materials.water),
+        ("skull", materials.cortical_bone),
+        ("brain", materials.brain),
+        ("tumor", materials.tumor),
+    ]
 
     def __init__(self, complexity="fast"):
         """Create a new instance of scenario 0."""
@@ -35,15 +40,6 @@ def __init__(self, complexity="fast"):
             complexity=complexity,
         )
 
-    @property
-    def _material_layers(self) -> list[tuple[str, Struct]]:
-        return [
-            ("water", materials.water),
-            ("skull", materials.cortical_bone),
-            ("brain", materials.brain),
-            ("tumor", materials.tumor),
-        ]
-
     @property
     def _material_outline_upsample_factor(self) -> int:
         return 16
diff --git a/src/neurotechdevkit/scenarios/_scenario_1.py b/src/neurotechdevkit/scenarios/_scenario_1.py
index c728d797..7022744b 100644
--- a/src/neurotechdevkit/scenarios/_scenario_1.py
+++ b/src/neurotechdevkit/scenarios/_scenario_1.py
@@ -7,8 +7,7 @@
 import stride
 from mosaic.types import Struct
 
-from .. import rendering, sources
-from . import materials
+from .. import materials, rendering, sources
 from ._base import Scenario2D, Scenario3D, Target
 from ._utils import add_material_fields_to_problem, make_grid
 
@@ -51,15 +50,13 @@ class _Scenario1Mixin:
         https://asa.scitation.org/doi/pdf/10.1121/10.0013426
     """
 
-    @property
-    def _material_layers(self: _Scenario1MixinProtocol) -> list[tuple[str, Struct]]:
-        return [
-            ("water", materials.water),
-            ("skin", materials.skin),
-            ("cortical bone", materials.cortical_bone),
-            ("trabecular bone", materials.trabecular_bone),
-            ("brain", materials.brain),
-        ]
+    _material_layers = [
+        ("water", materials.water),
+        ("skin", materials.skin),
+        ("cortical bone", materials.cortical_bone),
+        ("trabecular bone", materials.trabecular_bone),
+        ("brain", materials.brain),
+    ]
 
     @property
     def _material_outline_upsample_factor(self) -> int:
diff --git a/src/neurotechdevkit/scenarios/_scenario_2.py b/src/neurotechdevkit/scenarios/_scenario_2.py
index a0ee1406..4e9a0332 100644
--- a/src/neurotechdevkit/scenarios/_scenario_2.py
+++ b/src/neurotechdevkit/scenarios/_scenario_2.py
@@ -9,8 +9,7 @@
 import stride
 from mosaic.types import Struct
 
-from .. import rendering, sources
-from . import materials
+from .. import materials, rendering, sources
 from ._base import Scenario2D, Scenario3D, Target
 from ._utils import add_material_fields_to_problem, make_grid
 
@@ -51,14 +50,6 @@ class _Scenario2Mixin:
         https://asa.scitation.org/doi/pdf/10.1121/10.0013426
     """
 
-    @property
-    def _material_layers(self: _Scenario2MixinProtocol) -> list[tuple[str, Struct]]:
-        return [
-            ("water", materials.water),
-            ("skull", materials.cortical_bone),
-            ("brain", materials.brain),
-        ]
-
     @property
     def _material_outline_upsample_factor(self) -> int:
         return 4
@@ -138,6 +129,11 @@ class Scenario2_2D(_Scenario2Mixin, Scenario2D):
             ),
         ),
     }
+    _material_layers = [
+        ("water", materials.water),
+        ("skull", materials.cortical_bone),
+        ("brain", materials.brain),
+    ]
 
     def __init__(self, complexity="fast"):
         """Create a new instance of scenario 2."""
diff --git a/src/neurotechdevkit/scenarios/materials.py b/src/neurotechdevkit/scenarios/materials.py
deleted file mode 100644
index 78222f25..00000000
--- a/src/neurotechdevkit/scenarios/materials.py
+++ /dev/null
@@ -1,48 +0,0 @@
-"""Materials for the neurotechdevkit scenarios."""
-from mosaic.types import Struct
-
-# TODO: encapsulate mosaic struct behind an NDK materials type
-
-
-water = Struct()
-water.vp = 1500.0
-water.rho = 1000.0
-water.alpha = 0.0
-water.render_color = "#2E86AB"
-
-
-skin = Struct()
-skin.vp = 1610.0
-skin.rho = 1090.0
-skin.alpha = 0.2
-skin.render_color = "#FA8B53"
-
-
-cortical_bone = Struct()
-cortical_bone.vp = 2800.0
-cortical_bone.rho = 1850.0
-cortical_bone.alpha = 4.0
-cortical_bone.render_color = "#FAF0CA"
-
-
-trabecular_bone = Struct()
-trabecular_bone.vp = 2300.0
-trabecular_bone.rho = 1700.0
-trabecular_bone.alpha = 8.0
-trabecular_bone.render_color = "#EBD378"
-
-
-brain = Struct()
-brain.vp = 1560.0
-brain.rho = 1040.0
-brain.alpha = 0.3
-brain.render_color = "#DB504A"
-
-
-# these numbers are completely made up
-# TODO: research reasonable values
-tumor = Struct()
-tumor.vp = 1650.0
-tumor.rho = 1150.0
-tumor.alpha = 0.8
-tumor.render_color = "#94332F"
diff --git a/tests/neurotechdevkit/scenarios/test_base.py b/tests/neurotechdevkit/scenarios/test_base.py
index c4454024..60fbd46b 100644
--- a/tests/neurotechdevkit/scenarios/test_base.py
+++ b/tests/neurotechdevkit/scenarios/test_base.py
@@ -5,7 +5,6 @@
 import pytest
 import stride
 from frozenlist import FrozenList
-from mosaic.types import Struct
 
 from neurotechdevkit.results import PulsedResult, SteadyStateResult
 from neurotechdevkit.scenarios import materials
@@ -32,6 +31,7 @@ class ScenarioBaseTester(Scenario):
             description="bar",
         ),
     }
+    _material_layers = [("foo", materials.brain), ("bar", materials.skin)]
 
     default_source = PlanarSource3D(
         position=np.array([0.05, 0.1, 0.05]),
@@ -47,10 +47,6 @@ def __init__(self, complexity="fast"):
             complexity=complexity,
         )
 
-    @property
-    def _material_layers(self) -> list[tuple[str, Struct]]:
-        return [("foo", materials.brain), ("bar", materials.skin)]
-
     @property
     def _material_outline_upsample_factor(self) -> int:
         return 3

From 3ab0ac2c97cd3d43f8ad398eb56e44fb6a8b46e8 Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Thu, 15 Jun 2023 14:52:58 -0300
Subject: [PATCH 02/17] Removing docstring mentions of the fixed 500kHz
 frequency

---
 src/neurotechdevkit/scenarios/_base.py | 47 +++-----------------------
 1 file changed, 4 insertions(+), 43 deletions(-)

diff --git a/src/neurotechdevkit/scenarios/_base.py b/src/neurotechdevkit/scenarios/_base.py
index bb02ceac..5e6efeca 100644
--- a/src/neurotechdevkit/scenarios/_base.py
+++ b/src/neurotechdevkit/scenarios/_base.py
@@ -369,10 +369,6 @@ def simulate_steady_state(
         found by taking the Fourier transform of the last `n_cycles_steady_state` cycles
         of data and taking the amplitude of the component at the `center_frequency`.
 
-        !!! note
-            The only supported frequency currently supported is 500kHz. Any other
-            value will raise a NotImplementedError.
-
         !!! warning
             A poor choice of arguments to this function can lead to a failed
             simulation. Make sure you understand the impact of supplying parameter
@@ -380,8 +376,7 @@ def simulate_steady_state(
 
         Args:
             center_frequency: the center frequency (in hertz) to use for the
-                continuous-wave source output. No other value besides 500kHz (the
-                default) is currently supported.
+                continuous-wave source output.
             points_per_period: the number of points in time to simulate for each cycle
                 of the wave.
             n_cycles_steady_state: the number of complete cycles to use when calculating
@@ -397,9 +392,6 @@ def simulate_steady_state(
             n_jobs: the number of threads to be used for the computation. Use None to
                 leverage Devito automatic tuning.
 
-        Raises:
-            NotImplementedError: if a `center_frequency` other than 500kHz is provided.
-
         Returns:
             An object containing the result of the steady-state simulation.
         """
@@ -464,10 +456,6 @@ def simulate_pulse(
         In this simulation, the sources will emit a pulse containing a few cycles of
         oscillation and then let the pulse propagate out to all edges of the scenario.
 
-        !!! note
-            The only supported frequency currently supported is 500kHz. Any other
-            value will raise a NotImplementedError.
-
         !!! warning
             A poor choice of arguments to this function can lead to a failed
             simulation. Make sure you understand the impact of supplying parameter
@@ -475,8 +463,7 @@ def simulate_pulse(
 
         Args:
             center_frequency: the center frequency (in hertz) to use for the
-                continuous-wave source output. No other value besides
-                500kHz (the default) is currently supported.
+                continuous-wave source output.
             points_per_period: the number of points in time to simulate for each cycle
                 of the wave.
             simulation_time: the amount of time (in seconds) the simulation should run.
@@ -489,9 +476,6 @@ def simulate_pulse(
             n_jobs: the number of threads to be used for the computation. Use None to
                 leverage Devito automatic tuning.
 
-        Raises:
-            NotImplementedError: if a `center_frequency` other than 500kHz is provided.
-
         Returns:
             An object containing the result of the 2D pulsed simulation.
         """
@@ -520,18 +504,13 @@ def _simulate_pulse(
         In this simulation, the sources will emit a pulse containing a few cycles of
         oscillation and then let the pulse propagate out to all edges of the scenario.
 
-        !!! note
-            The only supported frequency currently supported is 500kHz. Any other
-            value will raise a NotImplementedError.
-
         Warning: A poor choice of arguments to this function can lead to a failed
         simulation. Make sure you understand the impact of supplying parameter values
         other than the default if you chose to do so.
 
         Args:
             center_frequency: the center frequency (in hertz) to use for the
-                continuous-wave source output. No other value besides
-                500kHz (the default) is currently supported.
+                continuous-wave source output.
             points_per_period: the number of points in time to simulate for each cycle
                 of the wave.
             simulation_time: the amount of time (in seconds) the simulation should run.
@@ -550,19 +529,9 @@ def _simulate_pulse(
                 which the slice of the 3D field should be made. Only valid if
                 `slice_axis` is not None.
 
-        Raises:
-            NotImplementedError: if a `center_frequency` other than 500kHz is provided.
-
         Returns:
             An object containing the result of the pulsed simulation.
         """
-        if center_frequency != 5.0e5:
-            raise NotImplementedError(
-                "500kHz is the only currently supported center frequency. Support for"
-                " other frequencies will be implemented once material properties as a"
-                " function of frequency has been implemented."
-            )
-
         problem = self.problem
         if simulation_time is None:
             simulation_time = select_simulation_time_for_pulsed(
@@ -998,10 +967,6 @@ def simulate_pulse(
         In this simulation, the sources will emit a pulse containing a few cycles of
         oscillation and then let the pulse propagate out to all edges of the scenario.
 
-        !!! note
-            The only supported frequency currently supported is 500kHz. Any
-            other value will raise a NotImplementedError.
-
         !!! warning
             A poor choice of arguments to this function can lead to a failed
             simulation. Make sure you understand the impact of supplying parameter
@@ -1009,8 +974,7 @@ def simulate_pulse(
 
         Args:
             center_frequency: the center frequency (in hertz) to use for the
-                continuous-wave source output. No other value besides
-                500kHz (the default) is currently supported.
+                continuous-wave source output.
             points_per_period: the number of points in time to simulate for each cycle
                 of the wave.
             simulation_time: the amount of time (in seconds) the simulation should run.
@@ -1029,9 +993,6 @@ def simulate_pulse(
                 which the slice of the 3D field should be made. Only valid if
                 `slice_axis` is not None.
 
-        Raises:
-            NotImplementedError: if a `center_frequency` other than 500kHz is provided.
-
         Returns:
             An object containing the result of the 3D pulsed simulation.
         """

From a0bfafeb94a17b667c0bea8c88c215db5f7f0626 Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Mon, 19 Jun 2023 16:45:27 +0100
Subject: [PATCH 03/17] Making materials dynamic, allowing the absorption
 calculation based on center frequency

---
 .vscode/settings.json                         |  5 ++
 .../plot_customized_center_frequency.py       | 13 ++--
 docs/examples/plot_full_scenario.py           | 22 +++---
 src/neurotechdevkit/materials.py              | 74 ++++++++++++++++---
 src/neurotechdevkit/results/_metrics.py       |  2 +-
 src/neurotechdevkit/scenarios/_base.py        | 63 ++++++++++------
 src/neurotechdevkit/scenarios/_scenario_0.py  | 17 ++---
 src/neurotechdevkit/scenarios/_scenario_1.py  | 72 ++++++------------
 src/neurotechdevkit/scenarios/_scenario_2.py  | 64 ++++++----------
 tests/neurotechdevkit/scenarios/test_base.py  |  7 +-
 .../neurotechdevkit/scenarios/test_metrics.py |  2 +-
 tests/neurotechdevkit/scenarios/test_time.py  | 16 ++--
 12 files changed, 197 insertions(+), 160 deletions(-)

diff --git a/.vscode/settings.json b/.vscode/settings.json
index 7a73a41b..1be6a53a 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -1,2 +1,7 @@
 {
+  "python.testing.pytestArgs": [
+    "tests"
+  ],
+  "python.testing.unittestEnabled": false,
+  "python.testing.pytestEnabled": true
 }
\ No newline at end of file
diff --git a/docs/examples/plot_customized_center_frequency.py b/docs/examples/plot_customized_center_frequency.py
index 4b16e4ba..bfca208d 100644
--- a/docs/examples/plot_customized_center_frequency.py
+++ b/docs/examples/plot_customized_center_frequency.py
@@ -14,12 +14,13 @@
 scenario = ndk.make('scenario-0-v0')
 
 # Customizing material layers with random values
-scenario._material_layers = [
-  ("water", ndk.materials.Material(1600.0, 1100.0, 0.0, "#2E86AB")),
-  ("skull", ndk.materials.Material(3000.0, 1850.0, 4.0, "#FAF0CA")),
-  ("brain", ndk.materials.Material(1660.0, 1140.0, 0.3, "#DB504A")),
-  ("tumor", ndk.materials.Material(1850.0, 1250.0, 0.8, "#94332F")),
-]
+scenario.material_layers = ["water", "cortical_bone", "brain", "tumor"]
+scenario.material_properties = {
+  "water": ndk.materials.Material(1600.0, 1100.0, 0.0, "#2E86AB"),
+  "cortical_bone": ndk.materials.Material(3000.0, 1850.0, 4.0, "#FAF0CA"),
+  "brain": ndk.materials.Material(1660.0, 1140.0, 0.3, "#DB504A"),
+  "tumor": ndk.materials.Material(1850.0, 1250.0, 0.8, "#94332F"),
+}
 
 result = scenario.simulate_steady_state(center_frequency=CENTER_FREQUENCY)
 result.render_steady_state_amplitudes(show_material_outlines=False)
diff --git a/docs/examples/plot_full_scenario.py b/docs/examples/plot_full_scenario.py
index 2774bd1c..05089908 100644
--- a/docs/examples/plot_full_scenario.py
+++ b/docs/examples/plot_full_scenario.py
@@ -13,10 +13,8 @@
 import numpy as np
 import stride
 
-from mosaic.types import Struct
-
 from neurotechdevkit import sources
-from neurotechdevkit.scenarios import Scenario2D, Target, make_grid, add_material_fields_to_problem, materials
+from neurotechdevkit.scenarios import Scenario2D, Target, make_grid, add_material_fields_to_problem
 
 
 class FullScenario(Scenario2D):
@@ -51,12 +49,12 @@ class FullScenario(Scenario2D):
     - render_color: the color used when rendering this material in the scenario layout
         plot.
     """
-    _material_layers = [
-        ("water", materials.water),
-        ("skin", materials.skin),
-        ("cortical bone", materials.cortical_bone),
-        ("trabecular bone", materials.trabecular_bone),
-        ("brain", materials.brain),
+    material_layers = [
+        "water",
+        "skin",
+        "cortical_bone",
+        "trabecular_bone",
+        "brain",
     ]
 
     def __init__(self, complexity="fast"):
@@ -103,7 +101,7 @@ def _compile_problem(self) -> stride.Problem:
             layer_ids=self.layer_ids,
             masks={
                 name: _create_scenario_1_mask(name, problem.grid)
-                for name in self.materials.keys()
+                for name in self.material_layers
             },
         )
         return problem
@@ -143,11 +141,11 @@ def _create_scenario_1_mask(material, grid):
     elif material == "skin":
         _fill_mask(mask, start=interfaces[0], end=interfaces[1], dx=dx)
 
-    elif material == "cortical bone":
+    elif material == "cortical_bone":
         _fill_mask(mask, start=interfaces[1], end=interfaces[2], dx=dx)
         _fill_mask(mask, start=interfaces[3], end=interfaces[4], dx=dx)
 
-    elif material == "trabecular bone":
+    elif material == "trabecular_bone":
         _fill_mask(mask, start=interfaces[2], end=interfaces[3], dx=dx)
 
     elif material == "brain":
diff --git a/src/neurotechdevkit/materials.py b/src/neurotechdevkit/materials.py
index b7db740c..cd09f31b 100644
--- a/src/neurotechdevkit/materials.py
+++ b/src/neurotechdevkit/materials.py
@@ -27,12 +27,68 @@ def to_struct(self) -> Struct:
         return struct
 
 
-water = Material(1500.0, 1000.0, 0.0, "#2E86AB")
-skin = Material(1610.0, 1090.0, 0.2, "#FA8B53")
-cortical_bone = Material(2800.0, 1850.0, 4.0, "#FAF0CA")
-trabecular_bone = Material(2300.0, 1700.0, 8.0, "#EBD378")
-brain = Material(1560.0, 1040.0, 0.3, "#DB504A")
-
-# these numbers are completely made up
-# TODO: research reasonable values
-tumor = Material(1650.0, 1150.0, 0.8, "#94332F")
+# The values below consider a center frequency of 500kHz
+SUPPORTED_MATERIALS = {
+    "water": Material(1500.0, 1000.0, 0.0, "#2E86AB"),
+    "skin": Material(1610.0, 1090.0, 0.2, "#FA8B53"),
+    "cortical_bone": Material(2800.0, 1850.0, 4.0, "#FAF0CA"),
+    "trabecular_bone": Material(2300.0, 1700.0, 8.0, "#EBD378"),
+    "brain": Material(1560.0, 1040.0, 0.3, "#DB504A"),
+    # these numbers are completely made up
+    # TODO: research reasonable values
+    "tumor": Material(1650.0, 1150.0, 0.8, "#94332F"),
+}
+
+
+def _calculate_absorption(material_name: str) -> float:
+    # TODO: calculate alpha
+    return SUPPORTED_MATERIALS[material_name].alpha
+
+
+def get_render_color(material_name: str) -> str:
+    """
+    Get the render color for a material.
+
+    Args:
+        material_name (str): the name of the material. E.g. "water"
+
+    Raises:
+        ValueError: raised if the material is not supported.
+
+    Returns:
+        str: the render color of the material.
+    """
+    if material_name not in SUPPORTED_MATERIALS:
+        raise ValueError(f"Unsupported material: {material_name}")
+
+    return SUPPORTED_MATERIALS[material_name].render_color
+
+
+def get_material(material_name: str, center_frequency: float = 5.0e5) -> Material:
+    """
+    Get a material with properties used in the neurotechdevkit scenarios.
+
+    Args:
+        material_name (str): the name of the material. E.g. "water"
+        center_frequency (float, optional): The center frequency of the
+            transducer. Defaults to 5.0e5.
+
+    Raises:
+        ValueError: raised if the material is not supported.
+
+    Returns:
+        Material: a material with properties used in the neurotechdevkit.
+    """
+    if material_name not in SUPPORTED_MATERIALS:
+        raise ValueError(f"Unsupported material: {material_name}")
+
+    base_material = SUPPORTED_MATERIALS[material_name]
+    if center_frequency == 5.0e5:
+        return base_material
+
+    return Material(
+        vp=base_material.vp,
+        rho=base_material.rho,
+        alpha=_calculate_absorption(material_name),
+        render_color=base_material.render_color,
+    )
diff --git a/src/neurotechdevkit/results/_metrics.py b/src/neurotechdevkit/results/_metrics.py
index 63ec8813..ece788d7 100644
--- a/src/neurotechdevkit/results/_metrics.py
+++ b/src/neurotechdevkit/results/_metrics.py
@@ -113,7 +113,7 @@ def calculate_all_metrics(
                     " intended use)."
                 ),
             }
-            for layer in result.scenario.ordered_layers
+            for layer in result.scenario.material_layers
         },
     }
 
diff --git a/src/neurotechdevkit/scenarios/_base.py b/src/neurotechdevkit/scenarios/_base.py
index 5e6efeca..bc2468f2 100644
--- a/src/neurotechdevkit/scenarios/_base.py
+++ b/src/neurotechdevkit/scenarios/_base.py
@@ -16,7 +16,7 @@
 from stride.problem import StructuredData
 
 from .. import rendering, results
-from ..materials import Material
+from ..materials import Material, get_material, get_render_color
 from ..sources import Source
 from ._resources import budget_time_and_memory_resources
 from ._shots import create_shot
@@ -60,7 +60,12 @@ class Scenario(abc.ABC):
 
     _SCENARIO_ID: str
     _TARGET_OPTIONS: dict[str, Target]
-    _material_layers: list[tuple[str, Material]]
+
+    # The list of material layers in the scenario.
+    material_layers: list[str]
+
+    # The customizations to the material layers.
+    material_properties: dict[str, Material] = {}
 
     def __init__(
         self,
@@ -216,9 +221,8 @@ def target_radius(self) -> float:
         """The radius of the target region (in meters)."""
         return self.target.radius
 
-    @property
-    def materials(self) -> Mapping[str, Struct]:
-        """A map between material name and material properties.
+    def get_materials(self, center_frequency: float = 5.0e5) -> Mapping[str, Struct]:
+        """Return a map between material name and material properties.
 
         - vp: the speed of sound (in m/s).
         - rho: the mass density (in kg/m³).
@@ -226,17 +230,36 @@ def materials(self) -> Mapping[str, Struct]:
         - render_color: the color used when rendering this material in the
         scenario layout plot.
         """
-        return {name: material.to_struct() for name, material in self._material_layers}
+        materials = {}
+        for layer in self.material_layers:
+            if layer not in self.material_properties:
+                material_properties = get_material(layer, center_frequency)
+            else:
+                material_properties = self.material_properties[layer]
+            materials[layer] = material_properties.to_struct()
+        return materials
 
     @property
-    def layer_ids(self) -> Mapping[str, int]:
-        """A map between material names and their layer id."""
-        return {name: n for n, (name, _) in enumerate(self._material_layers)}
+    def material_colors(self) -> dict[str, str]:
+        """
+        A map between material name and material render color.
+
+        Returns:
+            dict[str, str]: keys are material names and values are the hex color
+        """
+        material_colors = {}
+        for material in self.material_layers:
+            if material in self.material_properties:
+                color = self.material_properties[material].render_color
+            else:
+                color = get_render_color(material_name=material)
+            material_colors[material] = color
+        return material_colors
 
     @property
-    def ordered_layers(self) -> list[str]:
-        """A list of material names in order of their layer id."""
-        return [name for name, _ in self._material_layers]
+    def layer_ids(self) -> Mapping[str, int]:
+        """A map between material names and their layer id."""
+        return {name: n for n, name in enumerate(self.material_layers)}
 
     @property
     @abc.abstractmethod
@@ -398,7 +421,7 @@ def simulate_steady_state(
         problem = self.problem
         sim_time = select_simulation_time_for_steady_state(
             grid=problem.grid,
-            materials=self.materials,
+            materials=self.get_materials(center_frequency),
             freq_hz=center_frequency,
             time_to_steady_state=time_to_steady_state,
             n_cycles_steady_state=n_cycles_steady_state,
@@ -536,7 +559,7 @@ def _simulate_pulse(
         if simulation_time is None:
             simulation_time = select_simulation_time_for_pulsed(
                 grid=problem.grid,
-                materials=self.materials,
+                materials=self.get_materials(center_frequency),
                 delay=find_largest_delay_in_sources(self.sources),
             )
         problem.grid.time = create_time_grid(
@@ -863,9 +886,7 @@ def render_layout(
             show_material_outlines: whether or not to display a thin white outline of
                 the transition between different materials.
         """
-        color_sequence = [
-            self.materials[name].render_color for name in self.ordered_layers
-        ]
+        color_sequence = list(self.material_colors.values())
         field = self.get_field_data("layer").astype(int)
         fig, ax = rendering.create_layout_fig(
             self.extent, self.origin, color_sequence, field
@@ -898,7 +919,7 @@ def render_layout(
             fig=fig,
             ax=ax,
             color_sequence=color_sequence,
-            layer_labels=self.ordered_layers,
+            layer_labels=self.material_layers,
             show_sources=show_sources,
             show_target=show_target,
             extent=self.extent,
@@ -1038,9 +1059,7 @@ def render_layout(
         if slice_position is None:
             slice_position = self.get_default_slice_position(slice_axis)
 
-        color_sequence = [
-            self.materials[name].render_color for name in self.ordered_layers
-        ]
+        color_sequence = list(self.material_colors.values())
         field = self.get_field_data("layer").astype(int)
         field = slice_field(field, self, slice_axis, slice_position)
         extent = drop_element(self.extent, slice_axis)
@@ -1078,7 +1097,7 @@ def render_layout(
             fig=fig,
             ax=ax,
             color_sequence=color_sequence,
-            layer_labels=self.ordered_layers,
+            layer_labels=self.material_layers,
             show_sources=show_sources,
             show_target=show_target,
             extent=extent,
diff --git a/src/neurotechdevkit/scenarios/_scenario_0.py b/src/neurotechdevkit/scenarios/_scenario_0.py
index 29072edb..c8d70221 100644
--- a/src/neurotechdevkit/scenarios/_scenario_0.py
+++ b/src/neurotechdevkit/scenarios/_scenario_0.py
@@ -1,7 +1,6 @@
 import numpy as np
 import stride
 
-from .. import materials
 from ..sources import FocusedSource2D
 from ._base import Scenario2D, Target
 from ._utils import (
@@ -24,11 +23,11 @@ class Scenario0(Scenario2D):
             description="Represents a simulated tumor.",
         ),
     }
-    _material_layers = [
-        ("water", materials.water),
-        ("skull", materials.cortical_bone),
-        ("brain", materials.brain),
-        ("tumor", materials.tumor),
+    material_layers = [
+        "water",
+        "cortical_bone",
+        "brain",
+        "tumor",
     ]
 
     def __init__(self, complexity="fast"):
@@ -47,7 +46,7 @@ def _material_outline_upsample_factor(self) -> int:
     def _get_material_masks(self, problem):
         return {
             name: _create_scenario_0_mask(name, problem.grid, self._origin)
-            for name in self.materials.keys()
+            for name in self.material_layers
         }
 
     def _compile_problem(self) -> stride.Problem:
@@ -71,7 +70,7 @@ def _compile_problem(self) -> stride.Problem:
         )
         problem = add_material_fields_to_problem(
             problem=problem,
-            materials=self.materials,
+            materials=self.get_materials(c_freq),
             layer_ids=self.layer_ids,
             masks=self._get_material_masks(problem),
         )
@@ -94,7 +93,7 @@ def _create_scenario_0_mask(material, grid, origin):
         water_mask = ~outer_skull_mask
         return water_mask
 
-    elif material == "skull":
+    elif material == "cortical_bone":
         outer_skull_mask = _create_skull_interface_mask(grid, origin)
         outer_brain_mask = _create_brain_interface_mask(grid, origin)
         skull_mask = outer_skull_mask & ~outer_brain_mask
diff --git a/src/neurotechdevkit/scenarios/_scenario_1.py b/src/neurotechdevkit/scenarios/_scenario_1.py
index 7022744b..2c4c9d8a 100644
--- a/src/neurotechdevkit/scenarios/_scenario_1.py
+++ b/src/neurotechdevkit/scenarios/_scenario_1.py
@@ -1,44 +1,18 @@
 from __future__ import annotations
 
-from typing import Mapping, Protocol
+from typing import Mapping
 
 import numpy as np
 import numpy.typing as npt
 import stride
-from mosaic.types import Struct
 
-from .. import materials, rendering, sources
-from ._base import Scenario2D, Scenario3D, Target
+from .. import rendering, sources
+from ._base import Scenario, Scenario2D, Scenario3D, Target
 from ._utils import add_material_fields_to_problem, make_grid
 
 
-class _Scenario1MixinProtocol(Protocol):
-    """Provide type-hinting for Scenario 1 members used by mixins."""
-
-    @property
-    def scenario_id(self) -> str:
-        ...
-
-    @property
-    def complexity(self) -> str:
-        ...
-
-    @property
-    def materials(self) -> Mapping[str, Struct]:
-        ...
-
-    @property
-    def layer_ids(self) -> Mapping[str, int]:
-        ...
-
-    def _get_material_masks(
-        self, problem: stride.Problem
-    ) -> Mapping[str, npt.NDArray[np.bool_]]:
-        ...
-
-
-class _Scenario1Mixin:
-    """A mixin providing specific implementation detail for scenario 1.
+class Scenario1(Scenario):
+    """Specific implementation detail for scenario 1.
 
     Scenario 1 is based on benchmark 4 of the following paper:
 
@@ -50,12 +24,12 @@ class _Scenario1Mixin:
         https://asa.scitation.org/doi/pdf/10.1121/10.0013426
     """
 
-    _material_layers = [
-        ("water", materials.water),
-        ("skin", materials.skin),
-        ("cortical bone", materials.cortical_bone),
-        ("trabecular bone", materials.trabecular_bone),
-        ("brain", materials.brain),
+    material_layers = [
+        "water",
+        "skin",
+        "cortical_bone",
+        "trabecular_bone",
+        "brain",
     ]
 
     @property
@@ -63,15 +37,15 @@ def _material_outline_upsample_factor(self) -> int:
         return 8
 
     def _get_material_masks(
-        self: _Scenario1MixinProtocol, problem: stride.Problem
+        self, problem: stride.Problem
     ) -> Mapping[str, npt.NDArray[np.bool_]]:
         return {
             name: _create_scenario_1_mask(name, problem.grid)
-            for name in self.materials.keys()
+            for name in self.material_layers
         }
 
     def _compile_scenario_1_problem(
-        self: _Scenario1MixinProtocol, extent: npt.NDArray[np.float_]
+        self, extent: npt.NDArray[np.float_]
     ) -> stride.Problem:
         # scenario constants
         speed_water = 1500  # m/s
@@ -89,14 +63,14 @@ def _compile_scenario_1_problem(
         )
         problem = add_material_fields_to_problem(
             problem=problem,
-            materials=self.materials,
+            materials=self.get_materials(c_freq),
             layer_ids=self.layer_ids,
             masks=self._get_material_masks(problem),
         )
         return problem
 
 
-class Scenario1_2D(_Scenario1Mixin, Scenario2D):
+class Scenario1_2D(Scenario1, Scenario2D):
     """A 2D implementation of scenario 1.
 
     Scenario 1 is based on benchmark 4 of the following paper:
@@ -138,7 +112,7 @@ def get_default_source(self) -> sources.Source:
         )
 
 
-class Scenario1_3D(_Scenario1Mixin, Scenario3D):
+class Scenario1_3D(Scenario1, Scenario3D):
     """A 3D implementation of scenario 1.
 
     Scenario 1 is based on benchmark 4 of the following paper:
@@ -187,15 +161,15 @@ def viewer_config_3d(self) -> rendering.ViewerConfig3D:
             colormaps={
                 "water": "blue",
                 "skin": "viridis",
-                "cortical bone": "magma",
-                "trabecular bone": "inferno",
+                "cortical_bone": "magma",
+                "trabecular_bone": "inferno",
                 "brain": "bop orange",
             },
             opacities={
                 "water": 0.8,
                 "skin": 0.2,
-                "cortical bone": 0.2,
-                "trabecular bone": 0.2,
+                "cortical_bone": 0.2,
+                "trabecular_bone": 0.2,
                 "brain": 0.4,
             },
         )
@@ -249,11 +223,11 @@ def _create_scenario_1_mask(material, grid):
     elif material == "skin":
         _fill_mask(mask, start=interfaces[0], end=interfaces[1], dx=dx)
 
-    elif material == "cortical bone":
+    elif material == "cortical_bone":
         _fill_mask(mask, start=interfaces[1], end=interfaces[2], dx=dx)
         _fill_mask(mask, start=interfaces[3], end=interfaces[4], dx=dx)
 
-    elif material == "trabecular bone":
+    elif material == "trabecular_bone":
         _fill_mask(mask, start=interfaces[2], end=interfaces[3], dx=dx)
 
     elif material == "brain":
diff --git a/src/neurotechdevkit/scenarios/_scenario_2.py b/src/neurotechdevkit/scenarios/_scenario_2.py
index 4e9a0332..d3eab1ba 100644
--- a/src/neurotechdevkit/scenarios/_scenario_2.py
+++ b/src/neurotechdevkit/scenarios/_scenario_2.py
@@ -1,44 +1,20 @@
 from __future__ import annotations
 
 import pathlib
-from typing import Mapping, Protocol
+from typing import Mapping
 
 import hdf5storage
 import numpy as np
 import numpy.typing as npt
 import stride
-from mosaic.types import Struct
 
-from .. import materials, rendering, sources
-from ._base import Scenario2D, Scenario3D, Target
+from .. import rendering, sources
+from ._base import Scenario, Scenario2D, Scenario3D, Target
 from ._utils import add_material_fields_to_problem, make_grid
 
 
-class _Scenario2MixinProtocol(Protocol):
-    """Provide type-hinting for Scenario 2 members used by mixins."""
-
-    @property
-    def scenario_id(self) -> str:
-        ...
-
-    @property
-    def complexity(self) -> str:
-        ...
-
-    @property
-    def materials(self) -> Mapping[str, Struct]:
-        ...
-
-    @property
-    def layer_ids(self) -> Mapping[str, int]:
-        ...
-
-    def _get_material_masks(self) -> Mapping[str, npt.NDArray[np.bool_]]:
-        ...
-
-
-class _Scenario2Mixin:
-    """A mixin providing specific implementation detail for scenario 2.
+class Scenario2(Scenario):
+    """Specific implementation detail for scenario 2.
 
     Scenario 2 is based on benchmark 8 of the following paper:
 
@@ -54,8 +30,12 @@ class _Scenario2Mixin:
     def _material_outline_upsample_factor(self) -> int:
         return 4
 
+    def _get_material_masks(self) -> Mapping[str, npt.NDArray[np.bool_]]:
+        """Will be implemented by the subclass."""
+        raise NotImplementedError()
+
     def _compile_scenario_2_problem(
-        self: _Scenario2MixinProtocol, extent: npt.NDArray[np.float_]
+        self, extent: npt.NDArray[np.float_]
     ) -> stride.Problem:
         # scenario constants
         speed_water = 1500  # m/s
@@ -73,14 +53,14 @@ def _compile_scenario_2_problem(
         )
         problem = add_material_fields_to_problem(
             problem=problem,
-            materials=self.materials,
+            materials=self.get_materials(c_freq),
             layer_ids=self.layer_ids,
             masks=self._get_material_masks(),
         )
         return problem
 
 
-class Scenario2_2D(_Scenario2Mixin, Scenario2D):
+class Scenario2_2D(Scenario2, Scenario2D):
     """A 2D implementation of scenario 2.
 
     Scenario 2 is based on benchmark 8 of the following paper:
@@ -129,10 +109,10 @@ class Scenario2_2D(_Scenario2Mixin, Scenario2D):
             ),
         ),
     }
-    _material_layers = [
-        ("water", materials.water),
-        ("skull", materials.cortical_bone),
-        ("brain", materials.brain),
+    material_layers = [
+        "water",
+        "cortical_bone",
+        "brain",
     ]
 
     def __init__(self, complexity="fast"):
@@ -151,7 +131,7 @@ def _compile_problem(self) -> stride.Problem:
     def _get_material_masks(self) -> Mapping[str, npt.NDArray[np.bool_]]:
         return {
             name: _create_scenario_2_mask(name, convert_2d=True)
-            for name in self.materials.keys()
+            for name in self.material_layers
         }
 
     def get_default_source(self) -> sources.Source:
@@ -165,7 +145,7 @@ def get_default_source(self) -> sources.Source:
         )
 
 
-class Scenario2_3D(_Scenario2Mixin, Scenario3D):
+class Scenario2_3D(Scenario2, Scenario3D):
     """A 3D implementation of scenario 2.
 
     Scenario 2 is based on benchmark 8 of the following paper:
@@ -256,12 +236,12 @@ def viewer_config_3d(self) -> rendering.ViewerConfig3D:
             init_zoom=2.0,
             colormaps={
                 "water": "blue",
-                "skull": "magma",
+                "cortical_bone": "magma",
                 "brain": "bop orange",
             },
             opacities={
                 "water": 0.8,
-                "skull": 0.2,
+                "cortical_bone": 0.2,
                 "brain": 0.2,
             },
         )
@@ -282,7 +262,7 @@ def _compile_problem(self) -> stride.Problem:
     def _get_material_masks(self):
         return {
             name: _create_scenario_2_mask(name, convert_2d=False)
-            for name in self.materials.keys()
+            for name in self.material_layers
         }
 
     def get_default_source(self):
@@ -305,7 +285,7 @@ def _create_scenario_2_mask(material, convert_2d=False) -> npt.NDArray[np.bool_]
     skull_mask = mat_data["skull_mask"].astype(np.bool_)
     brain_mask = mat_data["brain_mask"].astype(np.bool_)
 
-    if material == "skull":
+    if material == "cortical_bone":
         mask = skull_mask
 
     elif material == "brain":
diff --git a/tests/neurotechdevkit/scenarios/test_base.py b/tests/neurotechdevkit/scenarios/test_base.py
index 60fbd46b..cb1977f5 100644
--- a/tests/neurotechdevkit/scenarios/test_base.py
+++ b/tests/neurotechdevkit/scenarios/test_base.py
@@ -7,7 +7,6 @@
 from frozenlist import FrozenList
 
 from neurotechdevkit.results import PulsedResult, SteadyStateResult
-from neurotechdevkit.scenarios import materials
 from neurotechdevkit.scenarios._base import Scenario, Target
 from neurotechdevkit.scenarios._utils import make_grid, wavelet_helper
 from neurotechdevkit.sources import FocusedSource3D, PlanarSource3D, Source
@@ -31,7 +30,7 @@ class ScenarioBaseTester(Scenario):
             description="bar",
         ),
     }
-    _material_layers = [("foo", materials.brain), ("bar", materials.skin)]
+    material_layers = ["brain", "skin"]
 
     default_source = PlanarSource3D(
         position=np.array([0.05, 0.1, 0.05]),
@@ -515,7 +514,7 @@ def test_get_layer_mask_with_wrong_layer_name(tester_with_layers):
 
 def test_get_layer_mask_with_first_layer(tester_with_layers):
     """Verify that get_layer_mask returns the expected mask for the first layer."""
-    mask = tester_with_layers.get_layer_mask("foo")
+    mask = tester_with_layers.get_layer_mask("brain")
     expected = np.zeros_like(mask, dtype=bool)
     expected[:5] = True
     np.testing.assert_allclose(mask, expected)
@@ -523,7 +522,7 @@ def test_get_layer_mask_with_first_layer(tester_with_layers):
 
 def test_get_layer_mask_with_last_layer(tester_with_layers):
     """Verify that get_layer_mask returns the expected mask for the last layer."""
-    mask = tester_with_layers.get_layer_mask("bar")
+    mask = tester_with_layers.get_layer_mask("skin")
     expected = np.zeros_like(mask, dtype=bool)
     expected[5:] = True
     np.testing.assert_allclose(mask, expected)
diff --git a/tests/neurotechdevkit/scenarios/test_metrics.py b/tests/neurotechdevkit/scenarios/test_metrics.py
index d6f56fef..baaf91fb 100644
--- a/tests/neurotechdevkit/scenarios/test_metrics.py
+++ b/tests/neurotechdevkit/scenarios/test_metrics.py
@@ -54,7 +54,7 @@ def fake_result(fake_steady_state_matrix):
     """Returns a fake SteadyStateResult2D that can be used for testing."""
     return SteadyStateResult2D(
         scenario=SimpleNamespace(
-            ordered_layers=["water", "brain"],
+            material_layers=["water", "brain"],
             get_layer_mask=fake_layer_masks,
             get_target_mask=lambda: fake_layer_masks("target"),
             get_field_data=fake_fields,
diff --git a/tests/neurotechdevkit/scenarios/test_time.py b/tests/neurotechdevkit/scenarios/test_time.py
index 30c4133a..d9f2d4cf 100644
--- a/tests/neurotechdevkit/scenarios/test_time.py
+++ b/tests/neurotechdevkit/scenarios/test_time.py
@@ -52,7 +52,7 @@ def test_select_simulation_time_for_steady_state_with_defined_time_to_ss(grid_fa
     When time_to_steady_state is provided, the function should use it to calculate the
     simulation time.
     """
-    test_materials = {"test": materials.brain}
+    test_materials = {"brain": materials.get_material("brain")}
     test_freq = 2.0e4
     n_cycles = 5
     test_time_ss = 1e-4
@@ -76,7 +76,10 @@ def test_select_simulation_time_for_steady_state_with_default_time_to_ss(grid_fa
     When time_to_steady_state is not provided, the function should estimate a value
     for it based on the size of the scenario and the lowest speed of sound.
     """
-    test_materials = {"fast": materials.cortical_bone, "slow": materials.brain}
+    test_materials = {
+        "fast": materials.get_material("cortical_bone"),
+        "slow": materials.get_material("brain"),
+    }
     test_freq = 2.0e4
     n_cycles = 5
     sim_time = select_simulation_time_for_steady_state(
@@ -89,14 +92,17 @@ def test_select_simulation_time_for_steady_state_with_default_time_to_ss(grid_fa
     )
     period = 1.0 / test_freq
     length = np.sqrt(0.3**2 + 0.4**2 + 0.5**2)
-    expected_time_ss = 2 * length / materials.brain.vp
+    expected_time_ss = 2 * length / materials.get_material("brain").vp
     expected_time = expected_time_ss + n_cycles * period
     np.testing.assert_allclose(sim_time, expected_time)
 
 
 def test_select_simulation_time_for_pulsed(grid_fake):
     """Verify that the returned simulation time is correct for a pulsed simulation."""
-    test_materials = {"fast": materials.cortical_bone, "slow": materials.brain}
+    test_materials = {
+        "fast": materials.get_material("cortical_bone"),
+        "slow": materials.get_material("brain"),
+    }
     source_delay = 31.0
     sim_time = select_simulation_time_for_pulsed(
         grid=grid_fake,
@@ -104,7 +110,7 @@ def test_select_simulation_time_for_pulsed(grid_fake):
         delay=source_delay,
     )
     length = np.sqrt(0.3**2 + 0.4**2 + 0.5**2)
-    expected_time = source_delay + length / materials.brain.vp
+    expected_time = source_delay + length / materials.get_material("brain").vp
     np.testing.assert_allclose(sim_time, expected_time)
 
 

From ac64c78a7a9a3c84c6c8c2bb205641075f59612a Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Thu, 22 Jun 2023 17:23:44 +0100
Subject: [PATCH 04/17] Adding algorithm to calculate the material absorption
 based on the frequency

---
 .../plot_customized_center_frequency.py       |  9 +--
 docs/examples/plot_full_scenario.py           |  7 --
 src/neurotechdevkit/materials.py              | 79 ++++++++++++++++---
 3 files changed, 73 insertions(+), 22 deletions(-)

diff --git a/docs/examples/plot_customized_center_frequency.py b/docs/examples/plot_customized_center_frequency.py
index bfca208d..42ead974 100644
--- a/docs/examples/plot_customized_center_frequency.py
+++ b/docs/examples/plot_customized_center_frequency.py
@@ -13,13 +13,12 @@
 
 scenario = ndk.make('scenario-0-v0')
 
-# Customizing material layers with random values
+# using default material layers
 scenario.material_layers = ["water", "cortical_bone", "brain", "tumor"]
+
+# Customizing material properties with random values
 scenario.material_properties = {
-  "water": ndk.materials.Material(1600.0, 1100.0, 0.0, "#2E86AB"),
-  "cortical_bone": ndk.materials.Material(3000.0, 1850.0, 4.0, "#FAF0CA"),
-  "brain": ndk.materials.Material(1660.0, 1140.0, 0.3, "#DB504A"),
-  "tumor": ndk.materials.Material(1850.0, 1250.0, 0.8, "#94332F"),
+  "tumor": ndk.materials.Material(vp=1850.0, rho=1250.0,alpha=0.8,render_color="#94332F"),
 }
 
 result = scenario.simulate_steady_state(center_frequency=CENTER_FREQUENCY)
diff --git a/docs/examples/plot_full_scenario.py b/docs/examples/plot_full_scenario.py
index 05089908..5b978a45 100644
--- a/docs/examples/plot_full_scenario.py
+++ b/docs/examples/plot_full_scenario.py
@@ -41,13 +41,6 @@ class FullScenario(Scenario2D):
 
     """
     The order of returned materials defines the layering of the scenario.
-
-    Each material has the following properties:
-    - vp: the speed of sound (in m/s).
-    - rho: the mass density (in kg/m³).
-    - alpha: the absorption (in dB/cm).
-    - render_color: the color used when rendering this material in the scenario layout
-        plot.
     """
     material_layers = [
         "water",
diff --git a/src/neurotechdevkit/materials.py b/src/neurotechdevkit/materials.py
index cd09f31b..46ecb4d0 100644
--- a/src/neurotechdevkit/materials.py
+++ b/src/neurotechdevkit/materials.py
@@ -29,20 +29,79 @@ def to_struct(self) -> Struct:
 
 # The values below consider a center frequency of 500kHz
 SUPPORTED_MATERIALS = {
-    "water": Material(1500.0, 1000.0, 0.0, "#2E86AB"),
-    "skin": Material(1610.0, 1090.0, 0.2, "#FA8B53"),
-    "cortical_bone": Material(2800.0, 1850.0, 4.0, "#FAF0CA"),
-    "trabecular_bone": Material(2300.0, 1700.0, 8.0, "#EBD378"),
-    "brain": Material(1560.0, 1040.0, 0.3, "#DB504A"),
+    "water": Material(vp=1500.0, rho=1000.0, alpha=0.0, render_color="#2E86AB"),
+    "skin": Material(vp=1610.0, rho=1090.0, alpha=0.2, render_color="#FA8B53"),
+    "cortical_bone": Material(vp=2800.0, rho=1850.0, alpha=4.0, render_color="#FAF0CA"),
+    "trabecular_bone": Material(
+        vp=2300.0, rho=1700.0, alpha=8.0, render_color="#EBD378"
+    ),
+    "brain": Material(vp=1560.0, rho=1040.0, alpha=0.3, render_color="#DB504A"),
     # these numbers are completely made up
     # TODO: research reasonable values
-    "tumor": Material(1650.0, 1150.0, 0.8, "#94332F"),
+    "tumor": Material(vp=1650.0, rho=1150.0, alpha=0.8, render_color="#94332F"),
 }
 
 
-def _calculate_absorption(material_name: str) -> float:
-    # TODO: calculate alpha
-    return SUPPORTED_MATERIALS[material_name].alpha
+@dataclass
+class AttenuationConstant:
+    """The parameters of the attenuation constant."""
+
+    a0: float  # α0 [Np/m/MHz]
+    b: float
+
+    def calculate_absorption(self, frequency: float) -> float:
+        """
+        Calculate the absorption coefficient for a given center frequency.
+
+        The absorption coefficient is calculated using the following formula:
+            α(f) = α0 * f^b
+
+        We convert the absorption coefficient from Np/m/Hz to dB/m/Hz
+        multiplying it by 8.6860000037.
+
+        Args:
+            frequency (float): in Hz
+
+        Returns:
+            float: attenuation in dB/cm/MHz
+        """
+        attenuation = self.a0 * (frequency / 1e6) ** self.b  # Np/m
+        db_attenuation = attenuation * 8.6860000037  # dB/m
+        return db_attenuation / 100  # convert from dB/m to dB/cm
+
+
+ATTENUATION_CONSTANTS = {
+    "trabecular_bone": AttenuationConstant(a0=47, b=1.2),
+    "cortical_bone": AttenuationConstant(a0=54.553, b=1),
+    "skin": AttenuationConstant(a0=21.158, b=1),
+    "water": AttenuationConstant(a0=0.025328436, b=1),
+    "brain": AttenuationConstant(a0=6.8032, b=1.3),
+    # these numbers are completely made up
+    # TODO: research reasonable values
+    "tumor": AttenuationConstant(a0=8, b=1),
+}
+
+
+def _calculate_absorption(material_name: str, center_frequency: float) -> float:
+    """
+    Calculate the absorption coefficient for a given center frequency.
+
+    Args:
+        material_name (str): the name of the material. E.g. "water"
+        center_frequency (float): the center frequency of the transducer in Hz
+
+    Raises:
+        ValueError: raised if the material is not supported.
+
+    Returns:
+        float: the absorption coefficient in dB/cm/Hz
+    """
+    try:
+        return ATTENUATION_CONSTANTS[material_name].calculate_absorption(
+            center_frequency
+        )
+    except KeyError:
+        raise ValueError(f"Unsupported material: {material_name}")
 
 
 def get_render_color(material_name: str) -> str:
@@ -89,6 +148,6 @@ def get_material(material_name: str, center_frequency: float = 5.0e5) -> Materia
     return Material(
         vp=base_material.vp,
         rho=base_material.rho,
-        alpha=_calculate_absorption(material_name),
+        alpha=_calculate_absorption(material_name, center_frequency),
         render_color=base_material.render_color,
     )

From f747f3e92a09dd316ef3cea36a622bc8472df3dc Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Thu, 22 Jun 2023 17:50:56 +0100
Subject: [PATCH 05/17] Removing .vscode settings

---
 .vscode/settings.json | 7 -------
 1 file changed, 7 deletions(-)
 delete mode 100644 .vscode/settings.json

diff --git a/.vscode/settings.json b/.vscode/settings.json
deleted file mode 100644
index 1be6a53a..00000000
--- a/.vscode/settings.json
+++ /dev/null
@@ -1,7 +0,0 @@
-{
-  "python.testing.pytestArgs": [
-    "tests"
-  ],
-  "python.testing.unittestEnabled": false,
-  "python.testing.pytestEnabled": true
-}
\ No newline at end of file

From e6167c0b6445bc241870f6a097b29fff4e4917c1 Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Thu, 22 Jun 2023 17:52:50 +0100
Subject: [PATCH 06/17] Adding tests for materials

---
 .../scenarios/test_materials.py               | 124 ++++++++++++++++++
 1 file changed, 124 insertions(+)
 create mode 100644 tests/neurotechdevkit/scenarios/test_materials.py

diff --git a/tests/neurotechdevkit/scenarios/test_materials.py b/tests/neurotechdevkit/scenarios/test_materials.py
new file mode 100644
index 00000000..0c706350
--- /dev/null
+++ b/tests/neurotechdevkit/scenarios/test_materials.py
@@ -0,0 +1,124 @@
+import numpy as np
+import pytest
+from mosaic.types import Struct
+
+from neurotechdevkit.materials import SUPPORTED_MATERIALS, Material
+from neurotechdevkit.scenarios import Scenario
+
+
+def compare_structs(struct1: Struct, struct2: Struct):
+    """Compare two Structs."""
+    assert struct1.vp == struct2.vp
+    assert struct1.rho == struct2.rho
+    assert struct1.alpha == struct2.alpha
+    assert struct1.render_color == struct2.render_color
+
+
+class BaseScenario(Scenario):
+    """A scenario for testing the materials module."""
+
+    _SCENARIO_ID = "scenario-tester"
+    _TARGET_OPTIONS = {}
+
+    def __init__(self, complexity="fast"):
+        self._target_id = "target_1"
+        super().__init__(
+            origin=np.array([-0.1, -0.1, 0.0]),
+            complexity=complexity,
+        )
+
+    def _compile_problem(self):
+        pass
+
+    def _material_outline_upsample_factor(self):
+        pass
+
+    def get_default_source(self):
+        pass
+
+    def get_target_mask(self):
+        pass
+
+
+def test_default_materials():
+    """Test that the default materials are used."""
+
+    class ScenarioWithDefaultMaterials(BaseScenario):
+        material_layers = ["brain", "skin"]
+
+    scenario = ScenarioWithDefaultMaterials()
+    assert scenario.material_colors == {"brain": "#DB504A", "skin": "#FA8B53"}
+
+    materials = scenario.get_materials()
+    assert list(materials.keys()) == ["brain", "skin"]
+    compare_structs(materials["brain"], SUPPORTED_MATERIALS["brain"].to_struct())
+    compare_structs(materials["skin"], SUPPORTED_MATERIALS["skin"].to_struct())
+
+
+def test_custom_material_property():
+    """Test that a custom material property is used."""
+
+    class ScenarioWithCustomMaterialProperties(BaseScenario):
+        material_layers = ["brain"]
+        material_properties = {
+            "brain": Material(vp=1600.0, rho=1100.0, alpha=0.0, render_color="#2E86AB")
+        }
+
+    scenario = ScenarioWithCustomMaterialProperties()
+    assert scenario.material_colors == {"brain": "#2E86AB"}
+
+    materials = scenario.get_materials()
+    assert list(materials.keys()) == ["brain"]
+    compare_structs(
+        materials["brain"],
+        Material(vp=1600.0, rho=1100.0, alpha=0.0, render_color="#2E86AB").to_struct(),
+    )
+
+
+def test_new_material():
+    """Test that a new material is used."""
+
+    class ScenarioWithCustomMaterial(BaseScenario):
+        material_layers = ["brain", "eye"]
+        material_properties = {
+            "eye": Material(vp=1600.0, rho=1100.0, alpha=0.0, render_color="#2E86AB")
+        }
+
+    scenario = ScenarioWithCustomMaterial()
+    assert scenario.material_colors == {"brain": "#DB504A", "eye": "#2E86AB"}
+
+    materials = scenario.get_materials()
+
+    assert list(materials.keys()) == ["brain", "eye"]
+    compare_structs(materials["brain"], SUPPORTED_MATERIALS["brain"].to_struct())
+    compare_structs(
+        materials["eye"],
+        Material(vp=1600.0, rho=1100.0, alpha=0.0, render_color="#2E86AB").to_struct(),
+    )
+
+
+def test_material_absorption_is_calculated():
+    """Test that the material absorption is calculated for a frequency !=500e3."""
+
+    class ScenarioWithBrainMaterial(BaseScenario):
+        material_layers = ["brain"]
+
+    scenario = ScenarioWithBrainMaterial()
+
+    materials = scenario.get_materials(600e3)
+
+    assert list(materials.keys()) == ["brain"]
+    assert materials["brain"].alpha == 0.3041793231753331
+
+
+def test_unknown_material_without_properties():
+    """Test that an unknown material without properties raises an error."""
+
+    class ScenarioWithCustomMaterial(BaseScenario):
+        material_layers = ["unknown_material"]
+
+    scenario = ScenarioWithCustomMaterial()
+    with pytest.raises(ValueError):
+        scenario.material_colors
+    with pytest.raises(ValueError):
+        scenario.get_materials()

From 238a935d36e629510a5877d369fc630282d68f13 Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Fri, 23 Jun 2023 08:24:54 +0100
Subject: [PATCH 07/17] Moving scenario2 layers definition to scenario2 class

---
 src/neurotechdevkit/scenarios/_scenario_2.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/neurotechdevkit/scenarios/_scenario_2.py b/src/neurotechdevkit/scenarios/_scenario_2.py
index d3eab1ba..d43558f7 100644
--- a/src/neurotechdevkit/scenarios/_scenario_2.py
+++ b/src/neurotechdevkit/scenarios/_scenario_2.py
@@ -26,6 +26,12 @@ class Scenario2(Scenario):
         https://asa.scitation.org/doi/pdf/10.1121/10.0013426
     """
 
+    material_layers = [
+        "water",
+        "cortical_bone",
+        "brain",
+    ]
+
     @property
     def _material_outline_upsample_factor(self) -> int:
         return 4
@@ -109,11 +115,6 @@ class Scenario2_2D(Scenario2, Scenario2D):
             ),
         ),
     }
-    material_layers = [
-        "water",
-        "cortical_bone",
-        "brain",
-    ]
 
     def __init__(self, complexity="fast"):
         """Create a new instance of scenario 2."""

From 25ac093cb886cd24685289225bfa3abdd907d7c4 Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Fri, 23 Jun 2023 08:36:32 +0100
Subject: [PATCH 08/17] Updating full scenario example

---
 docs/examples/plot_full_scenario.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/examples/plot_full_scenario.py b/docs/examples/plot_full_scenario.py
index 5b978a45..5210136d 100644
--- a/docs/examples/plot_full_scenario.py
+++ b/docs/examples/plot_full_scenario.py
@@ -90,7 +90,7 @@ def _compile_problem(self) -> stride.Problem:
         )
         problem = add_material_fields_to_problem(
             problem=problem,
-            materials=self.materials,
+            materials=self.get_materials(c_freq),
             layer_ids=self.layer_ids,
             masks={
                 name: _create_scenario_1_mask(name, problem.grid)

From 24a8696ebd02f2806b29a9ba39682c960a11cc04 Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Wed, 5 Jul 2023 10:36:57 +0200
Subject: [PATCH 09/17] Removing hard coded c_freq from methods

---
 docs/examples/plot_full_scenario.py               | 15 ++++++++++-----
 src/neurotechdevkit/scenarios/_base.py            |  7 ++++---
 src/neurotechdevkit/scenarios/_scenario_0.py      |  7 +++----
 src/neurotechdevkit/scenarios/_scenario_1.py      | 15 +++++++--------
 src/neurotechdevkit/scenarios/_scenario_2.py      | 15 +++++++--------
 tests/neurotechdevkit/scenarios/test_base.py      |  3 ++-
 tests/neurotechdevkit/scenarios/test_materials.py | 10 +++++-----
 7 files changed, 38 insertions(+), 34 deletions(-)

diff --git a/docs/examples/plot_full_scenario.py b/docs/examples/plot_full_scenario.py
index 5210136d..1acdea82 100644
--- a/docs/examples/plot_full_scenario.py
+++ b/docs/examples/plot_full_scenario.py
@@ -14,7 +14,12 @@
 import stride
 
 from neurotechdevkit import sources
-from neurotechdevkit.scenarios import Scenario2D, Target, make_grid, add_material_fields_to_problem
+from neurotechdevkit.scenarios import (
+    Scenario2D,
+    Target,
+    make_grid,
+    add_material_fields_to_problem,
+)
 
 
 class FullScenario(Scenario2D):
@@ -26,6 +31,7 @@ class FullScenario(Scenario2D):
     doi: 10.1121/10.0013426
     https://asa.scitation.org/doi/pdf/10.1121/10.0013426
     """
+
     _SCENARIO_ID = "the_id_for_this_scenario"
 
     """
@@ -71,18 +77,17 @@ def _material_outline_upsample_factor(self) -> int:
         """
         return 8
 
-    def _compile_problem(self) -> stride.Problem:
+    def _compile_problem(self, center_frequency) -> stride.Problem:
         """The problem definition for the scenario."""
         extent = np.array([0.12, 0.07])  # m
         # scenario constants
         speed_water = 1500  # m/s
-        c_freq = 500e3  # hz
 
         # desired resolution for complexity=fast
         ppw = 6
 
         # compute resolution
-        dx = speed_water / c_freq / ppw  # m
+        dx = speed_water / center_frequency / ppw  # m
 
         grid = make_grid(extent=extent, dx=dx)
         problem = stride.Problem(
@@ -90,7 +95,7 @@ def _compile_problem(self) -> stride.Problem:
         )
         problem = add_material_fields_to_problem(
             problem=problem,
-            materials=self.get_materials(c_freq),
+            materials=self.get_materials(center_frequency),
             layer_ids=self.layer_ids,
             masks={
                 name: _create_scenario_1_mask(name, problem.grid)
diff --git a/src/neurotechdevkit/scenarios/_base.py b/src/neurotechdevkit/scenarios/_base.py
index bc2468f2..4110dcca 100644
--- a/src/neurotechdevkit/scenarios/_base.py
+++ b/src/neurotechdevkit/scenarios/_base.py
@@ -78,7 +78,6 @@ def __init__(
             raise ValueError("the only complexity currently supported is 'fast'")
 
         self._origin = origin
-        self._problem = self._compile_problem()
         self._sources: FrozenList[Source] = FrozenList()
         self._target_id: str
 
@@ -221,7 +220,7 @@ def target_radius(self) -> float:
         """The radius of the target region (in meters)."""
         return self.target.radius
 
-    def get_materials(self, center_frequency: float = 5.0e5) -> Mapping[str, Struct]:
+    def get_materials(self, center_frequency=float) -> Mapping[str, Struct]:
         """Return a map between material name and material properties.
 
         - vp: the speed of sound (in m/s).
@@ -340,7 +339,7 @@ def get_field_data(self, field: str) -> npt.NDArray[np.float_]:
         return self.problem.medium.fields[field].data
 
     @abc.abstractmethod
-    def _compile_problem(self) -> stride.Problem:
+    def _compile_problem(self, center_frequency: float) -> stride.Problem:
         pass
 
     def reset(self) -> None:
@@ -418,6 +417,7 @@ def simulate_steady_state(
         Returns:
             An object containing the result of the steady-state simulation.
         """
+        self._problem = self._compile_problem(center_frequency)
         problem = self.problem
         sim_time = select_simulation_time_for_steady_state(
             grid=problem.grid,
@@ -555,6 +555,7 @@ def _simulate_pulse(
         Returns:
             An object containing the result of the pulsed simulation.
         """
+        self._problem = self._compile_problem(center_frequency)
         problem = self.problem
         if simulation_time is None:
             simulation_time = select_simulation_time_for_pulsed(
diff --git a/src/neurotechdevkit/scenarios/_scenario_0.py b/src/neurotechdevkit/scenarios/_scenario_0.py
index c8d70221..29c419c7 100644
--- a/src/neurotechdevkit/scenarios/_scenario_0.py
+++ b/src/neurotechdevkit/scenarios/_scenario_0.py
@@ -49,19 +49,18 @@ def _get_material_masks(self, problem):
             for name in self.material_layers
         }
 
-    def _compile_problem(self) -> stride.Problem:
+    def _compile_problem(self, center_frequency=float) -> stride.Problem:
         extent = np.array([0.05, 0.04])  # m
         origin = self.origin  # m
 
         # scenario constants
         speed_water = 1500  # m/s
-        c_freq = 500e3  # hz
 
         # desired resolution for complexity=fast
         ppw = 6
 
         # compute resolution
-        dx = speed_water / c_freq / ppw  # m
+        dx = speed_water / center_frequency / ppw  # m
 
         grid = make_grid(extent=extent, dx=dx)
         self._origin = origin
@@ -70,7 +69,7 @@ def _compile_problem(self) -> stride.Problem:
         )
         problem = add_material_fields_to_problem(
             problem=problem,
-            materials=self.get_materials(c_freq),
+            materials=self.get_materials(center_frequency),
             layer_ids=self.layer_ids,
             masks=self._get_material_masks(problem),
         )
diff --git a/src/neurotechdevkit/scenarios/_scenario_1.py b/src/neurotechdevkit/scenarios/_scenario_1.py
index 2c4c9d8a..87125000 100644
--- a/src/neurotechdevkit/scenarios/_scenario_1.py
+++ b/src/neurotechdevkit/scenarios/_scenario_1.py
@@ -45,17 +45,16 @@ def _get_material_masks(
         }
 
     def _compile_scenario_1_problem(
-        self, extent: npt.NDArray[np.float_]
+        self, extent: npt.NDArray[np.float_], center_frequency: float
     ) -> stride.Problem:
         # scenario constants
         speed_water = 1500  # m/s
-        c_freq = 500e3  # hz
 
         # desired resolution for complexity=fast
         ppw = 6
 
         # compute resolution
-        dx = speed_water / c_freq / ppw  # m
+        dx = speed_water / center_frequency / ppw  # m
 
         grid = make_grid(extent=extent, dx=dx)
         problem = stride.Problem(
@@ -63,7 +62,7 @@ def _compile_scenario_1_problem(
         )
         problem = add_material_fields_to_problem(
             problem=problem,
-            materials=self.get_materials(c_freq),
+            materials=self.get_materials(center_frequency),
             layer_ids=self.layer_ids,
             masks=self._get_material_masks(problem),
         )
@@ -97,9 +96,9 @@ def __init__(self, complexity="fast"):
             complexity=complexity,
         )
 
-    def _compile_problem(self) -> stride.Problem:
+    def _compile_problem(self, center_frequency: float) -> stride.Problem:
         extent = np.array([0.12, 0.07])  # m
-        return self._compile_scenario_1_problem(extent)
+        return self._compile_scenario_1_problem(extent, center_frequency)
 
     def get_default_source(self) -> sources.Source:
         """Return the default source for the scenario."""
@@ -183,9 +182,9 @@ def get_default_slice_position(self, axis: int) -> float:
         default_positions = np.array([0.064, 0.0, 0.0])
         return default_positions[axis]
 
-    def _compile_problem(self) -> stride.Problem:
+    def _compile_problem(self, center_frequency: float) -> stride.Problem:
         extent = np.array([0.12, 0.07, 0.07])  # m
-        return self._compile_scenario_1_problem(extent)
+        return self._compile_scenario_1_problem(extent, center_frequency)
 
     def get_default_source(self) -> sources.Source:
         """Return the default source for the scenario."""
diff --git a/src/neurotechdevkit/scenarios/_scenario_2.py b/src/neurotechdevkit/scenarios/_scenario_2.py
index d43558f7..ee8781b5 100644
--- a/src/neurotechdevkit/scenarios/_scenario_2.py
+++ b/src/neurotechdevkit/scenarios/_scenario_2.py
@@ -41,17 +41,16 @@ def _get_material_masks(self) -> Mapping[str, npt.NDArray[np.bool_]]:
         raise NotImplementedError()
 
     def _compile_scenario_2_problem(
-        self, extent: npt.NDArray[np.float_]
+        self, extent: npt.NDArray[np.float_], center_frequency: float
     ) -> stride.Problem:
         # scenario constants
         speed_water = 1500  # m/s
-        c_freq = 500e3  # hz
 
         # desired resolution for complexity=fast
         ppw = 6
 
         # compute resolution
-        dx = speed_water / c_freq / ppw  # m
+        dx = speed_water / center_frequency / ppw  # m
 
         grid = make_grid(extent=extent, dx=dx)
         problem = stride.Problem(
@@ -59,7 +58,7 @@ def _compile_scenario_2_problem(
         )
         problem = add_material_fields_to_problem(
             problem=problem,
-            materials=self.get_materials(c_freq),
+            materials=self.get_materials(center_frequency),
             layer_ids=self.layer_ids,
             masks=self._get_material_masks(),
         )
@@ -125,9 +124,9 @@ def __init__(self, complexity="fast"):
             complexity=complexity,
         )
 
-    def _compile_problem(self) -> stride.Problem:
+    def _compile_problem(self, center_frequency: float) -> stride.Problem:
         extent = np.array([0.225, 0.170])  # m
-        return self._compile_scenario_2_problem(extent)
+        return self._compile_scenario_2_problem(extent, center_frequency)
 
     def _get_material_masks(self) -> Mapping[str, npt.NDArray[np.bool_]]:
         return {
@@ -256,9 +255,9 @@ def get_default_slice_position(self, axis: int) -> float:
         default_positions = np.array([0.1, 0.0, 0.0])
         return default_positions[axis]
 
-    def _compile_problem(self) -> stride.Problem:
+    def _compile_problem(self, center_frequency: float) -> stride.Problem:
         extent = np.array([0.225, 0.170, 0.190])  # m
-        return self._compile_scenario_2_problem(extent)
+        return self._compile_scenario_2_problem(extent, center_frequency)
 
     def _get_material_masks(self):
         return {
diff --git a/tests/neurotechdevkit/scenarios/test_base.py b/tests/neurotechdevkit/scenarios/test_base.py
index cb1977f5..fc2af5a9 100644
--- a/tests/neurotechdevkit/scenarios/test_base.py
+++ b/tests/neurotechdevkit/scenarios/test_base.py
@@ -41,6 +41,7 @@ class ScenarioBaseTester(Scenario):
 
     def __init__(self, complexity="fast"):
         self._target_id = "target_1"
+        self._problem = self._compile_problem(center_frequency=5e5)
         super().__init__(
             origin=np.array([-0.1, -0.1, 0.0]),
             complexity=complexity,
@@ -50,7 +51,7 @@ def __init__(self, complexity="fast"):
     def _material_outline_upsample_factor(self) -> int:
         return 3
 
-    def _compile_problem(self) -> stride.Problem:
+    def _compile_problem(self, center_frequency: float) -> stride.Problem:
         extent = np.array([2.0, 3.0, 4.0])
         dx = 0.1
         grid = make_grid(extent=extent, dx=dx)
diff --git a/tests/neurotechdevkit/scenarios/test_materials.py b/tests/neurotechdevkit/scenarios/test_materials.py
index 0c706350..b2c03d02 100644
--- a/tests/neurotechdevkit/scenarios/test_materials.py
+++ b/tests/neurotechdevkit/scenarios/test_materials.py
@@ -27,7 +27,7 @@ def __init__(self, complexity="fast"):
             complexity=complexity,
         )
 
-    def _compile_problem(self):
+    def _compile_problem(self, center_frequency: float):
         pass
 
     def _material_outline_upsample_factor(self):
@@ -49,7 +49,7 @@ class ScenarioWithDefaultMaterials(BaseScenario):
     scenario = ScenarioWithDefaultMaterials()
     assert scenario.material_colors == {"brain": "#DB504A", "skin": "#FA8B53"}
 
-    materials = scenario.get_materials()
+    materials = scenario.get_materials(500e3)
     assert list(materials.keys()) == ["brain", "skin"]
     compare_structs(materials["brain"], SUPPORTED_MATERIALS["brain"].to_struct())
     compare_structs(materials["skin"], SUPPORTED_MATERIALS["skin"].to_struct())
@@ -67,7 +67,7 @@ class ScenarioWithCustomMaterialProperties(BaseScenario):
     scenario = ScenarioWithCustomMaterialProperties()
     assert scenario.material_colors == {"brain": "#2E86AB"}
 
-    materials = scenario.get_materials()
+    materials = scenario.get_materials(500e3)
     assert list(materials.keys()) == ["brain"]
     compare_structs(
         materials["brain"],
@@ -87,7 +87,7 @@ class ScenarioWithCustomMaterial(BaseScenario):
     scenario = ScenarioWithCustomMaterial()
     assert scenario.material_colors == {"brain": "#DB504A", "eye": "#2E86AB"}
 
-    materials = scenario.get_materials()
+    materials = scenario.get_materials(500e3)
 
     assert list(materials.keys()) == ["brain", "eye"]
     compare_structs(materials["brain"], SUPPORTED_MATERIALS["brain"].to_struct())
@@ -121,4 +121,4 @@ class ScenarioWithCustomMaterial(BaseScenario):
     with pytest.raises(ValueError):
         scenario.material_colors
     with pytest.raises(ValueError):
-        scenario.get_materials()
+        scenario.get_materials(500e3)

From 72c2105e54f78a9891665a49c333678cf50ceb89 Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Wed, 5 Jul 2023 11:17:11 +0200
Subject: [PATCH 10/17] Renaming variables to make the intention more clear

---
 src/neurotechdevkit/materials.py              | 20 ++++++++++---------
 .../scenarios/test_materials.py               |  8 ++++----
 2 files changed, 15 insertions(+), 13 deletions(-)

diff --git a/src/neurotechdevkit/materials.py b/src/neurotechdevkit/materials.py
index 46ecb4d0..0b741d81 100644
--- a/src/neurotechdevkit/materials.py
+++ b/src/neurotechdevkit/materials.py
@@ -3,6 +3,8 @@
 
 from mosaic.types import Struct
 
+NEPER_TO_DECIBEL = 8.6860000037
+
 
 @dataclass
 class Material:
@@ -28,7 +30,7 @@ def to_struct(self) -> Struct:
 
 
 # The values below consider a center frequency of 500kHz
-SUPPORTED_MATERIALS = {
+DEFAULT_MATERIALS = {
     "water": Material(vp=1500.0, rho=1000.0, alpha=0.0, render_color="#2E86AB"),
     "skin": Material(vp=1610.0, rho=1090.0, alpha=0.2, render_color="#FA8B53"),
     "cortical_bone": Material(vp=2800.0, rho=1850.0, alpha=4.0, render_color="#FAF0CA"),
@@ -66,7 +68,7 @@ def calculate_absorption(self, frequency: float) -> float:
             float: attenuation in dB/cm/MHz
         """
         attenuation = self.a0 * (frequency / 1e6) ** self.b  # Np/m
-        db_attenuation = attenuation * 8.6860000037  # dB/m
+        db_attenuation = attenuation * NEPER_TO_DECIBEL  # dB/m
         return db_attenuation / 100  # convert from dB/m to dB/cm
 
 
@@ -101,7 +103,7 @@ def _calculate_absorption(material_name: str, center_frequency: float) -> float:
             center_frequency
         )
     except KeyError:
-        raise ValueError(f"Unsupported material: {material_name}")
+        raise ValueError(f"Undefined material: {material_name}")
 
 
 def get_render_color(material_name: str) -> str:
@@ -117,10 +119,10 @@ def get_render_color(material_name: str) -> str:
     Returns:
         str: the render color of the material.
     """
-    if material_name not in SUPPORTED_MATERIALS:
-        raise ValueError(f"Unsupported material: {material_name}")
+    if material_name not in DEFAULT_MATERIALS:
+        raise ValueError(f"Undefined material: {material_name}")
 
-    return SUPPORTED_MATERIALS[material_name].render_color
+    return DEFAULT_MATERIALS[material_name].render_color
 
 
 def get_material(material_name: str, center_frequency: float = 5.0e5) -> Material:
@@ -138,10 +140,10 @@ def get_material(material_name: str, center_frequency: float = 5.0e5) -> Materia
     Returns:
         Material: a material with properties used in the neurotechdevkit.
     """
-    if material_name not in SUPPORTED_MATERIALS:
-        raise ValueError(f"Unsupported material: {material_name}")
+    if material_name not in DEFAULT_MATERIALS:
+        raise ValueError(f"Undefined material: {material_name}")
 
-    base_material = SUPPORTED_MATERIALS[material_name]
+    base_material = DEFAULT_MATERIALS[material_name]
     if center_frequency == 5.0e5:
         return base_material
 
diff --git a/tests/neurotechdevkit/scenarios/test_materials.py b/tests/neurotechdevkit/scenarios/test_materials.py
index b2c03d02..80dda432 100644
--- a/tests/neurotechdevkit/scenarios/test_materials.py
+++ b/tests/neurotechdevkit/scenarios/test_materials.py
@@ -2,7 +2,7 @@
 import pytest
 from mosaic.types import Struct
 
-from neurotechdevkit.materials import SUPPORTED_MATERIALS, Material
+from neurotechdevkit.materials import DEFAULT_MATERIALS, Material
 from neurotechdevkit.scenarios import Scenario
 
 
@@ -51,8 +51,8 @@ class ScenarioWithDefaultMaterials(BaseScenario):
 
     materials = scenario.get_materials(500e3)
     assert list(materials.keys()) == ["brain", "skin"]
-    compare_structs(materials["brain"], SUPPORTED_MATERIALS["brain"].to_struct())
-    compare_structs(materials["skin"], SUPPORTED_MATERIALS["skin"].to_struct())
+    compare_structs(materials["brain"], DEFAULT_MATERIALS["brain"].to_struct())
+    compare_structs(materials["skin"], DEFAULT_MATERIALS["skin"].to_struct())
 
 
 def test_custom_material_property():
@@ -90,7 +90,7 @@ class ScenarioWithCustomMaterial(BaseScenario):
     materials = scenario.get_materials(500e3)
 
     assert list(materials.keys()) == ["brain", "eye"]
-    compare_structs(materials["brain"], SUPPORTED_MATERIALS["brain"].to_struct())
+    compare_structs(materials["brain"], DEFAULT_MATERIALS["brain"].to_struct())
     compare_structs(
         materials["eye"],
         Material(vp=1600.0, rho=1100.0, alpha=0.0, render_color="#2E86AB").to_struct(),

From 6693fe74109274d4b90e40e6f174698335ff8b4d Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Wed, 5 Jul 2023 11:29:11 +0200
Subject: [PATCH 11/17] Fixing spellcheck error

---
 src/neurotechdevkit/scenarios/_base.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/neurotechdevkit/scenarios/_base.py b/src/neurotechdevkit/scenarios/_base.py
index 4110dcca..3949482c 100644
--- a/src/neurotechdevkit/scenarios/_base.py
+++ b/src/neurotechdevkit/scenarios/_base.py
@@ -64,7 +64,7 @@ class Scenario(abc.ABC):
     # The list of material layers in the scenario.
     material_layers: list[str]
 
-    # The customizations to the material layers.
+    # The customization to the material layers.
     material_properties: dict[str, Material] = {}
 
     def __init__(

From a6bf5474f18e9c1263b4f9e02d461311977b5d7f Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Wed, 5 Jul 2023 13:48:10 +0200
Subject: [PATCH 12/17] Calculating material alpha for all frequencies

---
 src/neurotechdevkit/materials.py              | 26 ++++++++++---------
 src/neurotechdevkit/rendering/napari.py       |  1 -
 src/neurotechdevkit/results/_results.py       |  1 +
 .../scenarios/test_materials.py               | 18 +------------
 4 files changed, 16 insertions(+), 30 deletions(-)

diff --git a/src/neurotechdevkit/materials.py b/src/neurotechdevkit/materials.py
index 0b741d81..805e74ec 100644
--- a/src/neurotechdevkit/materials.py
+++ b/src/neurotechdevkit/materials.py
@@ -7,14 +7,20 @@
 
 
 @dataclass
-class Material:
+class _BaseMaterial:
     """A material with properties used in the neurotechdevkit scenarios."""
 
     vp: float
     rho: float
-    alpha: float
     render_color: str
 
+
+@dataclass
+class Material(_BaseMaterial):
+    """A NDK Material with an attenuation coefficient."""
+
+    alpha: float
+
     def to_struct(self) -> Struct:
         """Return a Struct representation of the material.
 
@@ -31,16 +37,14 @@ def to_struct(self) -> Struct:
 
 # The values below consider a center frequency of 500kHz
 DEFAULT_MATERIALS = {
-    "water": Material(vp=1500.0, rho=1000.0, alpha=0.0, render_color="#2E86AB"),
-    "skin": Material(vp=1610.0, rho=1090.0, alpha=0.2, render_color="#FA8B53"),
-    "cortical_bone": Material(vp=2800.0, rho=1850.0, alpha=4.0, render_color="#FAF0CA"),
-    "trabecular_bone": Material(
-        vp=2300.0, rho=1700.0, alpha=8.0, render_color="#EBD378"
-    ),
-    "brain": Material(vp=1560.0, rho=1040.0, alpha=0.3, render_color="#DB504A"),
+    "water": _BaseMaterial(vp=1500.0, rho=1000.0, render_color="#2E86AB"),
+    "skin": _BaseMaterial(vp=1610.0, rho=1090.0, render_color="#FA8B53"),
+    "cortical_bone": _BaseMaterial(vp=2800.0, rho=1850.0, render_color="#FAF0CA"),
+    "trabecular_bone": _BaseMaterial(vp=2300.0, rho=1700.0, render_color="#EBD378"),
+    "brain": _BaseMaterial(vp=1560.0, rho=1040.0, render_color="#DB504A"),
     # these numbers are completely made up
     # TODO: research reasonable values
-    "tumor": Material(vp=1650.0, rho=1150.0, alpha=0.8, render_color="#94332F"),
+    "tumor": _BaseMaterial(vp=1650.0, rho=1150.0, render_color="#94332F"),
 }
 
 
@@ -144,8 +148,6 @@ def get_material(material_name: str, center_frequency: float = 5.0e5) -> Materia
         raise ValueError(f"Undefined material: {material_name}")
 
     base_material = DEFAULT_MATERIALS[material_name]
-    if center_frequency == 5.0e5:
-        return base_material
 
     return Material(
         vp=base_material.vp,
diff --git a/src/neurotechdevkit/rendering/napari.py b/src/neurotechdevkit/rendering/napari.py
index 2d049506..8c7af3b4 100644
--- a/src/neurotechdevkit/rendering/napari.py
+++ b/src/neurotechdevkit/rendering/napari.py
@@ -98,7 +98,6 @@ def render_amplitudes_3d_with_napari(result: "results.SteadyStateResult3D") -> N
     Raises:
         ImportError: If napari is not found.
     """
-    pass
     _create_napari_3d(scenario=result.scenario, amplitudes=result.get_steady_state())
 
 
diff --git a/src/neurotechdevkit/results/_results.py b/src/neurotechdevkit/results/_results.py
index 0f091e6b..b0a5b574 100644
--- a/src/neurotechdevkit/results/_results.py
+++ b/src/neurotechdevkit/results/_results.py
@@ -1297,6 +1297,7 @@ def load_result_from_disk(filepath: str | pathlib.Path) -> Result:
             wavefield=save_data.get("wavefield"),
             traces=None,
         )
+        scenario._problem = scenario._compile_problem(save_data["center_frequency"])
 
         if save_data.get("steady_state") is not None:
             fields_kwargs.update(steady_state=save_data["steady_state"])
diff --git a/tests/neurotechdevkit/scenarios/test_materials.py b/tests/neurotechdevkit/scenarios/test_materials.py
index 80dda432..e4728a6d 100644
--- a/tests/neurotechdevkit/scenarios/test_materials.py
+++ b/tests/neurotechdevkit/scenarios/test_materials.py
@@ -2,7 +2,7 @@
 import pytest
 from mosaic.types import Struct
 
-from neurotechdevkit.materials import DEFAULT_MATERIALS, Material
+from neurotechdevkit.materials import Material
 from neurotechdevkit.scenarios import Scenario
 
 
@@ -40,21 +40,6 @@ def get_target_mask(self):
         pass
 
 
-def test_default_materials():
-    """Test that the default materials are used."""
-
-    class ScenarioWithDefaultMaterials(BaseScenario):
-        material_layers = ["brain", "skin"]
-
-    scenario = ScenarioWithDefaultMaterials()
-    assert scenario.material_colors == {"brain": "#DB504A", "skin": "#FA8B53"}
-
-    materials = scenario.get_materials(500e3)
-    assert list(materials.keys()) == ["brain", "skin"]
-    compare_structs(materials["brain"], DEFAULT_MATERIALS["brain"].to_struct())
-    compare_structs(materials["skin"], DEFAULT_MATERIALS["skin"].to_struct())
-
-
 def test_custom_material_property():
     """Test that a custom material property is used."""
 
@@ -90,7 +75,6 @@ class ScenarioWithCustomMaterial(BaseScenario):
     materials = scenario.get_materials(500e3)
 
     assert list(materials.keys()) == ["brain", "eye"]
-    compare_structs(materials["brain"], DEFAULT_MATERIALS["brain"].to_struct())
     compare_structs(
         materials["eye"],
         Material(vp=1600.0, rho=1100.0, alpha=0.0, render_color="#2E86AB").to_struct(),

From 7dff8c41a8a40400328b94f94af5d48a6f4c98a1 Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Wed, 5 Jul 2023 13:57:54 +0200
Subject: [PATCH 13/17] Defining fixed material properties to the scenarios to
 match the paper

---
 src/neurotechdevkit/scenarios/_scenario_0.py |  9 +++++++++
 src/neurotechdevkit/scenarios/_scenario_1.py | 12 ++++++++++++
 src/neurotechdevkit/scenarios/_scenario_2.py |  8 ++++++++
 3 files changed, 29 insertions(+)

diff --git a/src/neurotechdevkit/scenarios/_scenario_0.py b/src/neurotechdevkit/scenarios/_scenario_0.py
index 29c419c7..ad0cfc9e 100644
--- a/src/neurotechdevkit/scenarios/_scenario_0.py
+++ b/src/neurotechdevkit/scenarios/_scenario_0.py
@@ -1,6 +1,7 @@
 import numpy as np
 import stride
 
+from ..materials import Material
 from ..sources import FocusedSource2D
 from ._base import Scenario2D, Target
 from ._utils import (
@@ -29,6 +30,14 @@ class Scenario0(Scenario2D):
         "brain",
         "tumor",
     ]
+    material_properties = {
+        "water": Material(vp=1500.0, rho=1000.0, alpha=0.0, render_color="#2E86AB"),
+        "cortical_bone": Material(
+            vp=2800.0, rho=1850.0, alpha=4.0, render_color="#FAF0CA"
+        ),
+        "brain": Material(vp=1560.0, rho=1040.0, alpha=0.3, render_color="#DB504A"),
+        "tumor": Material(vp=1650.0, rho=1150.0, alpha=0.8, render_color="#94332F"),
+    }
 
     def __init__(self, complexity="fast"):
         """Create a new instance of scenario 0."""
diff --git a/src/neurotechdevkit/scenarios/_scenario_1.py b/src/neurotechdevkit/scenarios/_scenario_1.py
index 87125000..b1fd15ff 100644
--- a/src/neurotechdevkit/scenarios/_scenario_1.py
+++ b/src/neurotechdevkit/scenarios/_scenario_1.py
@@ -7,6 +7,7 @@
 import stride
 
 from .. import rendering, sources
+from ..materials import Material
 from ._base import Scenario, Scenario2D, Scenario3D, Target
 from ._utils import add_material_fields_to_problem, make_grid
 
@@ -31,6 +32,17 @@ class Scenario1(Scenario):
         "trabecular_bone",
         "brain",
     ]
+    material_properties = {
+        "water": Material(vp=1500.0, rho=1000.0, alpha=0.0, render_color="#2E86AB"),
+        "skin": Material(vp=1610.0, rho=1090.0, alpha=0.2, render_color="#FA8B53"),
+        "cortical_bone": Material(
+            vp=2800.0, rho=1850.0, alpha=4.0, render_color="#FAF0CA"
+        ),
+        "trabecular_bone": Material(
+            vp=2300.0, rho=1700.0, alpha=8.0, render_color="#EBD378"
+        ),
+        "brain": Material(vp=1560.0, rho=1040.0, alpha=0.3, render_color="#DB504A"),
+    }
 
     @property
     def _material_outline_upsample_factor(self) -> int:
diff --git a/src/neurotechdevkit/scenarios/_scenario_2.py b/src/neurotechdevkit/scenarios/_scenario_2.py
index ee8781b5..647f81d9 100644
--- a/src/neurotechdevkit/scenarios/_scenario_2.py
+++ b/src/neurotechdevkit/scenarios/_scenario_2.py
@@ -9,6 +9,7 @@
 import stride
 
 from .. import rendering, sources
+from ..materials import Material
 from ._base import Scenario, Scenario2D, Scenario3D, Target
 from ._utils import add_material_fields_to_problem, make_grid
 
@@ -31,6 +32,13 @@ class Scenario2(Scenario):
         "cortical_bone",
         "brain",
     ]
+    material_properties = {
+        "water": Material(vp=1500.0, rho=1000.0, alpha=0.0, render_color="#2E86AB"),
+        "cortical_bone": Material(
+            vp=2800.0, rho=1850.0, alpha=4.0, render_color="#FAF0CA"
+        ),
+        "brain": Material(vp=1560.0, rho=1040.0, alpha=0.3, render_color="#DB504A"),
+    }
 
     @property
     def _material_outline_upsample_factor(self) -> int:

From 242fe7d9a34d884d8e8e51e5d7322f8ea33a1d90 Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Wed, 5 Jul 2023 14:03:33 +0200
Subject: [PATCH 14/17] Linting gallery examples

---
 Makefile                                      |  20 +--
 docs/examples/plot_3d.py                      |  29 ++--
 docs/examples/plot_adding_custom_source.py    |   7 +-
 .../plot_customized_center_frequency.py       |   7 +-
 docs/examples/plot_full_scenario.py           |   7 +-
 docs/examples/plot_metrics.py                 |  11 +-
 docs/examples/plot_multiple_sources.py        |  34 ++--
 docs/examples/plot_phased_array_source.py     | 160 +++++++++++-------
 docs/examples/plot_pulsed_simulation.py       |   7 +-
 docs/examples/plot_scenarios.py               |  15 +-
 docs/examples/plot_store_results.py           |  17 +-
 11 files changed, 191 insertions(+), 123 deletions(-)

diff --git a/Makefile b/Makefile
index 41778a41..5910cdca 100644
--- a/Makefile
+++ b/Makefile
@@ -4,11 +4,11 @@ help:
 	@echo "Available commands are: \n*lint, lint-check, spellcheck, test, test-unit, test-integration docs"
 
 lint:
-	poetry run isort src tests
-	poetry run black src tests
-	poetry run flake8 src tests
-	poetry run mypy src
-	poetry run codespell src
+	poetry run isort src tests docs/examples
+	poetry run black src tests docs/examples
+	poetry run flake8 src tests docs/examples
+	poetry run mypy src docs/examples
+	poetry run codespell src docs/examples
 	poetry run pydocstyle src
 	poetry run pyright
 
@@ -16,11 +16,11 @@ spellcheck:
 	poetry run pylint --disable all --enable spelling --spelling-dict en_US --spelling-private-dict-file=whitelist.txt src
 
 lint-check:
-	poetry run isort --check src tests
-	poetry run black --check src tests
-	poetry run flake8  src tests
-	poetry run mypy src
-	poetry run codespell src
+	poetry run isort --check src tests docs/examples
+	poetry run black --check src tests docs/examples
+	poetry run flake8  src tests docs/examples
+	poetry run mypy src docs/examples
+	poetry run codespell src docs/examples
 	poetry run pydocstyle src
 	poetry run pyright --warnings
 
diff --git a/docs/examples/plot_3d.py b/docs/examples/plot_3d.py
index 5289d2b1..a90bba62 100644
--- a/docs/examples/plot_3d.py
+++ b/docs/examples/plot_3d.py
@@ -3,7 +3,8 @@
 Visualizing 3D results with Napari
 ====================================
 
-This example demonstrates how to render a steady state result for a 3D scenario using [napari](https://napari.org/stable/).
+This example demonstrates how to render a steady state result for
+a 3D scenario using [napari](https://napari.org/stable/).
 
 Running such simulations is computationally expensive and can take a long time to
 complete. For this reason, we recommend running this simulation on an external machine,
@@ -15,13 +16,12 @@
 # %%
 # The following step downloads and loads a simulation executed on an external machine.
 import pooch
+
 import neurotechdevkit as ndk
 
-URL = 'https://neurotechdevkit.s3.us-west-2.amazonaws.com/result-scenario-2-3d.tz'
+URL = "https://neurotechdevkit.s3.us-west-2.amazonaws.com/result-scenario-2-3d.tz"
 known_hash = "6a5de26466028c673d253ca014c75c719467ec6c28d7178baf9287b44ad15191"
-downloaded_file_path = pooch.retrieve(
-    url=URL, known_hash=known_hash, progressbar=True
-)
+downloaded_file_path = pooch.retrieve(url=URL, known_hash=known_hash, progressbar=True)
 result = ndk.load_result_from_disk(downloaded_file_path)
 
 # %%
@@ -32,21 +32,28 @@
 # pip install "napari[all]"
 # ```
 #
-# or by following the `napari` [installation instructions](https://napari.org/stable/tutorials/fundamentals/installation.html).
+# or by following the `napari` installation instructions
+# [link](https://napari.org/stable/tutorials/fundamentals/installation.html).
 
 try:
-    import napari
+    import napari  # noqa: F401
+
+    assert isinstance(result, ndk.results.SteadyStateResult3D)
     result.render_steady_state_amplitudes_3d()
 except ImportError:
-    print("napari has not been installed. Please install it with: pip install napari[all]")
+    print(
+        "napari has not been installed. Please install it with: pip install napari[all]"
+    )
 
 # %%
 # If you have napari installed you should see an output like the following:
 #
 # ```
-# Opening the napari viewer. The window might not show up on top of your notebook; look through your open applications if it does not."""
+# Opening the napari viewer. The window might not show up on top of your notebook;
+# look through your open applications if it does not.
 # ```
 
 # %%
-# If you have napari installed you should have been able to see an image like the following:
-# ![3d-visualization](../../images/3d_visualization.gif)
\ No newline at end of file
+# If you have napari installed you should have been able to see an image like
+# the following:
+# ![3d-visualization](../../images/3d_visualization.gif)
diff --git a/docs/examples/plot_adding_custom_source.py b/docs/examples/plot_adding_custom_source.py
index 22615213..39a7bd1e 100644
--- a/docs/examples/plot_adding_custom_source.py
+++ b/docs/examples/plot_adding_custom_source.py
@@ -5,7 +5,8 @@
 ====================================
 
 !!! note
-    NDK and its examples are under constant development, more information and content will be added to this example soon!
+    NDK and its examples are under constant development, more information
+    and content will be added to this example soon!
 
 This example demonstrates how to add a source to the simulation.
 """
@@ -27,6 +28,7 @@
 #   emitting. Defaults to 0.0.
 
 import numpy as np
+
 import neurotechdevkit as ndk
 
 source = ndk.sources.FocusedSource2D(
@@ -37,9 +39,10 @@
     num_points=1000,
 )
 
-scenario = ndk.make('scenario-0-v0')
+scenario = ndk.make("scenario-0-v0")
 scenario.add_source(source)
 result = scenario.simulate_steady_state()
+assert isinstance(result, ndk.results.SteadyStateResult2D)
 result.render_steady_state_amplitudes()
 
 # %%
diff --git a/docs/examples/plot_customized_center_frequency.py b/docs/examples/plot_customized_center_frequency.py
index 42ead974..5901ba86 100644
--- a/docs/examples/plot_customized_center_frequency.py
+++ b/docs/examples/plot_customized_center_frequency.py
@@ -11,17 +11,20 @@
 
 CENTER_FREQUENCY = 6e5
 
-scenario = ndk.make('scenario-0-v0')
+scenario = ndk.make("scenario-0-v0")
 
 # using default material layers
 scenario.material_layers = ["water", "cortical_bone", "brain", "tumor"]
 
 # Customizing material properties with random values
 scenario.material_properties = {
-  "tumor": ndk.materials.Material(vp=1850.0, rho=1250.0,alpha=0.8,render_color="#94332F"),
+    "tumor": ndk.materials.Material(
+        vp=1850.0, rho=1250.0, alpha=0.8, render_color="#94332F"
+    ),
 }
 
 result = scenario.simulate_steady_state(center_frequency=CENTER_FREQUENCY)
+assert isinstance(result, ndk.results.SteadyStateResult2D)
 result.render_steady_state_amplitudes(show_material_outlines=False)
 
 # %%
diff --git a/docs/examples/plot_full_scenario.py b/docs/examples/plot_full_scenario.py
index 4f07be3c..6cc4e0d3 100644
--- a/docs/examples/plot_full_scenario.py
+++ b/docs/examples/plot_full_scenario.py
@@ -4,7 +4,8 @@
 ====================================
 
 !!! note
-    NDK and its examples are under constant development, more information and content will be added to this example soon!
+    NDK and its examples are under constant development, more information
+    and content will be added to this example soon!
 
 The following code is a simplified implementation of NDK's Scenario 1.
 """
@@ -14,11 +15,12 @@
 import stride
 
 from neurotechdevkit import sources
+from neurotechdevkit.results import SteadyStateResult2D
 from neurotechdevkit.scenarios import (
     Scenario2D,
     Target,
-    make_grid,
     add_material_fields_to_problem,
+    make_grid,
 )
 
 
@@ -171,4 +173,5 @@ def _fill_mask(mask, start, end, dx):
 
 scenario = FullScenario()
 result = scenario.simulate_steady_state()
+assert isinstance(result, SteadyStateResult2D)
 result.render_steady_state_amplitudes(show_material_outlines=False)
diff --git a/docs/examples/plot_metrics.py b/docs/examples/plot_metrics.py
index 7682a727..7ccc4481 100644
--- a/docs/examples/plot_metrics.py
+++ b/docs/examples/plot_metrics.py
@@ -5,7 +5,8 @@
 ====================================
 
 !!! note
-    NDK and its examples are under constant development, more information and content will be added to this example soon!
+    NDK and its examples are under constant development, more information
+    and content will be added to this example soon!
 
 This example demonstrates how to display the metrics collected from the simulation.
 """
@@ -13,9 +14,9 @@
 # ## Rendering scenario
 import neurotechdevkit as ndk
 
-
-scenario = ndk.make('scenario-0-v0')
+scenario = ndk.make("scenario-0-v0")
 result = scenario.simulate_steady_state()
+assert isinstance(result, ndk.results.SteadyStateResult2D)
 result.render_steady_state_amplitudes(show_material_outlines=False)
 
 # %%
@@ -23,5 +24,7 @@
 for metric, metric_value in result.metrics.items():
     print(f"{metric}:")
     print(f"\t {metric_value['description']}")
-    print(f"\t Unit: {metric_value['unit-of-measurement']} Value: {metric_value['value']}")
+    print(
+        f"\t Unit: {metric_value['unit-of-measurement']} Value: {metric_value['value']}"
+    )
     print()
diff --git a/docs/examples/plot_multiple_sources.py b/docs/examples/plot_multiple_sources.py
index 3bd074af..f4a0def1 100644
--- a/docs/examples/plot_multiple_sources.py
+++ b/docs/examples/plot_multiple_sources.py
@@ -4,7 +4,8 @@
 ====================================
 
 !!! note
-    NDK and its examples are under constant development, more information and content will be added to this example soon!
+    NDK and its examples are under constant development, more information
+    and content will be added to this example soon!
 
 Adding multiple sources for transcranial ultrasound stimulation enables
 greater precision and control in targeting specific areas of the brain.
@@ -15,31 +16,32 @@
 """
 # %%
 import numpy as np
+
 import neurotechdevkit as ndk
 
-scenario = ndk.make('scenario-0-v0')
+scenario = ndk.make("scenario-0-v0")
 
 s1 = ndk.sources.FocusedSource2D(
-            position=np.array([0.01, 0.0]),
-            direction=np.array([0.92, 0.25]),
-            aperture=0.01,
-            focal_length=0.022,
-            num_points=1000,
-        )
+    position=np.array([0.01, 0.0]),
+    direction=np.array([0.92, 0.25]),
+    aperture=0.01,
+    focal_length=0.022,
+    num_points=1000,
+)
 
 s2 = ndk.sources.FocusedSource2D(
-            position=np.array([0.04, -0.002]),
-            direction=np.array([-0.85, 0.35]),
-            aperture=0.01,
-            focal_length=0.011,
-            num_points=1000,
-            delay=5.1e-6,
-        )
+    position=np.array([0.04, -0.002]),
+    direction=np.array([-0.85, 0.35]),
+    aperture=0.01,
+    focal_length=0.011,
+    num_points=1000,
+    delay=5.1e-6,
+)
 
 scenario.add_source(s1)
 scenario.add_source(s2)
 result = scenario.simulate_pulse()
-
+assert isinstance(result, ndk.results.PulsedResult2D)
 result.render_pulsed_simulation_animation()
 
 # %%
diff --git a/docs/examples/plot_phased_array_source.py b/docs/examples/plot_phased_array_source.py
index acc9ce8b..ac44d3eb 100644
--- a/docs/examples/plot_phased_array_source.py
+++ b/docs/examples/plot_phased_array_source.py
@@ -3,26 +3,55 @@
 Phased array source
 ====================================
 
-Phased arrays consist of multiple individual transducer elements arranged in a specific pattern, such as a linear or circle. Each element can be controlled independently, allowing for precise manipulation of the generated ultrasound waves. By adjusting the timing or phase of the signals sent to these elements, the ultrasound waves can be focused, steered, and shaped without mechanically moving the transducer.
-
-They are becoming increasingly popular in the field of transcranial ultrasound stimulation as they offer several advantages over traditional ultrasound transducers -- among others are precise targeting and better control over steering and shaping:
-
-* **Precise targeting:** as individual transducer elements can be individually controlled, it allows for the generation of a focused ultrasound beam with high spatial accuracy. This enables the stimulation of specific brain regions without affecting the surrounding healthy tissue and minimizes the risk of potential side effects.
-
-* **Steering and shaping:** The phased array technology allows the ultrasound beam to be electronically steered and shaped in real-time without mechanically moving the transducers. This enables targeting of different brain regions or adjusting the stimulation pattern as needed during a session, making the procedure more versatile and adaptable.
-
-These features allow the stimulation to be tailored to suit to more specific requirements. With ongoing research and development, they have the potential to revolutionize the field of brain stimulation and offer new treatment options for a range of neurological and psychiatric disorders.
-
-This notebook will show how to define a phased array within NDK and experiment with some of the available features. For more details, checkout the [NDK documentation](https://agencyenterprise.github.io/neurotechdevkit/).
+Phased arrays consist of multiple individual transducer elements arranged in a
+specific pattern, such as a linear or circle. Each element can be controlled
+independently, allowing for precise manipulation of the generated ultrasound
+waves. By adjusting the timing or phase of the signals sent to these elements,
+the ultrasound waves can be focused, steered, and shaped without mechanically
+moving the transducer.
+
+They are becoming increasingly popular in the field of transcranial ultrasound
+stimulation as they offer several advantages over traditional ultrasound
+transducers -- among others are precise targeting and better control over steering
+and shaping:
+
+* **Precise targeting:** as individual transducer elements can be individually
+controlled, it allows for the generation of a focused ultrasound beam with high
+spatial accuracy. This enables the stimulation of specific brain regions without
+affecting the surrounding healthy tissue and minimizes the risk of potential side
+effects.
+
+* **Steering and shaping:** The phased array technology allows the ultrasound beam
+to be electronically steered and shaped in real-time without mechanically moving the
+transducers. This enables targeting of different brain regions or adjusting the
+stimulation pattern as needed during a session, making the procedure more versatile
+and adaptable.
+
+These features allow the stimulation to be tailored to suit to more specific
+requirements. With ongoing research and development, they have the potential to
+revolutionize the field of brain stimulation and offer new treatment options for
+a range of neurological and psychiatric disorders.
+
+This notebook will show how to define a phased array within NDK and experiment with
+some of the available features. For more details, checkout the
+[NDK documentation](https://agencyenterprise.github.io/neurotechdevkit/).
 
 ## Phased Arrays Components
 
-* **Elements:** A phased array consists of multiple transducer elements that can be individually controlled in terms of phase.
-* **Pitch:** The pitch of a phased array refers to the distance between adjacent transducer elements, and affects the resolution of the array.
-* **Spacing:** The spacing between transducer elements affects the width and shape of the ultrasound beam that is produced.
-* **Aperture:** The aperture of a phased array refers to the total size of the array, and affects the power and focus of the ultrasound beam.
-* **Element width:** The element width of a phased array refers to the size of each individual transducer element. Wider elements generally result in a more powerful and focused beam, while narrower elements can increase the resolution of the array.
-* **Tilt angle**: The tilt angle of a phased array refers to the angle at which the ultrasound beam is directed relative to the normal of plane of the array. Tilt angle can be adjusted by controlling the phase (delays) of the individual transducer elements.
+* **Elements:** A phased array consists of multiple transducer elements that can be
+individually controlled in terms of phase.
+* **Pitch:** The pitch of a phased array refers to the distance between adjacent
+transducer elements, and affects the resolution of the array.
+* **Spacing:** The spacing between transducer elements affects the width and shape
+of the ultrasound beam that is produced.
+* **Aperture:** The aperture of a phased array refers to the total size of the array,
+and affects the power and focus of the ultrasound beam.
+* **Element width:** The element width of a phased array refers to the size of each
+individual transducer element. Wider elements generally result in a more powerful
+and focused beam, while narrower elements can increase the resolution of the array.
+* **Tilt angle**: The tilt angle of a phased array refers to the angle at which the
+ultrasound beam is directed relative to the normal of plane of the array. Tilt angle
+can be adjusted by controlling the phase (delays) of the individual transducer elements.
 
 ![phased-array-schematic](../../images/phased-array-schematic.png)
 
@@ -31,10 +60,13 @@
 
 ![phased-array-axis-definition](../../images/phased-array-axis-definition.png)
 
-The image above shows the axis definition and the look of a phased array in the real world. Image source [link](https://www.biomecardio.com/MUST/functions/html/txdelay_doc.html).
+The image above shows the axis definition and the look of a phased array in the real
+world. Image source
+[link](https://www.biomecardio.com/MUST/functions/html/txdelay_doc.html).
 
 
-Below we will go through some examples of how to use the NDK to API to simulate Phased Arrays in the most common settings. Examples we wil cover:
+Below we will go through some examples of how to use the NDK to API to simulate
+Phased Arrays in the most common settings. Examples we will cover:
 
 * Setting up a tilted wavefront
 * Focusing phased arrays
@@ -58,10 +90,11 @@
 - num_elements `(int)`: the number of elements of the phased array.
 - pitch `(float)`: the distance (in meters) between the centers of neighboring
     elements in the phased array.
-- element_width `(float)`: the width (in meters) of each individual element of the array.
-- tilt_angle `(float)`: the desired tilt angle (in degrees) of the wavefront. The angle is
-    measured between the direction the wavefront travels and the normal to the
-    surface of the transducer, with positive angles resulting in a
+- element_width `(float)`: the width (in meters) of each individual element of
+    the array.
+- tilt_angle `(float)`: the desired tilt angle (in degrees) of the wavefront.
+    The angle is measured between the direction the wavefront travels and the
+    normal to the surface of the transducer, with positive angles resulting in a
     counter-clockwise tilt away from the normal.
 - focal_length `(float)`: the distance (in meters) from `position` to the focal
     point.
@@ -81,75 +114,80 @@
 # argument different than zero. Positive angles result in counter-clockwise rotations.
 
 import numpy as np
-import neurotechdevkit as ndk
 
+import neurotechdevkit as ndk
 
 # define the source
 source = ndk.sources.PhasedArraySource2D(
-            position=np.array([0.0, 0.0]),
-            direction=np.array([1., 0.]),
-            num_elements=20,
-            pitch=1.5e-3,
-            element_width=1.2e-3,
-            tilt_angle=30,
-            num_points=1000,
+    position=np.array([0.0, 0.0]),
+    direction=np.array([1.0, 0.0]),
+    num_elements=20,
+    pitch=1.5e-3,
+    element_width=1.2e-3,
+    tilt_angle=30,
+    num_points=1000,
 )
 
-scenario = ndk.make('scenario-1-2d-v0')
+scenario = ndk.make("scenario-1-2d-v0")
 scenario.add_source(source)
 result = scenario.simulate_steady_state()
+assert isinstance(result, ndk.results.SteadyStateResult2D)
 result.render_steady_state_amplitudes()
 
 # %%
 # ## Focusing
-# Similarly, a phased array can be focused by applying delays in a specific way. Such delays are automatically computed when `focal_length` is defined. Let's explore how the API looks like:
+# Similarly, a phased array can be focused by applying delays in a specific way. Such
+# delays are automatically computed when `focal_length` is defined. Let's explore how
+# the API looks like:
 
-scenario = ndk.make('scenario-1-2d-v0')
+scenario = ndk.make("scenario-1-2d-v0")
 
 phased_array = ndk.sources.PhasedArraySource2D(
-            position=np.array([0.0, 0.0]),
-            direction=np.array([1.0, 0.0]),
-            num_elements=20,
-            pitch=1.5e-3,
-            element_width=1.2e-3,
-            focal_length=0.03,
-            num_points=1000,
+    position=np.array([0.0, 0.0]),
+    direction=np.array([1.0, 0.0]),
+    num_elements=20,
+    pitch=1.5e-3,
+    element_width=1.2e-3,
+    focal_length=0.03,
+    num_points=1000,
 )
 
 scenario.add_source(phased_array)
-print(f'Focal point is: {phased_array.focal_point}')
+print(f"Focal point is: {phased_array.focal_point}")
 # %%
 # `focal_point` shows the coordinates (in meters) where the array focuses.
 
 result = scenario.simulate_steady_state()
+assert isinstance(result, ndk.results.SteadyStateResult2D)
 result.render_steady_state_amplitudes()
 # %%
 # We can see that the arrays focuses at a distance equal to `focal_length`.
 
 
-
 # %%
 # ## Custom Wavefront Shapes
-# In the examples shown so far we specified high level arguments and the required delays were
-# automatically computed. The NDK API offers also the possibility to specify delays directly
-# to create arbitrary wavefront shapes.
+# In the examples shown so far we specified high level arguments and the required
+# delays were automatically computed. The NDK API offers also the possibility to
+# specify delays directly to create arbitrary wavefront shapes.
 #
-# In the example below we apply monotonically increasing delays to mimic a counter-clockwise angle.
+# In the example below we apply monotonically increasing delays to mimic a
+# counter-clockwise angle.
 
-scenario = ndk.make('scenario-1-2d-v0')
+scenario = ndk.make("scenario-1-2d-v0")
 
 phased_array = ndk.sources.PhasedArraySource2D(
-            position=np.array([0.0, 0.0]),
-            direction=np.array([1., 0.]),
-            num_elements=20,
-            pitch=1.5e-3,
-            element_width=1.2e-3,
-            element_delays=np.linspace(0, 1e-5, 20),
-            num_points=1000,
+    position=np.array([0.0, 0.0]),
+    direction=np.array([1.0, 0.0]),
+    num_elements=20,
+    pitch=1.5e-3,
+    element_width=1.2e-3,
+    element_delays=np.linspace(0, 1e-5, 20),
+    num_points=1000,
 )
 
 scenario.add_source(phased_array)
 result = scenario.simulate_steady_state()
+assert isinstance(result, ndk.results.SteadyStateResult2D)
 result.render_steady_state_amplitudes()
 
 # %%
@@ -166,21 +204,23 @@
 #
 # The rest of the API is identical for both 2D and 3D scenarios.
 
-scenario = ndk.make('scenario-1-3d-v0')
+scenario = ndk.make("scenario-1-3d-v0")
 
 phased_3d = ndk.sources.PhasedArraySource3D(
-    position=np.array([0.,0.,0.]),
-    direction=np.array([1.,0.,0.]),
-    center_line=np.array([0.,0.,1.]),
+    position=np.array([0.0, 0.0, 0.0]),
+    direction=np.array([1.0, 0.0, 0.0]),
+    center_line=np.array([0.0, 0.0, 1.0]),
     num_points=20_000,
     num_elements=16,
     pitch=1.5e-3,
     element_width=1.2e-3,
     tilt_angle=30,
-    height=5.e-3)
+    height=5.0e-3,
+)
 
 scenario.add_source(phased_3d)
 results = scenario.simulate_steady_state()
+assert isinstance(results, ndk.results.SteadyStateResult3D)
 results.render_steady_state_amplitudes(slice_axis=1, slice_position=0.0)
 
 # %%
diff --git a/docs/examples/plot_pulsed_simulation.py b/docs/examples/plot_pulsed_simulation.py
index 067d2f59..dee1c005 100644
--- a/docs/examples/plot_pulsed_simulation.py
+++ b/docs/examples/plot_pulsed_simulation.py
@@ -11,13 +11,16 @@
 
 scenario = ndk.make("scenario-0-v0")
 result = scenario.simulate_pulse()
+assert isinstance(result, ndk.results.PulsedResult2D)
 result.render_pulsed_simulation_animation()
 
 # %%
 #
 # Generating a video file
 # =======================
-# You can also generate a video file of the simulation (which requires [ffmpeg](https://ffmpeg.org/download.html) installed).
+# You can also generate a video file of the simulation (which requires
+# [ffmpeg](https://ffmpeg.org/download.html) installed).
 #
-# To create and save the video as `animation.mp4` in the current folder, all you need is the execute following command:
+# To create and save the video as `animation.mp4` in the current folder, all you
+# need is the execute following command:
 result.create_video_file("animation.mp4", fps=25, overwrite=True)
diff --git a/docs/examples/plot_scenarios.py b/docs/examples/plot_scenarios.py
index 1711704c..fcdbadb5 100644
--- a/docs/examples/plot_scenarios.py
+++ b/docs/examples/plot_scenarios.py
@@ -11,22 +11,23 @@
 
 
 def plot_scenario(scenario_id):
-  print(f"Simulating scenario: {scenario_id}")
-  scenario = ndk.make(scenario_id)
-  result = scenario.simulate_steady_state()
-  result.render_steady_state_amplitudes(show_material_outlines=False)
+    print(f"Simulating scenario: {scenario_id}")
+    scenario = ndk.make(scenario_id)
+    result = scenario.simulate_steady_state()
+    result.render_steady_state_amplitudes(show_material_outlines=False)
+
 
 # %%
 # Simulating scenario: scenario-0-v0
 # ===================================
-plot_scenario('scenario-0-v0')
+plot_scenario("scenario-0-v0")
 
 # %%
 # Simulating scenario: scenario-1-2d-v0
 # ===================================
-plot_scenario('scenario-1-2d-v0')
+plot_scenario("scenario-1-2d-v0")
 
 # %%
 # Simulating scenario: scenario-2-2d-v0
 # ===================================
-plot_scenario('scenario-2-2d-v0')
+plot_scenario("scenario-2-2d-v0")
diff --git a/docs/examples/plot_store_results.py b/docs/examples/plot_store_results.py
index 51fc0403..bafdec6a 100644
--- a/docs/examples/plot_store_results.py
+++ b/docs/examples/plot_store_results.py
@@ -4,26 +4,29 @@
 ====================================
 
 !!! note
-    NDK and its examples are under constant development, more information and content will be added to this example soon!
+    NDK and its examples are under constant development, more information and
+    content will be added to this example soon!
 
 This example demonstrates how a simulation can be executed and stored in one
-computer and the results can be loaded and rendered later in the same computer or another one.
+computer and the results can be loaded and rendered later in the same computer
+or another one.
 """
 # %%
 # Execute the following code in a computer with ndk installed
 import neurotechdevkit as ndk
 
-scenario = ndk.make('scenario-0-v0')
+scenario = ndk.make("scenario-0-v0")
 result = scenario.simulate_steady_state()
-result.save_to_disk('scenario-0-v0-results.tar.gz')
+result.save_to_disk("scenario-0-v0-results.tar.gz")
 
 
 # %%
 # The output file from the previous step could be copied to
 # another computer with ndk installed or stored for future use. The results can
 # be loaded running the following code:
-import neurotechdevkit as ndk
+import neurotechdevkit as ndk  # noqa: E402
 
-result = ndk.load_result_from_disk('scenario-0-v0-results.tar.gz')
-result.render_steady_state_amplitudes(show_material_outlines=False)
+result2 = ndk.load_result_from_disk("scenario-0-v0-results.tar.gz")
+assert isinstance(result2, ndk.results.SteadyStateResult2D)
+result2.render_steady_state_amplitudes(show_material_outlines=False)
 # %%

From 4aebaf714892d28649ca9b1d01e45e6bf9bcf9dc Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Wed, 5 Jul 2023 14:10:23 +0200
Subject: [PATCH 15/17] Adding alpha symbol to spellcheck whitelist

---
 whitelist.txt | 1 +
 1 file changed, 1 insertion(+)

diff --git a/whitelist.txt b/whitelist.txt
index 899953c8..84f3b8f6 100644
--- a/whitelist.txt
+++ b/whitelist.txt
@@ -1,3 +1,4 @@
+α
 al
 antiparallel
 args

From d2f654d3ea20f402b03cf7410f4bf4c399b219d3 Mon Sep 17 00:00:00 2001
From: Newton Sander <sandernewtonsander@gmail.com>
Date: Wed, 5 Jul 2023 14:29:56 +0200
Subject: [PATCH 16/17] Updating gallery example title

---
 docs/examples/plot_customized_center_frequency.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/examples/plot_customized_center_frequency.py b/docs/examples/plot_customized_center_frequency.py
index 5901ba86..24488782 100644
--- a/docs/examples/plot_customized_center_frequency.py
+++ b/docs/examples/plot_customized_center_frequency.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 """
-Plot customized center frequency
+Custom center frequency
 ====================================
 
 This example demonstrates how to use a customized center frequency

From 22d842658d500a49475b8d18591f1cbe9002b713 Mon Sep 17 00:00:00 2001
From: Diogo de Lucena <90583560+d-lucena@users.noreply.github.com>
Date: Sun, 16 Jul 2023 20:43:24 -0400
Subject: [PATCH 17/17] Update
 docs/examples/plot_customized_center_frequency.py

---
 docs/examples/plot_customized_center_frequency.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/examples/plot_customized_center_frequency.py b/docs/examples/plot_customized_center_frequency.py
index 24488782..a77bad68 100644
--- a/docs/examples/plot_customized_center_frequency.py
+++ b/docs/examples/plot_customized_center_frequency.py
@@ -16,7 +16,7 @@
 # using default material layers
 scenario.material_layers = ["water", "cortical_bone", "brain", "tumor"]
 
-# Customizing material properties with random values
+# Customizing material properties
 scenario.material_properties = {
     "tumor": ndk.materials.Material(
         vp=1850.0, rho=1250.0, alpha=0.8, render_color="#94332F"