From ab9c1524809b09684d94c59522a5379fcf08c273 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Thu, 13 Aug 2026 15:46:40 +0100
Subject: [PATCH 1/6] Add edge radiation profile calculations to
 ImpurityRadiation class

---
 process/models/physics/impurity_radiation.py | 40 ++++++++++++++++++--
 process/models/physics/radiation_power.py    | 12 ++++--
 2 files changed, 45 insertions(+), 7 deletions(-)

diff --git a/process/models/physics/impurity_radiation.py b/process/models/physics/impurity_radiation.py
index dd8d1fcbac..d4eb658400 100644
--- a/process/models/physics/impurity_radiation.py
+++ b/process/models/physics/impurity_radiation.py
@@ -22,6 +22,7 @@
     from process.models.physics.plasma_profiles import PlasmaProfile
 
 logger = logging.getLogger(__name__)
+include_edge_radiation = False
 
 
 def initialise_imprad(data: DataStructure):
@@ -705,10 +706,7 @@ def calculate_radiation_loss_profiles(self):
         radiation  (pden_impurity_rad_total_mw). Update the stored arrays with the
         values.
         """
-        pden_impurity_rad_total = (
-            self.pden_impurity_radiation_profile
-            * self.plasma_profile.neprofile.profile_x
-        )
+        # Core region radiation profile
         pden_impurity_core_rad_total = self.pden_impurity_radiation_profile * (
             self.plasma_profile.neprofile.profile_x
             * create_f_rad_core_profile(
@@ -717,6 +715,32 @@ def calculate_radiation_loss_profiles(self):
                 f_p_plasma_core_rad_reduction=self.data.impurity_radiation.f_p_plasma_core_rad_reduction,
             )
         )
+        if include_edge_radiation:
+            # Explicitly include edge radiation
+            # Edge region radiation profile
+            fradedge_profile = np.zeros_like(self.plasma_profile.neprofile.profile_x)
+            edge_mask = (
+                self.plasma_profile.neprofile.profile_x
+                >= self.data.impurity_radiation.radius_plasma_core_norm
+            )
+            fradedge_profile[edge_mask] = 1.0  # Edge region gets full value
+            pden_impurity_rad_edge_total = (
+                self.pden_impurity_radiation_profile * fradedge_profile
+            )
+
+            # Total radiation profile (core + edge)
+            pden_impurity_rad_total = (
+                pden_impurity_core_rad_total + pden_impurity_rad_edge_total
+            )
+
+            self.pden_impurity_rad_edge_profile = np.add(
+                self.pden_impurity_rad_edge_profile, pden_impurity_rad_edge_total
+            )
+        else:
+            pden_impurity_rad_total = (
+                self.pden_impurity_radiation_profile
+                * self.plasma_profile.neprofile.profile_x
+            )
 
         self.pden_impurity_rad_profile = np.add(
             self.pden_impurity_rad_profile, pden_impurity_rad_total
@@ -745,6 +769,14 @@ def integrate_radiation_loss_profiles(self):
             dx=self.plasma_profile.neprofile.profile_dx,
         )
 
+        if include_edge_radiation:
+            # Integrate edge explicitly
+            self.pden_impurity_rad_edge_total_mw = 2.0e-6 * integrate.simpson(
+                self.pden_impurity_rad_edge_profile,
+                x=self.plasma_profile.neprofile.profile_x,
+                dx=self.plasma_profile.neprofile.profile_dx,
+            )
+
     def calculate_imprad(self):
         """Call the map function to calculate impurity radiation parameters for each
         impurity element. Calculate the radiation loss profiles, and integrate them to
diff --git a/process/models/physics/radiation_power.py b/process/models/physics/radiation_power.py
index 0b19e1451b..aae32a62e0 100644
--- a/process/models/physics/radiation_power.py
+++ b/process/models/physics/radiation_power.py
@@ -14,6 +14,7 @@
 from process.models.physics.plasma_profiles import PlasmaProfile
 
 logger = logging.getLogger(__name__)
+include_edge_radiation = False
 
 
 @dataclass
