Dump sources

pyproject.toml

@@ -26,6 +26,7 @@ dependencies = [
     "pandas",
     "jsonargparse",
     "scipy",
+    "actigamma",
 ]
 [project.optional-dependencies]
 tests = [

src/f4epurity/decay_chain_calc.py

@@ -1,6 +1,6 @@
+import math
 from copy import deepcopy
 
-import math
 import numpy as np
 
 from f4epurity.utilities import add_user_irrad_scenario, convert_names
@@ -42,6 +42,10 @@
         + [3920, 400] * 3,
         "fluxes": [0.00536, 0.0412, 0, 0.083] + [0, 1] * 17 + [0, 1.4] * 3,
     },
+    "Y1": {
+        "times": [365 * DAY_TO_SEC],
+        "fluxes": [1],
+    },
 }
 
 
@@ -105,7 +109,6 @@ def irradiate(
         decay_data_dic[nuclide]["lambda"][1:],
         decay_data_dic[nuclide]["Decay_daughter_names"][1:],
     ):
-
         # Tool cannot handle fission or very long lived >1e19 s half-lives, therefore treat these as stable.
         if decay_type == "f" or ("BR" in decay_data_dic[nuclide] and lambd < 1e-20):
             continue
@@ -168,14 +171,12 @@ def get_nuclides(scenario: dict, parent: str, decay_data_dic: dict, atoms: int)
 
     # Loop over each time step in the irradiation and decay
     for time, flux in zip(scenario["times"], scenario["fluxes"]):
-
         # initialise a dictionary of nuclides which will occur
         decay_data_dic[parent]["lambda"] = [flux * x for x in initial_lambda]
         sum_nuclides = {}
 
         # Loop over each of the nuclides that is present at the start of the pulse
         for nuclide in init_nuclides:
-
             # Set the first lambda in the decay chain to the lambda of this nuclide
             lambda_temp[0] = decay_data_dic[nuclide]["lambda"][0]
 
@@ -268,7 +269,7 @@ def create_dictionary(decay_data: dict, parent: str, daughters: list) -> dict:
 
 def calculate_total_activity(
     nuclide_dict: dict, irrad_scenario: dict, decay_time: int, decay_data: dict
-) -> dict:
+) -> dict[str, list[np.ndarray]]:
     """Function to calculate the total activity of a given nuclide dictionary
 
     Parameters
@@ -284,7 +285,7 @@ def calculate_total_activity(
 
     Returns
     -------
-    dict
+    dict[str, list[np.ndarray]]
         dictionary containing the total activity for each nuclide in the nuclide dictionary
     """
 

src/f4epurity/dose.py

@@ -1,10 +1,10 @@
+import logging
 import os
 import sys
+from importlib.resources import as_file, files
 from math import pi
 from pathlib import Path
 
-from importlib.resources import files, as_file
-import logging
 import matplotlib.pyplot as plt
 import numpy as np
 import pyevtk
@@ -102,7 +102,6 @@ def dose_from_line_source(dose, x1, y1, z1, x2, y2, z2, x, y, z):
 def is_within_bounds(
     x1, y1, z1, flux_grid_file, stl_files_path=None, x2=None, y2=None, z2=None
 ):
-
     # If a path to the STL files is provided, use it
     if stl_files_path is not None:
         folder_path = Path(stl_files_path)
@@ -178,7 +177,6 @@ def write_vtk_file(
     z2=None,
     output_all_vtr=False,
 ):
