make hadrons continous runable, add progressbars, improve create_trac…

…k_distribution

electrons/common/generic_electron_solver.py

@@ -6,6 +6,7 @@
 
 import numpy as np
 from numpy.typing import NDArray
+from tqdm import tqdm
 
 
 def create_sc_gradients(
@@ -174,7 +175,7 @@ def calculate(self):
         Start the simulation by evolving the distribution one step at a time
         """
 
-        for time_step in range(self.computation_time_steps):
+        for time_step in tqdm(range(self.computation_time_steps)):
 
             """
             Refill the array with the electron density each time step

hadrons/common/continous_beam.py

@@ -1,10 +1,11 @@
 from dataclasses import dataclass
 from math import exp, sqrt
-from random import Random
 
 import numpy as np
+from numpy.random import Generator, default_rng
 
 from hadrons.common.generic_hadron_solver import GenericHadronSolver
+from hadrons.utils.common import doserate_to_fluence
 from hadrons.utils.track_distribution import create_track_distribution
 
 
@@ -13,9 +14,9 @@ def get_continous_beam_pde_solver(base_solver_class: GenericHadronSolver):
     @dataclass
     class ContinousHadronSolver(base_solver_class):
         # TODO: dataclasses do not allow for non-default properties if the inherided class has default properties - not sure how to fix this yet
-        fluence_rate: float = None  # [cm-^2/s]
+        doserate: float = None  # [cm-^2/s]
         # TODO: add a seed param
-        _random_generator: Random = Random(2137)
+        default_random_generator: Generator = default_rng(2137)
 
         @property
         def buffer_radius(self):
@@ -30,7 +31,7 @@ def no_xy(self) -> int:
 
         @property
         def random_generator(self):
-            return self._random_generator
+            return self.default_random_generator
 
         @property
         def xy_middle_idx(self):
@@ -42,23 +43,32 @@ def inner_radius(self):
 
             return outer_radius - self.buffer_radius
 
+        @property
+        def fluence_rate(self):
+            return doserate_to_fluence(
+                self.doserate, self.energy, particle=self.particle
+            )
+
         def __post_init__(self):
             super().__post_init__()
             self.simulation_time = self.computation_time_steps * self.dt
             self.number_of_tracks = int(
                 self.fluence_rate * self.simulation_time * self.track_area
             )
+
             self.number_of_tracks = max(1, self.number_of_tracks)
 
             self.track_distribution = create_track_distribution(
                 self.computation_time_steps,
                 self.dt,
                 self.number_of_tracks,
                 self.separation_time_steps,
-                self._random_generator,
+                self.random_generator,
             )
 
         def get_number_of_tracks(self, time_step: int) -> bool:
+            if time_step == 0:
+                print(self.track_distribution)
             return self.track_distribution[time_step]
 
         def get_random_xy_coordinate(self):

hadrons/common/generic_hadron_solver.py

@@ -5,6 +5,7 @@
 
 import numpy as np
 from numpy.typing import NDArray
+from tqdm import tqdm
 
 from hadrons.utils.common import calculate_track_radius, get_LET_per_um
 
@@ -109,12 +110,15 @@ def electric_field(self) -> float:
         return self.voltage / self.electrode_gap
 
     @property
-    def computation_time_steps(self) -> int:
-        separation_time_steps = int(
+    def separation_time_steps(self) -> int:
+        return int(
             self.electrode_gap / (2.0 * ion_mobility * self.electric_field * self.dt)
         )
 
-        return separation_time_steps * 2
+    @property
+    def computation_time_steps(self) -> int:
+
+        return self.separation_time_steps * 2
 
     @property
     def RDD_function(self):
@@ -203,7 +207,9 @@ def calculate(self):
 
         sc_pos, sc_neg, sc_center = create_sc_gradients(self.s, self.c)
 
-        for time_step in range(self.computation_time_steps):
+        for time_step in tqdm(
+            range(self.computation_time_steps), desc="Calculating..."
+        ):
             # calculate the new densities and store them in temporary arrays
 
             for _ in range(self.get_number_of_tracks(time_step)):
@@ -269,8 +275,14 @@ def calculate(self):
             positive_array[:] = positive_array_temp[:]
             negative_array[:] = negative_array_temp[:]
 
-        f_steps_list[time_step] = (
-            no_initialised_charge_carriers - no_recombined_charge_carriers
-        ) / no_initialised_charge_carriers
+            print(no_initialised_charge_carriers, no_recombined_charge_carriers)
+
+            # calculate the fraction of charge carriers that have not recombined, if no charge carriers have been initialised, set the fraction to 1
+            if no_initialised_charge_carriers != 0:
+                f_steps_list[time_step] = (
+                    no_initialised_charge_carriers - no_recombined_charge_carriers
+                ) / no_initialised_charge_carriers
+            else:
+                f_steps_list[time_step] = 1
 
         return f_steps_list

hadrons/python/run_simulation.py

@@ -14,6 +14,7 @@ def run_simulation(
     electrode_gap=0.1,
     energy=200,
     electron_density_per_cm3=1e9,
+    doserate=100,
     verbose=False,
 ):
     if solver_name not in SOLVER_MAP.keys():
@@ -28,13 +29,16 @@ def run_simulation(
         print(f"Voltage: {voltage} [V]")
         print(f"Electrode gap: {electrode_gap} [cm]")
         print(f"Electron density per cm3: {electron_density_per_cm3}")
+        if solver_name == "continous":
+            print(f"doserate: {doserate} [Gy/s]")
 
     solver = Solver(
         # LET=LET,
         voltage=voltage,  # [V/cm] magnitude of the electric field
         IC_angle=IC_angle,
         electrode_gap=electrode_gap,
         energy=energy,
+        doserate=doserate,
     )
 
     return solver.calculate()