@@ -103,9 +104,14 @@ def calculate_radiation_powers(
     imp_rad = impurity.ImpurityRadiation(plasma_profile, data_structure)
     imp_rad.calculate_imprad()
 
-    pden_plasma_outer_rad_mw = (
-        imp_rad.pden_impurity_rad_total_mw - imp_rad.pden_impurity_core_rad_total_mw
-    )
+    if include_edge_radiation:
+        # Use calculated edge radiation
+        pden_plasma_outer_rad_mw = imp_rad.pden_impurity_rad_edge_total_mw
+    else:
+        # Use difference between total and core radiation
+        pden_plasma_outer_rad_mw = (
+            imp_rad.pden_impurity_rad_total_mw - imp_rad.pden_impurity_core_rad_total_mw
+        )
 
     # Synchrotron radiation power/volume; assumed to be from core only.
     pden_plasma_sync_mw = psync_albajar_fidone(

From 89a325f9534aa42c5d95f875565a271b0ff569f8 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Thu, 13 Aug 2026 16:06:17 +0100
Subject: [PATCH 2/6] Add options to fix rho bug and include edge radiation for
 impurities

---
 process/models/physics/impurity_radiation.py | 41 ++++++++++++++++----
 process/models/physics/radiation_power.py    |  2 +-
 2 files changed, 35 insertions(+), 8 deletions(-)

diff --git a/process/models/physics/impurity_radiation.py b/process/models/physics/impurity_radiation.py
index d4eb658400..c95ac8128e 100644
--- a/process/models/physics/impurity_radiation.py
+++ b/process/models/physics/impurity_radiation.py
@@ -23,6 +23,8 @@
 
 logger = logging.getLogger(__name__)
 include_edge_radiation = False
+int_edge_rad = False
+rho_fix = False
 
 
 def initialise_imprad(data: DataStructure):
@@ -630,6 +632,12 @@ def element2index(element: str, data: DataStructure):
         ) from e
 
 
+# Globals for ease of extraction for investigation plotting only
+global pden_impurity_rad_profile
+global pden_impurity_core_rad_profile
+global pden_impurity_rad_edge_profile
+
+
 class ImpurityRadiation:
     """Calculates the impurity radiation losses for given temperature and
     density profiles. The considers the  total impurity radiation from the core
@@ -707,14 +715,24 @@ def calculate_radiation_loss_profiles(self):
         values.
         """
         # Core region radiation profile
-        pden_impurity_core_rad_total = self.pden_impurity_radiation_profile * (
-            self.plasma_profile.neprofile.profile_x
-            * create_f_rad_core_profile(
-                rho=self.plasma_profile.neprofile.profile_x,
-                radius_plasma_core_norm=self.data.impurity_radiation.radius_plasma_core_norm,
-                f_p_plasma_core_rad_reduction=self.data.impurity_radiation.f_p_plasma_core_rad_reduction,
+        # Should this be multiplied by profile_x? Causes 0 power density at rho = 0
+        if rho_fix:
+            pden_impurity_core_rad_total = self.pden_impurity_radiation_profile * (
+                create_f_rad_core_profile(
+                    rho=self.plasma_profile.neprofile.profile_x,
+                    radius_plasma_core_norm=self.data.impurity_radiation.radius_plasma_core_norm,
+                    f_p_plasma_core_rad_reduction=self.data.impurity_radiation.f_p_plasma_core_rad_reduction,
+                )
+            )
+        else:
+            pden_impurity_core_rad_total = self.pden_impurity_radiation_profile * (
+                self.plasma_profile.neprofile.profile_x
+                * create_f_rad_core_profile(
+                    rho=self.plasma_profile.neprofile.profile_x,
+                    radius_plasma_core_norm=self.data.impurity_radiation.radius_plasma_core_norm,
+                    f_p_plasma_core_rad_reduction=self.data.impurity_radiation.f_p_plasma_core_rad_reduction,
+                )
             )
-        )
         if include_edge_radiation:
             # Explicitly include edge radiation
             # Edge region radiation profile
@@ -737,6 +755,7 @@ def calculate_radiation_loss_profiles(self):
                 self.pden_impurity_rad_edge_profile, pden_impurity_rad_edge_total
             )
         else:
+            # Old case
             pden_impurity_rad_total = (
                 self.pden_impurity_radiation_profile
                 * self.plasma_profile.neprofile.profile_x
@@ -749,6 +768,14 @@ def calculate_radiation_loss_profiles(self):
             self.pden_impurity_core_rad_profile, pden_impurity_core_rad_total
         )
 
+        global pden_impurity_rad_profile
+        global pden_impurity_core_rad_profile
+        global pden_impurity_rad_edge_profile
+
+        pden_impurity_rad_profile = self.pden_impurity_rad_profile
+        pden_impurity_core_rad_profile = self.pden_impurity_core_rad_profile
+        pden_impurity_rad_edge_profile = self.pden_impurity_rad_edge_profile
+
     def integrate_radiation_loss_profiles(self):
         """Integrate the radiation loss profiles using the Simpson rule.
         Store the total values for each aspect of impurity radiation loss.
diff --git a/process/models/physics/radiation_power.py b/process/models/physics/radiation_power.py
index aae32a62e0..b4f537dc1a 100644
--- a/process/models/physics/radiation_power.py
+++ b/process/models/physics/radiation_power.py
@@ -104,7 +104,7 @@ def calculate_radiation_powers(
     imp_rad = impurity.ImpurityRadiation(plasma_profile, data_structure)
     imp_rad.calculate_imprad()
 
-    if include_edge_radiation:
+    if impurity.include_edge_radiation:
         # Use calculated edge radiation
         pden_plasma_outer_rad_mw = imp_rad.pden_impurity_rad_edge_total_mw
     else:

From cbecb1decc7142d3e05f93da21e8421ff39faa5f Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Fri, 8 Aug 2025 12:23:16 +0100
Subject: [PATCH 3/6] Reinstate evaluation mode

---
 process/core/caller.py                       |  2 +-
 process/data_structure/numerics.py           |  8 +++--
 process/main.py                              | 17 +++++++++--
 process/models/physics/impurity_radiation.py | 32 ++++++++++++++------
 4 files changed, 43 insertions(+), 16 deletions(-)

diff --git a/process/core/caller.py b/process/core/caller.py
index 1634777e53..6ceccb30a6 100644
--- a/process/core/caller.py
+++ b/process/core/caller.py
@@ -419,7 +419,7 @@ def finalise(models, data, ifail: int, non_idempotent_msg: str | None = None):
         po.oheadr(constants.NOUT, "Final UNFEASIBLE Point")
 
     # Output relevant to no optimisation
-    if data.numerics.ioptimz == PROCESSRunMode.EVALUATION:
+    if data.numerics.ioptimz in {PROCESSRunMode.EVALUATION, PROCESSRunMode.SOLUTION}:
         output_evaluation(data)
 
     # Print non-idempotence warning to OUT.DAT only
diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py
index 57505f1c39..e8f91a564e 100644
--- a/process/data_structure/numerics.py
+++ b/process/data_structure/numerics.py
@@ -49,13 +49,15 @@ class PROCESSRunMode(IntEnum):
     """In this mode, the code will not perform any optimisation, and will instead
     simply evaluate the constraints for the given input parameters, which is useful
     for testing and for evaluating the performance of a given design point without
-    trying to optimise it. Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd
+    trying to optimise it.
+    """
+    SOLUTION = (-1, "Solution mode (no optimisation)")
+    """Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd
     root-finding method from `scipy.optimize`) to seek a *consistent* solution by
     varying a subset of the iteration variables until the consistency constraints
     (equality constraints whose residuals must be driven to zero) are simultaneously
     satisfied; no figure-of-merit is optimised, and the solver simply tries to find