-
     # Get the bounds in which to make the plot
     plot_bounds = is_within_bounds(
         x1, y1, z1, flux_grid_file, stl_files_path, x2, y2, z2
@@ -227,28 +225,40 @@ def write_vtk_file(
                     )
                     dose_array[i, j, k] = dose_point
 
-    # Create a 3D numpy array filled with the dose values based on 1/r^2
+    # Create a 3D numpy array filled with the dose values based on 1/r^2 ( original bottom-left point source model)
     else:
         for i in range(len(x) - 1):
             for j in range(len(y) - 1):
                 for k in range(len(z) - 1):
                     distance = np.sqrt(
                         (x[i] - x1) ** 2 + (y[j] - y1) ** 2 + (z[k] - z1) ** 2
                     )
+                    ## get midpoints
+                    # x_mid = (x[1:] + x[:-1]) / 2
+                    # y_mid = (y[1:] + y[:-1]) / 2
+                    # z_mid = (z[1:] + z[:-1]) / 2
+                    # for i in range(len(x_mid)):
+                    #     for j in range(len(y_mid)):
+                    #         for k in range(len(z_mid)):
+                    #             distance = np.sqrt(
+                    #                 (x_mid[i] - x1) ** 2 + (y_mid[j] - y1) ** 2 + (z_mid[k] - z1) ** 2
+                    #             )
                     # Avoid division by zero
                     if distance == 0:
                         dose_array[i, j, k] = dose[0]
+                        print("Dose at source point is:", dose[0])
                     else:
                         dose_array[i, j, k] = dose[0] / (4 * pi * distance**2)
-
     # Get the indices of the max value
     indices = np.unravel_index(np.argmax(dose_array, axis=None), dose_array.shape)
 
     # Get the coordinates of the max value
     x_max = x[indices[0]]
     y_max = y[indices[1]]
     z_max = z[indices[2]]
-
+    # print(
+    #     f"Max dose value: {dose_array[indices]} at coordinates: ({x_max}, {y_max}, {z_max})"
+    # )
     # Create the output directory if it doesn't exist
     os.makedirs("output", exist_ok=True)
 
@@ -265,7 +275,6 @@ def write_vtk_file(
 
 
 def plot_slice(dose_array, x, y, z, slice_axis, slice_location, plot_bounds):
-
     # Map axis names to indices
     axis_dict = {"x": 0, "y": 1, "z": 2}
 
@@ -352,8 +361,9 @@ def format_func(value, tick_number):
         # Loop over each solid in the stl file and plot the slice
         for solid in stl.solids:
             slice = solid.slice(plane=slice_axis, intercept=slice_location)
-            pairs = getattr(slice, axis_pairs[slice_axis][0] + "_pairs"), getattr(
-                slice, axis_pairs[slice_axis][1] + "_pairs"
+            pairs = (
+                getattr(slice, axis_pairs[slice_axis][0] + "_pairs"),
+                getattr(slice, axis_pairs[slice_axis][1] + "_pairs"),
             )
             for pair in zip(*pairs):
                 ax.plot(*pair, color="black", linewidth=1)

src/f4epurity/main.py

@@ -1,30 +1,32 @@
-import logging
-from jsonargparse import Namespace
 import csv
 import datetime
 import json
-import numpy as np
+import logging
 import os
-from importlib.resources import files, as_file
+from importlib.resources import as_file, files
+
+import numpy as np
 import pandas as pd
+from jsonargparse import Namespace
 
-from f4epurity.decay_chain_calc import calculate_total_activity
 from f4epurity.collapse import collapse_flux, extract_xs
-from f4epurity.dose import convert_to_dose, write_vtk_file, plot_slice
+from f4epurity.decay_chain_calc import calculate_total_activity
+from f4epurity.dose import convert_to_dose, plot_slice, write_vtk_file
 from f4epurity.maintenance import (
     dose_within_workstation,
     get_dose_at_workstation,
     read_maintenance_locations,
 )
+from f4epurity.parsing import parse_arguments, parse_isotopes_activities_file
+from f4epurity.psource import GlobalPointSource, PointSource
 from f4epurity.reaction_rate import calculate_reaction_rate
 from f4epurity.utilities import (
     calculate_number_of_atoms,
     get_isotopes,
-    sum_vtr_files,
-    normalise_nuclide_name,
     get_reactions_from_file,
+    normalise_nuclide_name,
+    sum_vtr_files,
 )
-from f4epurity.parsing import parse_arguments, parse_isotopes_activities_file
 
 F4Epurity_TITLE = """
   _____ _  _   _____                  _ _         
@@ -76,7 +78,14 @@ def calculate_dose_for_source(
     x2: np.ndarray = None,
     y2: np.ndarray = None,
     z2: np.ndarray = None,
-) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, list[float]]:
+) -> tuple[
+    np.ndarray,
+    np.ndarray,
+    np.ndarray,
+    np.ndarray,
+    list[float],
+    dict[str, list[np.ndarray]],
+]:
     """Calculate the dose for a given source
 
     Parameters
@@ -141,7 +150,6 @@ def calculate_dose_for_source(
             sigma_eff, flux_spectrum = collapse_flux(
                 xs_values, args.input_flux, x1, y1, z1, x2, y2, z2
             )
-
             # Calculate the reaction rate based on the flux and effective cross section
             reaction_rate = calculate_reaction_rate(
                 args.delta_impurity, sigma_eff, flux_spectrum
@@ -161,7 +169,6 @@ def calculate_dose_for_source(
     logging.info("Calculating the Dose...")
     # Determine the Dose for each nuclide
     for nuclide, nuclide_activity in activities.items():
-
         # Convert to format in dose conversion spreadsheet
         nuclide = normalise_nuclide_name(nuclide)
 
@@ -210,7 +217,7 @@ def calculate_dose_for_source(
         else:
             plt.savefig(f"{run_dir}/dose_{x1}_{y1}_{z1}.png")
 
-    return dose_array, x, y, z, total_dose
+    return dose_array, x, y, z, total_dose, activities
 
 
 def calculate_dose_at_workstations(
@@ -301,6 +308,7 @@ def process_sources(args: Namespace) -> None:
         dose_factors_df = pd.read_excel(fp)
 
     dose_arrays = []
+    sources = []
     # Check if a second point was provided - line source
     if args.x2 is not None and args.y2 is not None and args.z2 is not None:
         logging.info("Line source(s) selected")
@@ -309,7 +317,7 @@ def process_sources(args: Namespace) -> None:
         for x1, y1, z1, x2, y2, z2 in zip(
             args.x1, args.y1, args.z1, args.x2, args.y2, args.z2
         ):
-            dose_array, x, y, z, dose = calculate_dose_for_source(
+            dose_array, x, y, z, dose, activities = calculate_dose_for_source(
                 args,
                 x1,
                 y1,
@@ -339,8 +347,8 @@ def process_sources(args: Namespace) -> None:
         logging.info("Point source(s) selected")
 
         # Handle multiple coordinates being provided
-        for x1, y1, z1 in zip(args.x1, args.y1, args.z1):
-            dose_array, x, y, z, dose = calculate_dose_for_source(
+        for i, (x1, y1, z1) in enumerate(zip(args.x1, args.y1, args.z1)):
+            dose_array, x, y, z, dose, activities = calculate_dose_for_source(
                 args,
                 x1,
                 y1,
@@ -352,8 +360,18 @@ def process_sources(args: Namespace) -> None:
                 dose_factors_df,
             )
             dose_arrays.append(dose_array)
+            if args.dump_source:
+                if args.m:
+                    mass = args.m[i]
+                else:
+                    mass = 1
+                sources.append(PointSource(activities, [x1, y1, z1], mass=mass))
             calculate_dose_at_workstations(args, dose, x1, y1, z1, run_dir)
 
+    if args.dump_source:
+        global_source = GlobalPointSource(sources)
+        global_source.to_sdef(f"{run_dir}/source.sdef")  # TODO
+
     # If more than one dose array is present, sum the dose arrays (multiple sources)
     if len(dose_arrays) > 1:
         sum_vtr_files(dose_arrays, x, y, z, run_dir, masses=args.m)

src/f4epurity/parsing.py

@@ -1,9 +1,9 @@
+import os
 from typing import List, Optional
 
-import os
 import numpy as np
-from jsonargparse import ArgumentParser, Namespace, ActionConfigFile
 import pandas as pd
+from jsonargparse import ActionConfigFile, ArgumentParser, Namespace
 
 
 def parse_arguments(args_list: Optional[List[str]] = None) -> Namespace:
@@ -138,7 +138,12 @@ def parse_arguments(args_list: Optional[List[str]] = None) -> Namespace:
         type=str,
         help="Path to STL files of ITER rooms that can be used for plotting",
     )
-
+    # Add the optional "mcnp" argument
+    parser.add_argument(
+        "--dump_source",
+        action="store_true",
+        help="Optional 'dump_source' argument to generate source.sdef.file for MCNP",
+    )
     # Parse the arguments
     args = parser.parse_args(args_list)