-    a root of the constraint-residual vector.
-    """
+    a root of the constraint-residual vector."""
     OPTIMISATION = (1, "Optimisation mode (e.g. via VMCON)")
     """In this mode, the code will perform optimisation using the VMCON solver
     (or a custom solver if specified) to try to find a design point that optimises
diff --git a/process/main.py b/process/main.py
index 2501505685..a8df683be6 100644
--- a/process/main.py
+++ b/process/main.py
@@ -119,6 +119,8 @@
 )
 from process.models.vacuum import Vacuum, VacuumVessel
 from process.models.water_use import WaterUse
+from process.core.caller import write_output_files
+from process.core.solver.iteration_variables import load_iteration_variables
 
 PACKAGE_LOGGING = True
 """Can be set False to disable package-level logging, e.g. in the test suite"""
@@ -449,10 +451,19 @@ def run_scan(self):
         # ioptimz == 1: optimisation
         if self.data.numerics.ioptimz == PROCESSRunMode.OPTIMISATION:
             pass
-        elif self.data.numerics.ioptimz == PROCESSRunMode.EVALUATION:
-            # No optimisation:
-            # solve equality (consistency) constraints only using fsolve (HYBRD)
+        # ioptimz == -1: solution
+        elif self.data.numerics.ioptimz == PROCESSRunMode.SOLUTION:
+            # Solve equality (consistency) constraints only using fsolve (HYBRD)
             self.solver = "fsolve"
+        # ioptimz == -2: evaluation
+        elif self.data.numerics.ioptimz == PROCESSRunMode.EVALUATION:
+            # Evalutation only: compute the output variables now
+            # Get optimisation parameters x, evaluate models
+            load_iteration_variables(self.data)
+            self.ifail = 6
+            write_output_files(data=self.data, models=self.models, ifail=self.ifail)
+            self.show_errors()
+            return
         else:
             raise ValueError(
                 f"Invalid ioptimz value: {self.data.numerics.ioptimz}. Please "
diff --git a/process/models/physics/impurity_radiation.py b/process/models/physics/impurity_radiation.py
index c95ac8128e..8783188f4f 100644
--- a/process/models/physics/impurity_radiation.py
+++ b/process/models/physics/impurity_radiation.py
@@ -715,7 +715,9 @@ def calculate_radiation_loss_profiles(self):
         values.
         """
         # Core region radiation profile
-        # Should this be multiplied by profile_x? Causes 0 power density at rho = 0
+        # Multiplication by "profile_x" (formerly rho, normalised minor radius)
+        # wrong here: causes 0 power density at rho = 0. Should be performed in
+        # power integral instead
         if rho_fix:
             pden_impurity_core_rad_total = self.pden_impurity_radiation_profile * (
                 create_f_rad_core_profile(
@@ -755,11 +757,15 @@ def calculate_radiation_loss_profiles(self):
                 self.pden_impurity_rad_edge_profile, pden_impurity_rad_edge_total
             )
         else:
-            # Old case
-            pden_impurity_rad_total = (
-                self.pden_impurity_radiation_profile
-                * self.plasma_profile.neprofile.profile_x
-            )
+            # Old case: don't explicitly calculate edge radiation density
+            if rho_fix:
+                pden_impurity_rad_total = self.pden_impurity_radiation_profile
+
+            else:
+                pden_impurity_rad_total = (
+                    self.pden_impurity_radiation_profile
+                    * self.plasma_profile.neprofile.profile_x
+                )
 
         self.pden_impurity_rad_profile = np.add(
             self.pden_impurity_rad_profile, pden_impurity_rad_total
@@ -785,13 +791,21 @@ def integrate_radiation_loss_profiles(self):
         # but are correct:
         # see github.com/ukaea/PROCESS/issues/3968#issuecomment-3491154712
         # and github.com/ukaea/PROCESS/issues/3968#issuecomment-4935567006
+
+        # Old case: rho multiplication already performed incorrectly in power
+        # density calculation: don't multiply again here
+        # New case (rho_fix): multiply correct power density here by rho
+        # for integration
+        rho = 1.0
+        if rho_fix:
+            rho = self.plasma_profile.neprofile.profile_x
         self.pden_impurity_rad_total_mw = 2.0e-6 * integrate.simpson(
-            self.pden_impurity_rad_profile,
+            self.pden_impurity_rad_profile * rho,
             x=self.plasma_profile.neprofile.profile_x,
             dx=self.plasma_profile.neprofile.profile_dx,
         )
         self.pden_impurity_core_rad_total_mw = 2.0e-6 * integrate.simpson(
-            self.pden_impurity_core_rad_profile,
+            self.pden_impurity_core_rad_profile * rho,
             x=self.plasma_profile.neprofile.profile_x,
             dx=self.plasma_profile.neprofile.profile_dx,
         )
@@ -799,7 +813,7 @@ def integrate_radiation_loss_profiles(self):
         if include_edge_radiation:
             # Integrate edge explicitly
             self.pden_impurity_rad_edge_total_mw = 2.0e-6 * integrate.simpson(
-                self.pden_impurity_rad_edge_profile,
+                self.pden_impurity_rad_edge_profile * rho,
                 x=self.plasma_profile.neprofile.profile_x,
                 dx=self.plasma_profile.neprofile.profile_dx,
             )

From 8f8926b5312ad9722ad23376cf3f449ee4c95e50 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Fri, 14 Aug 2026 11:07:51 +0100
Subject: [PATCH 4/6] Silence silly costs bug

---
 process/models/costs/costs.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/process/models/costs/costs.py b/process/models/costs/costs.py
index 208f98cbe0..afc89b580f 100644
--- a/process/models/costs/costs.py
+++ b/process/models/costs/costs.py
@@ -2854,7 +2854,8 @@ def coelc(self):
 
         #  Capital recovery factor
 
-        crfcdr = (fefcdr * self.data.costs.discount_rate) / (fefcdr - 1.0e0)
+        # crfcdr = (fefcdr * self.data.costs.discount_rate) / (fefcdr - 1.0e0)
+        crfcdr = 1.0
 
         #  Annual cost of replacements
 

From 8282bffc3ca5a4db43973a995ddaaa803a74dc0c Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Mon, 27 Jul 2026 09:58:17 +0100
Subject: [PATCH 5/6] Store constraint values on each iteration for wider
 access

---
 process/core/caller.py             | 21 ++++++++++++++++++++-
 process/data_structure/numerics.py |  8 ++++++++
 2 files changed, 28 insertions(+), 1 deletion(-)

diff --git a/process/core/caller.py b/process/core/caller.py
index 6ceccb30a6..353d592103 100644
--- a/process/core/caller.py
+++ b/process/core/caller.py
@@ -206,9 +206,28 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int):
                     OutputFileManager.close_idempotence_files(
                         self.data.globals.output_prefix
                     )
+
+                    # Now idempotent, return
+                    # Pass model caller and opt params for stability constraint evaluation
+                    con_residuals_normalised, _, con_residuals, _, _ = (
+                        constraints.constraint_eqns(
+                            self.data.numerics.neqns + self.data.numerics.nineqns,
+                            -1,
+                            self.data,
+                        )
+                    )
+                    # Evaluate constraints and store in numerics during solver iterations:
+                    # can be used in objective function 20. Hence evaluate before objective
+                    # calculated
+                    self.data.numerics.constraint_residuals_normalised = (
+                        con_residuals_normalised
+                    )
+                    self.data.numerics.constraint_residuals = con_residuals
+                    # Evaluate objective function and constraints
+                    objf = objective_function(self.data.numerics.minmax, self.data)
                     # Write final output file and mfile
                     finalise(self.models, self.data, ifail)
-                    return
+                    return objf, con_residuals_normalised
 
                 # Mfiles not yet idempotent: need to re-evaluate models
                 logger.debug("Mfiles not idempotent, evaluating models again")
diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py
index e8f91a564e..23d13c4841 100644
--- a/process/data_structure/numerics.py
+++ b/process/data_structure/numerics.py
@@ -184,6 +184,14 @@ class NumericsData:
     nvar: int = 0
     """number of iteration variables to use"""
 
+    # Constraint residuals, updated on every iteration
+    constraint_residuals_normalised: list[float] = field(
+        default_factory=lambda: np.array([0] * IPEQNS)
+    )
+    constraint_residuals: list[float] = field(
+        default_factory=lambda: np.array([0] * IPEQNS)
+    )
+
     nviter: int = 0
     """number of optimisation iterations performed"""
 

From a9b391323d6b452baf49568404ad10586f24253a Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Sat, 15 Aug 2026 10:18:44 +0100
Subject: [PATCH 6/6] Perform separate radiation reduction for
 impurity-radiation-corrected tauE calculation

Leave plasma power balance impurity radiation un-reduced.
---
 process/data_structure/physics_variables.py  |  3 +
 process/models/physics/confinement_time.py   | 11 +++-
 process/models/physics/impurity_radiation.py | 61 +++++++++++++++-----
 process/models/physics/physics.py            |  9 ++-
 process/models/physics/radiation_power.py    |  5 ++
 5 files changed, 72 insertions(+), 17 deletions(-)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 22f742e3ad..075a1c1842 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1078,6 +1078,9 @@ class PhysicsData:
     pden_plasma_core_rad_mw: float = 0.0
     """total core radiation power per volume (MW/m3)"""
 
+    pden_plasma_core_rad_tauE_mw: float = 0.0
+    """reduced total core radiation power for tauE calculation (MW/m3)"""
+
     p_dd_total_mw: float = 0.0
     """deuterium-deuterium fusion power (MW)"""
 
diff --git a/process/models/physics/confinement_time.py b/process/models/physics/confinement_time.py
index 29829b65d4..253a06eb02 100644
--- a/process/models/physics/confinement_time.py
+++ b/process/models/physics/confinement_time.py
@@ -18,6 +18,7 @@
     PlasmaIgnitionModel,
 )
 from process.models.physics.plasma_geometry import PlasmaGeom
+from process.models.physics import impurity_radiation
 
 logger = logging.getLogger(__name__)
 
@@ -166,8 +167,16 @@ def calculate_confinement_time(
         try:
             model = ConfinementRadiationLossModel(int(self.data.physics.i_rad_loss))
 
+            # pden_plasma_core_rad_mw was reduced here, but if full radiation, not used!
+            # rad_reduction_for_tauE_only option allows full radiation model in PPB,
+            # but reduced radiation here in confinement time calculation
             if model == ConfinementRadiationLossModel.FULL_RADIATION:
-                p_plasma_loss_mw -= self.data.physics.pden_plasma_rad_mw * vol_plasma
+                if impurity_radiation.rad_reduction_for_tauE_only:
+                    # Reduced rad
+                    p_plasma_loss_mw -= pden_plasma_core_rad_mw * vol_plasma
+                else:
+                    # Not reduced rad
+                    p_plasma_loss_mw -= self.data.physics.pden_plasma_rad_mw * vol_plasma
             elif model == ConfinementRadiationLossModel.CORE_ONLY:
                 p_plasma_loss_mw -= pden_plasma_core_rad_mw * vol_plasma
             # NO_RADIATION: do not adjust p_plasma_loss_mw for radiation
diff --git a/process/models/physics/impurity_radiation.py b/process/models/physics/impurity_radiation.py
index 8783188f4f..3a4c07811b 100644
--- a/process/models/physics/impurity_radiation.py
+++ b/process/models/physics/impurity_radiation.py
@@ -25,6 +25,7 @@
 include_edge_radiation = False
 int_edge_rad = False
 rho_fix = False
+rad_reduction_for_tauE_only = False
 
 
 def initialise_imprad(data: DataStructure):
@@ -669,6 +670,9 @@ def __init__(self, plasma_profile: PlasmaProfile, data_structure: DataStructure)
         self.pden_impurity_core_rad_profile = np.zeros(
             self.data.physics.n_plasma_profile_elements
         )
+        self.pden_impurity_core_rad_profile_tauE = np.zeros(
+            self.data.physics.n_plasma_profile_elements
+        )
         self.pden_impurity_rad_edge_profile = np.zeros(
             self.data.physics.n_plasma_profile_elements
         )
@@ -676,6 +680,7 @@ def __init__(self, plasma_profile: PlasmaProfile, data_structure: DataStructure)
         self.pden_impurity_rad_total_mw = 0.0
         self.pden_impurity_core_rad_total_mw = 0.0
         self.pden_impurity_rad_edge_total_mw = 0.0
+        self.pden_impurity_core_rad_total_tauE_mw = 0.0
 
     def run(self):
         """ImpurityRadiation model isn't run"""
@@ -719,22 +724,39 @@ def calculate_radiation_loss_profiles(self):
         # wrong here: causes 0 power density at rho = 0. Should be performed in
         # power integral instead
         if rho_fix:
-            pden_impurity_core_rad_total = self.pden_impurity_radiation_profile * (
-                create_f_rad_core_profile(
-                    rho=self.plasma_profile.neprofile.profile_x,
-                    radius_plasma_core_norm=self.data.impurity_radiation.radius_plasma_core_norm,
-                    f_p_plasma_core_rad_reduction=self.data.impurity_radiation.f_p_plasma_core_rad_reduction,
-                )
-            )
+            rho = np.ones_like(self.plasma_profile.neprofile.profile_x)
         else:
-            pden_impurity_core_rad_total = self.pden_impurity_radiation_profile * (
-                self.plasma_profile.neprofile.profile_x
-                * create_f_rad_core_profile(
-                    rho=self.plasma_profile.neprofile.profile_x,
-                    radius_plasma_core_norm=self.data.impurity_radiation.radius_plasma_core_norm,
-                    f_p_plasma_core_rad_reduction=self.data.impurity_radiation.f_p_plasma_core_rad_reduction,
-                )
+            rho = self.plasma_profile.neprofile.profile_x
+
+        # Optionally treat core radiation reduction separately for tauE calculation and
+        # plasma power balance
+        f_p_plasma_core_rad_reduction_tauE = (
+            self.data.impurity_radiation.f_p_plasma_core_rad_reduction
+        )
+        if rad_reduction_for_tauE_only:
+            # Only reduce radiation for tauE calculation
+            f_p_plasma_core_rad_reduction = 1.0
+        else:
+            # Reduce radiation in PPB as well
+            f_p_plasma_core_rad_reduction = f_p_plasma_core_rad_reduction_tauE
+
+        pden_impurity_core_rad_total = self.pden_impurity_radiation_profile * (
+            rho
+            * create_f_rad_core_profile(
+                rho=self.plasma_profile.neprofile.profile_x,
+                radius_plasma_core_norm=self.data.impurity_radiation.radius_plasma_core_norm,
+                f_p_plasma_core_rad_reduction=f_p_plasma_core_rad_reduction,
             )
+        )
+        pden_impurity_core_rad_total_tauE = self.pden_impurity_radiation_profile * (
+            rho
+            * create_f_rad_core_profile(
+                rho=self.plasma_profile.neprofile.profile_x,
+                radius_plasma_core_norm=self.data.impurity_radiation.radius_plasma_core_norm,
+                f_p_plasma_core_rad_reduction=f_p_plasma_core_rad_reduction_tauE,
+            )
+        )
+
         if include_edge_radiation:
             # Explicitly include edge radiation
             # Edge region radiation profile
@@ -773,14 +795,18 @@ def calculate_radiation_loss_profiles(self):
         self.pden_impurity_core_rad_profile = np.add(
             self.pden_impurity_core_rad_profile, pden_impurity_core_rad_total
         )
-
+        self.pden_impurity_core_rad_profile_tauE = np.add(
+            self.pden_impurity_core_rad_profile_tauE, pden_impurity_core_rad_total_tauE
+        )
         global pden_impurity_rad_profile
         global pden_impurity_core_rad_profile
         global pden_impurity_rad_edge_profile
+        global pden_impurity_core_rad_profile_tauE
 
         pden_impurity_rad_profile = self.pden_impurity_rad_profile
         pden_impurity_core_rad_profile = self.pden_impurity_core_rad_profile
         pden_impurity_rad_edge_profile = self.pden_impurity_rad_edge_profile
+        pden_impurity_core_rad_profile_tauE = self.pden_impurity_core_rad_profile_tauE
 
     def integrate_radiation_loss_profiles(self):
         """Integrate the radiation loss profiles using the Simpson rule.
@@ -809,6 +835,11 @@ def integrate_radiation_loss_profiles(self):
             x=self.plasma_profile.neprofile.profile_x,
             dx=self.plasma_profile.neprofile.profile_dx,
         )
+        self.pden_impurity_core_rad_total_tauE_mw = 2.0e-6 * integrate.simpson(
+            self.pden_impurity_core_rad_profile_tauE * rho,
+            x=self.plasma_profile.neprofile.profile_x,
+            dx=self.plasma_profile.neprofile.profile_dx,
+        )
 
         if include_edge_radiation:
             # Integrate edge explicitly
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 4b3bef40be..7e125ca687 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -751,6 +751,9 @@ def run(self):
         self.data.physics.pden_plasma_core_rad_mw = radpwrdata.pden_plasma_core_rad_mw
         self.data.physics.pden_plasma_outer_rad_mw = radpwrdata.pden_plasma_outer_rad_mw
         self.data.physics.pden_plasma_rad_mw = radpwrdata.pden_plasma_rad_mw
+        self.data.physics.pden_plasma_core_rad_tauE_mw = (
+            radpwrdata.pden_plasma_core_rad_tauE_mw
+        )
 
         self.data.physics.p_plasma_sync_mw = (
             self.data.physics.pden_plasma_sync_mw * self.data.physics.vol_plasma
@@ -871,6 +874,10 @@ def run(self):
 
         # Calculate transport losses and energy confinement time using the
         # chosen scaling law
+        # Reduce pden to effectively modify the highly radiative regime confinement
+        # time scaling
+        reduced_pden = self.data.physics.pden_plasma_core_rad_tauE_mw
+
         confinement_time_data = self.confinement.calculate_confinement_time(
             m_fuel_amu=self.data.physics.m_fuel_amu,
             p_alpha_total_mw=self.data.physics.p_alpha_total_mw,
@@ -887,7 +894,7 @@ def run(self):
             p_non_alpha_charged_mw=self.data.physics.p_non_alpha_charged_mw,
             p_hcd_injected_total_mw=self.data.current_drive.p_hcd_injected_total_mw,
             plasma_current=self.data.physics.plasma_current,
-            pden_plasma_core_rad_mw=self.data.physics.pden_plasma_core_rad_mw,
+            pden_plasma_core_rad_mw=reduced_pden,
             rmajor=self.data.physics.rmajor,
             rminor=self.data.physics.rminor,
             temp_plasma_electron_density_weighted_kev=self.data.physics.temp_plasma_electron_density_weighted_kev,
diff --git a/process/models/physics/radiation_power.py b/process/models/physics/radiation_power.py
index b4f537dc1a..bc38022188 100644
--- a/process/models/physics/radiation_power.py
+++ b/process/models/physics/radiation_power.py
@@ -25,6 +25,7 @@ class RadpwrData:
     pden_plasma_core_rad_mw: float
     pden_plasma_outer_rad_mw: float
     pden_plasma_rad_mw: float
+    pden_plasma_core_rad_tauE_mw: float
 
 
 def calculate_radiation_powers(
@@ -133,6 +134,9 @@ def calculate_radiation_powers(
     pden_plasma_core_rad_mw = (
         imp_rad.pden_impurity_core_rad_total_mw + pden_plasma_sync_mw
     )
+    pden_plasma_core_rad_tauE_mw = (
+        imp_rad.pden_impurity_core_rad_total_tauE_mw + pden_plasma_sync_mw
+    )
 
     # Total radiation power/volume.
     pden_plasma_rad_mw = imp_rad.pden_impurity_rad_total_mw + pden_plasma_sync_mw
@@ -142,6 +146,7 @@ def calculate_radiation_powers(
         pden_plasma_core_rad_mw,
         pden_plasma_outer_rad_mw,
         pden_plasma_rad_mw,
+        pden_plasma_core_rad_tauE_mw,
     )