From c3ff25d78dff2e4a59e5a3e19de59300c497dc0d Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 16:16:29 +0000 Subject: [PATCH 01/20] Add deuteron-triton Coulomb logarithm profile to physics variables --- process/data_structure/physics_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 7f6adaa12..2dce159e2 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1329,6 +1329,9 @@ plasma_coulomb_log_electron_triton_profile: list[float] = None """Profile of electron-triton Coulomb logarithm in plasma""" +plasma_coulomb_log_deuteron_triton_profile: list[float] = None +"""Profile of deuteron-triton Coulomb logarithm in plasma""" + freq_plasma_electron_profile: list[float] = None """Electron plasma frequency profile (Hz)""" @@ -1666,6 +1669,7 @@ def init_physics_variables(): plasma_coulomb_log_electron_electron_profile, \ plasma_coulomb_log_electron_deuteron_profile, \ plasma_coulomb_log_electron_triton_profile, \ + plasma_coulomb_log_deuteron_triton_profile, \ freq_plasma_electron_profile, \ freq_plasma_larmor_toroidal_electron_profile, \ freq_plasma_larmor_toroidal_deuteron_profile, \ @@ -1935,6 +1939,7 @@ def init_physics_variables(): plasma_coulomb_log_electron_electron_profile = [] plasma_coulomb_log_electron_deuteron_profile = [] plasma_coulomb_log_electron_triton_profile = [] + plasma_coulomb_log_deuteron_triton_profile = [] freq_plasma_electron_profile = [] freq_plasma_larmor_toroidal_electron_profile = [] freq_plasma_larmor_toroidal_deuteron_profile = [] From d52244a7840917cbdce2f5cd4f7a055ac9aca443 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 16:27:56 +0000 Subject: [PATCH 02/20] Add deuteron-triton Coulomb logarithm calculations and plotting --- process/io/plot_proc.py | 13 ++++++++++ process/models/physics/physics.py | 41 +++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index d97d0bbcb..c95fcf707 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12771,6 +12771,11 @@ def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) ] + plasma_coulomb_log_deuteron_triton_profile = [ + mfile_data.data[f"plasma_coulomb_log_deuteron_triton_profile{i}"].get_scan(scan) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + axis.plot( np.linspace(0, 1, len(plasma_coulomb_log_electron_electron_profile)), plasma_coulomb_log_electron_electron_profile, @@ -12795,6 +12800,14 @@ def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): label=r"$ln \Lambda_{e-T}$", ) + axis.plot( + np.linspace(0, 1, len(plasma_coulomb_log_deuteron_triton_profile)), + plasma_coulomb_log_deuteron_triton_profile, + color="orange", + linestyle="-", + label=r"$ln \Lambda_{D-T}$", + ) + axis.set_ylabel("Coulomb Logarithm") axis.set_xlabel("$\\rho \\ [r/a]$") axis.grid(True, which="both", linestyle="--", alpha=0.5) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 91dd3bef0..6f717fe4d 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -10056,6 +10056,37 @@ def run(self): for i in range(len(physics_variables.len_plasma_debye_electron_profile)) ]) + physics_variables.plasma_coulomb_log_deuteron_triton_profile = np.array([ + self.calculate_coulomb_log_from_impact( + impact_param_max=physics_variables.len_plasma_debye_electron_profile[i], + impact_param_min=max( + self.calculate_classical_distance_of_closest_approach( + charge1=1, + charge2=1, + m_reduced=self.calculate_reduced_mass( + mass1=constants.TRITON_MASS, + mass2=constants.DEUTERON_MASS, + ), + vel_relative=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_triton_profile[i], + velocity_2=physics_variables.vel_plasma_deuteron_profile[i], + ), + ), + self.calculate_debroglie_wavelength( + mass=self.calculate_reduced_mass( + mass1=constants.TRITON_MASS, + mass2=constants.DEUTERON_MASS, + ), + velocity=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_triton_profile[i], + velocity_2=physics_variables.vel_plasma_deuteron_profile[i], + ), + ), + ), + ) + for i in range(len(physics_variables.len_plasma_debye_electron_profile)) + ]) + @staticmethod def calculate_debye_length( temp_plasma_species_kev: float | np.ndarray, @@ -10402,3 +10433,13 @@ def output_detailed_physics(self): f"(plasma_coulomb_log_electron_triton_profile{i})", physics_variables.plasma_coulomb_log_electron_triton_profile[i], ) + + for i in range( + len(physics_variables.plasma_coulomb_log_deuteron_triton_profile) + ): + po.ovarre( + self.mfile, + f"Deuteron-triton Coulomb log at point {i}", + f"(plasma_coulomb_log_deuteron_triton_profile{i})", + physics_variables.plasma_coulomb_log_deuteron_triton_profile[i], + ) From cf01f15f22af767a5118558f1ccab8902dc29379 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 16:31:30 +0000 Subject: [PATCH 03/20] Add deuteron plasma frequency profile to physics variables --- process/data_structure/physics_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 2dce159e2..844fb0675 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1335,6 +1335,9 @@ freq_plasma_electron_profile: list[float] = None """Electron plasma frequency profile (Hz)""" +freq_plasma_deuteron_profile: list[float] = None +"""Deuteron plasma frequency profile (Hz)""" + freq_plasma_larmor_toroidal_electron_profile: list[float] = None """Profile of electron Larmor frequency in plasma due to toroidal magnetic field (Hz)""" @@ -1671,6 +1674,7 @@ def init_physics_variables(): plasma_coulomb_log_electron_triton_profile, \ plasma_coulomb_log_deuteron_triton_profile, \ freq_plasma_electron_profile, \ + freq_plasma_deuteron_profile, \ freq_plasma_larmor_toroidal_electron_profile, \ freq_plasma_larmor_toroidal_deuteron_profile, \ freq_plasma_larmor_toroidal_triton_profile @@ -1941,6 +1945,7 @@ def init_physics_variables(): plasma_coulomb_log_electron_triton_profile = [] plasma_coulomb_log_deuteron_triton_profile = [] freq_plasma_electron_profile = [] + freq_plasma_deuteron_profile = [] freq_plasma_larmor_toroidal_electron_profile = [] freq_plasma_larmor_toroidal_deuteron_profile = [] freq_plasma_larmor_toroidal_triton_profile = [] From 1140eafbf600e1d78f37b30c97d57fc58c8641e1 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 17:00:41 +0000 Subject: [PATCH 04/20] Add alpha particle thermal velocity profile to physics variables --- process/data_structure/physics_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 844fb0675..ea427f73d 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1320,6 +1320,9 @@ vel_plasma_triton_profile: list[float] = None """Profile of triton thermal velocity in plasma (m/s)""" +vel_plasma_alpha_thermal_profile: list[float] = None +"""Profile of thermal alpha particle velocity in plasma (m/s)""" + plasma_coulomb_log_electron_electron_profile: list[float] = None """Profile of electron-electron Coulomb logarithm in plasma""" @@ -1669,6 +1672,7 @@ def init_physics_variables(): vel_plasma_electron_profile, \ vel_plasma_deuteron_profile, \ vel_plasma_triton_profile, \ + vel_plasma_alpha_thermal_profile, \ plasma_coulomb_log_electron_electron_profile, \ plasma_coulomb_log_electron_deuteron_profile, \ plasma_coulomb_log_electron_triton_profile, \ @@ -1940,6 +1944,7 @@ def init_physics_variables(): vel_plasma_electron_profile = [] vel_plasma_deuteron_profile = [] vel_plasma_triton_profile = [] + vel_plasma_alpha_thermal_profile = [] plasma_coulomb_log_electron_electron_profile = [] plasma_coulomb_log_electron_deuteron_profile = [] plasma_coulomb_log_electron_triton_profile = [] From 1a470de4c96b3aebf910034540756d8af992d170 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 17:05:23 +0000 Subject: [PATCH 05/20] Add alpha thermal velocity profile calculations and plotting --- process/io/plot_proc.py | 11 +++++++++++ process/models/physics/physics.py | 16 ++++++++++++++++ 2 files changed, 27 insertions(+) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index c95fcf707..f4f449e32 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12623,6 +12623,10 @@ def plot_velocity_profile(axis, mfile_data, scan): mfile_data.data[f"vel_plasma_triton_profile{i}"].get_scan(scan) for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) ] + vel_plasma_alpha_thermal_profile = [ + mfile_data.data[f"vel_plasma_alpha_thermal_profile{i}"].get_scan(scan) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] axis.plot( np.linspace(0, 1, len(vel_plasma_electron_profile)), @@ -12645,6 +12649,13 @@ def plot_velocity_profile(axis, mfile_data, scan): linestyle="-", label=r"$v_{T}$", ) + axis.plot( + np.linspace(0, 1, len(vel_plasma_alpha_thermal_profile)), + vel_plasma_alpha_thermal_profile, + color="orange", + linestyle="-", + label=r"$v_{\alpha,thermal}$", + ) axis.set_yscale("log") axis.set_ylabel("Velocity [m/s]") diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 6f717fe4d..f8e232616 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9921,6 +9921,15 @@ def run(self): ) ) + physics_variables.vel_plasma_alpha_thermal_profile = ( + self.calculate_relativistic_particle_speed( + e_kinetic=self.plasma_profile.teprofile.profile_y + * constants.KILOELECTRON_VOLT + * physics_variables.f_temp_plasma_ion_electron, + mass=constants.ALPHA_MASS, + ) + ) + # ============================ # Plasma frequencies # ============================ @@ -10364,6 +10373,13 @@ def output_detailed_physics(self): f"(vel_plasma_triton_profile{i})", physics_variables.vel_plasma_triton_profile[i], ) + for i in range(len(physics_variables.vel_plasma_alpha_thermal_profile)): + po.ovarre( + self.mfile, + f"Plasma alpha thermal velocity at point {i}", + f"(vel_plasma_alpha_thermal_profile{i})", + physics_variables.vel_plasma_alpha_thermal_profile[i], + ) po.osubhd(self.outfile, "Frequencies:") From fb8f337f80650861863745b63f633397b7c36149 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 17:10:43 +0000 Subject: [PATCH 06/20] Add electron-alpha Coulomb logarithm profile to physics variables --- process/data_structure/physics_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index ea427f73d..3566758c3 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1335,6 +1335,9 @@ plasma_coulomb_log_deuteron_triton_profile: list[float] = None """Profile of deuteron-triton Coulomb logarithm in plasma""" +plasma_coulomb_log_electron_alpha_thermal_profile: list[float] = None +"""Profile of electron-alpha Coulomb logarithm in plasma""" + freq_plasma_electron_profile: list[float] = None """Electron plasma frequency profile (Hz)""" @@ -1677,6 +1680,7 @@ def init_physics_variables(): plasma_coulomb_log_electron_deuteron_profile, \ plasma_coulomb_log_electron_triton_profile, \ plasma_coulomb_log_deuteron_triton_profile, \ + plasma_coulomb_log_electron_alpha_thermal_profile, \ freq_plasma_electron_profile, \ freq_plasma_deuteron_profile, \ freq_plasma_larmor_toroidal_electron_profile, \ @@ -1949,6 +1953,7 @@ def init_physics_variables(): plasma_coulomb_log_electron_deuteron_profile = [] plasma_coulomb_log_electron_triton_profile = [] plasma_coulomb_log_deuteron_triton_profile = [] + plasma_coulomb_log_electron_alpha_thermal_profile = [] freq_plasma_electron_profile = [] freq_plasma_deuteron_profile = [] freq_plasma_larmor_toroidal_electron_profile = [] From 098ee6f7f9dcab31107f9bb5bc1df5e470b40f1a Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 17:25:45 +0000 Subject: [PATCH 07/20] Add electron-alpha thermal Coulomb logarithm calculations and plotting --- process/io/plot_proc.py | 15 +++++++++++ process/models/physics/physics.py | 45 +++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index f4f449e32..7c65e0c6e 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12787,6 +12787,13 @@ def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) ] + plasma_coulomb_log_electron_alpha_thermal_profile = [ + mfile_data.data[ + f"plasma_coulomb_log_electron_alpha_thermal_profile{i}" + ].get_scan(scan) + for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) + ] + axis.plot( np.linspace(0, 1, len(plasma_coulomb_log_electron_electron_profile)), plasma_coulomb_log_electron_electron_profile, @@ -12819,6 +12826,14 @@ def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): label=r"$ln \Lambda_{D-T}$", ) + axis.plot( + np.linspace(0, 1, len(plasma_coulomb_log_electron_alpha_thermal_profile)), + plasma_coulomb_log_electron_alpha_thermal_profile, + color="purple", + linestyle="-", + label=r"$ln \Lambda_{e-\alpha,thermal}$", + ) + axis.set_ylabel("Coulomb Logarithm") axis.set_xlabel("$\\rho \\ [r/a]$") axis.grid(True, which="both", linestyle="--", alpha=0.5) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index f8e232616..de578223a 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -10096,6 +10096,41 @@ def run(self): for i in range(len(physics_variables.len_plasma_debye_electron_profile)) ]) + physics_variables.plasma_coulomb_log_electron_alpha_thermal_profile = np.array([ + self.calculate_coulomb_log_from_impact( + impact_param_max=physics_variables.len_plasma_debye_electron_profile[i], + impact_param_min=max( + self.calculate_classical_distance_of_closest_approach( + charge1=1, + charge2=2, + m_reduced=self.calculate_reduced_mass( + mass1=constants.ELECTRON_MASS, + mass2=constants.ALPHA_MASS, + ), + vel_relative=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_electron_profile[i], + velocity_2=physics_variables.vel_plasma_alpha_thermal_profile[ + i + ], + ), + ), + self.calculate_debroglie_wavelength( + mass=self.calculate_reduced_mass( + mass1=constants.ELECTRON_MASS, + mass2=constants.ALPHA_MASS, + ), + velocity=self.calculate_average_relative_velocity( + velocity_1=physics_variables.vel_plasma_electron_profile[i], + velocity_2=physics_variables.vel_plasma_alpha_thermal_profile[ + i + ], + ), + ), + ), + ) + for i in range(len(physics_variables.len_plasma_debye_electron_profile)) + ]) + @staticmethod def calculate_debye_length( temp_plasma_species_kev: float | np.ndarray, @@ -10459,3 +10494,13 @@ def output_detailed_physics(self): f"(plasma_coulomb_log_deuteron_triton_profile{i})", physics_variables.plasma_coulomb_log_deuteron_triton_profile[i], ) + + for i in range( + len(physics_variables.plasma_coulomb_log_electron_alpha_thermal_profile) + ): + po.ovarre( + self.mfile, + f"Electron-alpha thermal Coulomb log at point {i}", + f"(plasma_coulomb_log_electron_alpha_thermal_profile{i})", + physics_variables.plasma_coulomb_log_electron_alpha_thermal_profile[i], + ) From 06fd11d83af72ad8c17403aab5fc73ca641398e9 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 17:31:33 +0000 Subject: [PATCH 08/20] Add alpha particle birth velocity profile to physics variables --- process/data_structure/physics_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 3566758c3..ae9398621 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1323,6 +1323,9 @@ vel_plasma_alpha_thermal_profile: list[float] = None """Profile of thermal alpha particle velocity in plasma (m/s)""" +vel_plasma_alpha_birth: float = None +"""Birth velocity of alpha particles in plasma (m/s)""" + plasma_coulomb_log_electron_electron_profile: list[float] = None """Profile of electron-electron Coulomb logarithm in plasma""" @@ -1676,6 +1679,7 @@ def init_physics_variables(): vel_plasma_deuteron_profile, \ vel_plasma_triton_profile, \ vel_plasma_alpha_thermal_profile, \ + vel_plasma_alpha_birth, \ plasma_coulomb_log_electron_electron_profile, \ plasma_coulomb_log_electron_deuteron_profile, \ plasma_coulomb_log_electron_triton_profile, \ @@ -1949,6 +1953,7 @@ def init_physics_variables(): vel_plasma_deuteron_profile = [] vel_plasma_triton_profile = [] vel_plasma_alpha_thermal_profile = [] + vel_plasma_alpha_birth = 0.0 plasma_coulomb_log_electron_electron_profile = [] plasma_coulomb_log_electron_deuteron_profile = [] plasma_coulomb_log_electron_triton_profile = [] From 95bdfd199fa5f47e5eabe55eb594e049f23f493b Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 17:36:45 +0000 Subject: [PATCH 09/20] Add alpha particle birth velocity calculation and plotting --- process/io/plot_proc.py | 9 +++++++++ process/models/physics/physics.py | 13 +++++++++++++ 2 files changed, 22 insertions(+) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 7c65e0c6e..0f0f5ae01 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12628,6 +12628,8 @@ def plot_velocity_profile(axis, mfile_data, scan): for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) ] + vel_plasma_alpha_birth = mfile_data.data["vel_plasma_alpha_birth"].get_scan(scan) + axis.plot( np.linspace(0, 1, len(vel_plasma_electron_profile)), vel_plasma_electron_profile, @@ -12656,6 +12658,13 @@ def plot_velocity_profile(axis, mfile_data, scan): linestyle="-", label=r"$v_{\alpha,thermal}$", ) + axis.axhline( + vel_plasma_alpha_birth, + color="red", + linestyle="--", + linewidth=1.5, + label=r"$v_{\alpha,birth}$", + ) axis.set_yscale("log") axis.set_ylabel("Velocity [m/s]") diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index de578223a..2acfc8c0d 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9930,6 +9930,13 @@ def run(self): ) ) + physics_variables.vel_plasma_alpha_birth = ( + self.calculate_relativistic_particle_speed( + e_kinetic=constants.DT_ALPHA_ENERGY, + mass=constants.ALPHA_MASS, + ) + ) + # ============================ # Plasma frequencies # ============================ @@ -10415,6 +10422,12 @@ def output_detailed_physics(self): f"(vel_plasma_alpha_thermal_profile{i})", physics_variables.vel_plasma_alpha_thermal_profile[i], ) + po.ovarre( + self.outfile, + "Plasma alpha birth velocity (m/s)", + "(vel_plasma_alpha_birth)", + physics_variables.vel_plasma_alpha_birth, + ) po.osubhd(self.outfile, "Frequencies:") From af44da4494f4a6b701ef5d4b7277529c291a8c27 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 17:55:01 +0000 Subject: [PATCH 10/20] Add deuteron Larmor radius profile to physics variables --- process/data_structure/physics_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index ae9398621..6313b459e 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1308,6 +1308,9 @@ len_plasma_debye_electron_profile: list[float] = None """Profile of electron Debye length in plasma (m)""" +radius_plasma_deuteron_larmor_isotropic_profile: list[float] = None +"""Profile of deuteron Larmor radius in plasma, assuming equal speeds in all directions (m)""" + len_plasma_debye_electron_vol_avg: float = None """Volume averaged electron Debye length in plasma (m)""" @@ -1674,6 +1677,7 @@ def init_physics_variables(): n_charge_plasma_effective_mass_weighted_vol_avg, \ j_plasma_bootstrap_sauter_profile, \ len_plasma_debye_electron_profile, \ + radius_plasma_deuteron_larmor_isotropic_profile, \ len_plasma_debye_electron_vol_avg, \ vel_plasma_electron_profile, \ vel_plasma_deuteron_profile, \ @@ -1948,6 +1952,7 @@ def init_physics_variables(): n_charge_plasma_effective_profile = [] n_charge_plasma_effective_mass_weighted_vol_avg = 0.0 len_plasma_debye_electron_profile = [] + radius_plasma_deuteron_larmor_isotropic_profile = [] len_plasma_debye_electron_vol_avg = 0.0 vel_plasma_electron_profile = [] vel_plasma_deuteron_profile = [] From 539edc2ae2e90f0454f4e725eb45c0d0593a370e Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 17:59:36 +0000 Subject: [PATCH 11/20] Add Larmor radius calculation method to DetailedPhysics class --- process/models/physics/physics.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 2acfc8c0d..56afd8c06 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -10332,6 +10332,22 @@ def calculate_larmor_frequency( 2 * np.pi * m_particle ) + @staticmethod + def calculate_larmor_radius( + vel_perp: float | np.ndarray, + freq_larmor: float, + ) -> float | np.ndarray: + """ + Calculate the Larmor radius for a particle species. + :param vel_perp: Perpendicular velocity of the particle to the magnetic field (m/s). + :type vel_perp: float | np.ndarray + :param freq_larmor: Larmor frequency of the particle (Hz). + :type freq_larmor: float + :returns: Larmor radius in meters. + :rtype: float | np.ndarray + """ + return vel_perp / (freq_larmor) + @staticmethod def calculate_reduced_mass(mass1: float, mass2: float) -> float: """ From 71b41e2ae982177a6f7429779435a5df2b73836f Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 19:00:57 +0000 Subject: [PATCH 12/20] Add toroidal Larmor radius calculation and plotting to physics variables --- process/data_structure/physics_variables.py | 8 +- process/io/plot_proc.py | 93 ++++++++++++++------- process/models/physics/physics.py | 31 +++++++ 3 files changed, 98 insertions(+), 34 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 6313b459e..978307244 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1308,8 +1308,8 @@ len_plasma_debye_electron_profile: list[float] = None """Profile of electron Debye length in plasma (m)""" -radius_plasma_deuteron_larmor_isotropic_profile: list[float] = None -"""Profile of deuteron Larmor radius in plasma, assuming equal speeds in all directions (m)""" +radius_plasma_deuteron_toroidal_larmor_isotropic_profile: list[float] = None +"""Profile of deuteron toroidal Larmor radius in plasma, assuming equal speeds in all directions (m)""" len_plasma_debye_electron_vol_avg: float = None """Volume averaged electron Debye length in plasma (m)""" @@ -1677,7 +1677,7 @@ def init_physics_variables(): n_charge_plasma_effective_mass_weighted_vol_avg, \ j_plasma_bootstrap_sauter_profile, \ len_plasma_debye_electron_profile, \ - radius_plasma_deuteron_larmor_isotropic_profile, \ + radius_plasma_deuteron_toroidal_larmor_isotropic_profile, \ len_plasma_debye_electron_vol_avg, \ vel_plasma_electron_profile, \ vel_plasma_deuteron_profile, \ @@ -1952,7 +1952,7 @@ def init_physics_variables(): n_charge_plasma_effective_profile = [] n_charge_plasma_effective_mass_weighted_vol_avg = 0.0 len_plasma_debye_electron_profile = [] - radius_plasma_deuteron_larmor_isotropic_profile = [] + radius_plasma_deuteron_toroidal_larmor_isotropic_profile = [] len_plasma_debye_electron_vol_avg = 0.0 vel_plasma_electron_profile = [] vel_plasma_deuteron_profile = [] diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 0f0f5ae01..7bcbbcfa5 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12560,6 +12560,37 @@ def plot_ebw_ecrh_coupling_graph(axis: plt.Axes, mfile: mf.MFile, scan: int): axis.minorticks_on() +def plot_larmor_radius_profile(axis: plt.Axes, mfile_data: mf.MFile, scan: int): + """Plot the Larmor radius profile on the given axis.""" + radius_plasma_deuteron_larmor_profile = [ + mfile_data.data[ + f"radius_plasma_deuteron_toroidal_larmor_isotropic_profile{i}" + ].get_scan(scan) + for i in range( + 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + radius_plasma_deuteron_larmor_profile_mm = [ + radius * 1e3 for radius in radius_plasma_deuteron_larmor_profile + ] + + axis.plot( + np.linspace(-1, 1, len(radius_plasma_deuteron_larmor_profile_mm)), + radius_plasma_deuteron_larmor_profile_mm, + color="red", + linestyle="-", + label=r"$\rho_{Larmor,toroidal,D}$", + ) + + axis.set_ylabel(r"Larmor Radii [mm]") + axis.set_title(r"Larmor Radii ($v_{\perp}^2 = 2v_{th}^2$)") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + + def plot_debye_length_profile(axis: plt.Axes, mfile_data: mf.MFile, scan: int): """Plot the Debye length profile on the given axis. @@ -13310,9 +13341,11 @@ def main_plot( plot_plasma_coloumb_logarithms(figs[15].add_subplot(231), m_file, scan) + plot_larmor_radius_profile(figs[16].add_subplot(311), m_file, scan) + # Plot poloidal cross-section poloidal_cross_section( - figs[16].add_subplot(121, aspect="equal"), + figs[17].add_subplot(121, aspect="equal"), m_file, scan, demo_ranges, @@ -13322,7 +13355,7 @@ def main_plot( # Plot toroidal cross-section toroidal_cross_section( - figs[16].add_subplot(122, aspect="equal"), + figs[17].add_subplot(122, aspect="equal"), m_file, scan, demo_ranges, @@ -13330,19 +13363,19 @@ def main_plot( ) # Plot color key - ax17 = figs[16].add_subplot(222) + ax17 = figs[17].add_subplot(222) ax17.set_position([0.5, 0.5, 0.5, 0.5]) color_key(ax17, m_file, scan, colour_scheme) - ax18 = figs[17].add_subplot(211) + ax18 = figs[18].add_subplot(211) ax18.set_position([0.1, 0.33, 0.8, 0.6]) plot_radial_build(ax18, m_file, colour_scheme) # Make each axes smaller vertically to leave room for the legend - ax185 = figs[18].add_subplot(211) + ax185 = figs[19].add_subplot(211) ax185.set_position([0.1, 0.61, 0.8, 0.32]) - ax18b = figs[18].add_subplot(212) + ax18b = figs[19].add_subplot(212) ax18b.set_position([0.1, 0.13, 0.8, 0.32]) plot_upper_vertical_build(ax185, m_file, colour_scheme) plot_lower_vertical_build(ax18b, m_file, colour_scheme) @@ -13350,52 +13383,52 @@ def main_plot( # Can only plot WP and turn structure if superconducting coil at the moment if m_file.get("i_tf_sup", scan=scan) == 1: # TF coil with WP - ax19 = figs[19].add_subplot(221, aspect="equal") + ax19 = figs[20].add_subplot(221, aspect="equal") ax19.set_position([ 0.025, 0.45, 0.5, 0.5, ]) # Half height, a bit wider, top left - plot_superconducting_tf_wp(ax19, m_file, scan, figs[19]) + plot_superconducting_tf_wp(ax19, m_file, scan, figs[20]) # TF coil turn structure - ax20 = figs[20].add_subplot(325, aspect="equal") + ax20 = figs[21].add_subplot(325, aspect="equal") ax20.set_position([0.025, 0.5, 0.4, 0.4]) - plot_tf_cable_in_conduit_turn(ax20, figs[20], m_file, scan) - plot_205 = figs[20].add_subplot(223, aspect="equal") + plot_tf_cable_in_conduit_turn(ax20, figs[21], m_file, scan) + plot_205 = figs[21].add_subplot(223, aspect="equal") plot_205.set_position([0.075, 0.1, 0.3, 0.3]) - plot_cable_in_conduit_cable(plot_205, figs[20], m_file, scan) + plot_cable_in_conduit_cable(plot_205, figs[21], m_file, scan) else: - ax19 = figs[19].add_subplot(211, aspect="equal") + ax19 = figs[20].add_subplot(211, aspect="equal") ax19.set_position([0.06, 0.55, 0.675, 0.4]) - plot_resistive_tf_wp(ax19, m_file, scan, figs[19]) - plot_resistive_tf_info(ax19, m_file, scan, figs[19]) + plot_resistive_tf_wp(ax19, m_file, scan, figs[20]) + plot_resistive_tf_info(ax19, m_file, scan, figs[20]) plot_tf_coil_structure( - figs[21].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme + figs[22].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme ) - plot_plasma_outboard_toroidal_ripple_map(figs[22], m_file, scan) + plot_plasma_outboard_toroidal_ripple_map(figs[23], m_file, scan) - plot_tf_stress(figs[23].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) + plot_tf_stress(figs[24].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) - plot_current_profiles_over_time(figs[24].add_subplot(111), m_file, scan) + plot_current_profiles_over_time(figs[25].add_subplot(111), m_file, scan) plot_cs_coil_structure( - figs[25].add_subplot(121, aspect="equal"), figs[25], m_file, scan + figs[26].add_subplot(121, aspect="equal"), figs[26], m_file, scan ) plot_cs_turn_structure( - figs[25].add_subplot(224, aspect="equal"), figs[25], m_file, scan + figs[26].add_subplot(224, aspect="equal"), figs[26], m_file, scan ) plot_first_wall_top_down_cross_section( - figs[26].add_subplot(221, aspect="equal"), m_file, scan + figs[27].add_subplot(221, aspect="equal"), m_file, scan ) - plot_first_wall_poloidal_cross_section(figs[26].add_subplot(122), m_file, scan) - plot_fw_90_deg_pipe_bend(figs[26].add_subplot(337), m_file, scan) + plot_first_wall_poloidal_cross_section(figs[27].add_subplot(122), m_file, scan) + plot_fw_90_deg_pipe_bend(figs[27].add_subplot(337), m_file, scan) - plot_blkt_pipe_bends(figs[27], m_file, scan) - ax_blanket = figs[27].add_subplot(122, aspect="equal") + plot_blkt_pipe_bends(figs[28], m_file, scan) + ax_blanket = figs[28].add_subplot(122, aspect="equal") plot_blanket(ax_blanket, m_file, scan, radial_build, colour_scheme) plot_firstwall(ax_blanket, m_file, scan, radial_build, colour_scheme) ax_blanket.set_xlabel("Radial position [m]") @@ -13438,13 +13471,13 @@ def main_plot( ) plot_main_power_flow( - figs[28].add_subplot(111, aspect="equal"), m_file, scan, figs[28] + figs[29].add_subplot(111, aspect="equal"), m_file, scan, figs[29] ) - ax24 = figs[29].add_subplot(111) + ax24 = figs[30].add_subplot(111) # set_position([left, bottom, width, height]) -> height ~ 0.66 => ~2/3 of page height ax24.set_position([0.08, 0.35, 0.84, 0.57]) - plot_system_power_profiles_over_time(ax24, m_file, scan, figs[29]) + plot_system_power_profiles_over_time(ax24, m_file, scan, figs[30]) def create_thickness_builds(m_file, scan: int): @@ -13521,7 +13554,7 @@ def main(args=None): # create main plot # Increase range when adding new page - pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(30)] + pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(31)] # run main_plot main_plot( diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 56afd8c06..a2b7c4bd6 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9975,6 +9975,21 @@ def run(self): ) ) + # ============================ + # Larmor radii + # ============================ + + # Since isotropic (v⟂)² = 2(v)² for a Maxwellian distribution, + # we can use the total velocity to calculate the Larmor radius for an isotropic profile + physics_variables.radius_plasma_deuteron_toroidal_larmor_isotropic_profile = self.calculate_larmor_radius( + vel_perp=np.sqrt(2) + * np.concatenate([ + physics_variables.vel_plasma_deuteron_profile[::-1], + physics_variables.vel_plasma_deuteron_profile, + ]), + freq_larmor=physics_variables.freq_plasma_larmor_toroidal_deuteron_profile, + ) + # ============================ # Coulomb logarithm # ============================ @@ -10408,6 +10423,22 @@ def output_detailed_physics(self): physics_variables.len_plasma_debye_electron_profile[i], ) + po.osubhd(self.outfile, "Larmor radii:") + + for i in range( + len( + physics_variables.radius_plasma_deuteron_toroidal_larmor_isotropic_profile + ) + ): + po.ovarre( + self.mfile, + f"Plasma deuteron isotropic Larmor radius at point {i}", + f"(radius_plasma_deuteron_toroidal_larmor_isotropic_profile{i})", + physics_variables.radius_plasma_deuteron_toroidal_larmor_isotropic_profile[ + i + ], + ) + po.osubhd(self.outfile, "Velocities:") for i in range(len(physics_variables.vel_plasma_electron_profile)): From c5754d9483644c5b4b39a8db026658831ee40be7 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 19:06:20 +0000 Subject: [PATCH 13/20] Add triton Larmor radius profile to physics variables --- process/data_structure/physics_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 978307244..29c1e8916 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1311,6 +1311,9 @@ radius_plasma_deuteron_toroidal_larmor_isotropic_profile: list[float] = None """Profile of deuteron toroidal Larmor radius in plasma, assuming equal speeds in all directions (m)""" +radius_plasma_triton_toroidal_larmor_isotropic_profile: list[float] = None +"""Profile of triton toroidal Larmor radius in plasma, assuming equal speeds in all directions (m)""" + len_plasma_debye_electron_vol_avg: float = None """Volume averaged electron Debye length in plasma (m)""" @@ -1678,6 +1681,7 @@ def init_physics_variables(): j_plasma_bootstrap_sauter_profile, \ len_plasma_debye_electron_profile, \ radius_plasma_deuteron_toroidal_larmor_isotropic_profile, \ + radius_plasma_triton_toroidal_larmor_isotropic_profile, \ len_plasma_debye_electron_vol_avg, \ vel_plasma_electron_profile, \ vel_plasma_deuteron_profile, \ @@ -1953,6 +1957,7 @@ def init_physics_variables(): n_charge_plasma_effective_mass_weighted_vol_avg = 0.0 len_plasma_debye_electron_profile = [] radius_plasma_deuteron_toroidal_larmor_isotropic_profile = [] + radius_plasma_triton_toroidal_larmor_isotropic_profile = [] len_plasma_debye_electron_vol_avg = 0.0 vel_plasma_electron_profile = [] vel_plasma_deuteron_profile = [] From 75da07603d5de9ebb507d629268cf8c04eacf4bd Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sat, 14 Feb 2026 19:10:54 +0000 Subject: [PATCH 14/20] Add triton isotropic Larmor radius calculation to DetailedPhysics class --- process/io/plot_proc.py | 23 +++++++++++++++- process/models/physics/physics.py | 44 ++++++++++++++++++++++++++----- 2 files changed, 60 insertions(+), 7 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 7bcbbcfa5..d818b80d3 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12571,10 +12571,23 @@ def plot_larmor_radius_profile(axis: plt.Axes, mfile_data: mf.MFile, scan: int): ) ] + radius_plasma_triton_larmor_profile = [ + mfile_data.data[ + f"radius_plasma_triton_toroidal_larmor_isotropic_profile{i}" + ].get_scan(scan) + for i in range( + 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + radius_plasma_deuteron_larmor_profile_mm = [ radius * 1e3 for radius in radius_plasma_deuteron_larmor_profile ] + radius_plasma_triton_larmor_profile_mm = [ + radius * 1e3 for radius in radius_plasma_triton_larmor_profile + ] + axis.plot( np.linspace(-1, 1, len(radius_plasma_deuteron_larmor_profile_mm)), radius_plasma_deuteron_larmor_profile_mm, @@ -12583,8 +12596,16 @@ def plot_larmor_radius_profile(axis: plt.Axes, mfile_data: mf.MFile, scan: int): label=r"$\rho_{Larmor,toroidal,D}$", ) + axis.plot( + np.linspace(-1, 1, len(radius_plasma_triton_larmor_profile_mm)), + radius_plasma_triton_larmor_profile_mm, + color="green", + linestyle="-", + label=r"$\rho_{Larmor,toroidal,T}$", + ) + axis.set_ylabel(r"Larmor Radii [mm]") - axis.set_title(r"Larmor Radii ($v_{\perp}^2 = 2v_{th}^2$)") + axis.set_title(r" Toroidal Larmor Radii ($v_{\perp}^2 = 2v_{th}^2$)") axis.set_xlabel("$\\rho \\ [r/a]$") axis.grid(True, which="both", linestyle="--", alpha=0.5) axis.minorticks_on() diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index a2b7c4bd6..3c557da71 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9982,12 +9982,33 @@ def run(self): # Since isotropic (v⟂)² = 2(v)² for a Maxwellian distribution, # we can use the total velocity to calculate the Larmor radius for an isotropic profile physics_variables.radius_plasma_deuteron_toroidal_larmor_isotropic_profile = self.calculate_larmor_radius( - vel_perp=np.sqrt(2) - * np.concatenate([ - physics_variables.vel_plasma_deuteron_profile[::-1], - physics_variables.vel_plasma_deuteron_profile, - ]), - freq_larmor=physics_variables.freq_plasma_larmor_toroidal_deuteron_profile, + vel_perp=np.sqrt( + 2 + * np.concatenate([ + physics_variables.vel_plasma_deuteron_profile[::-1], + physics_variables.vel_plasma_deuteron_profile, + ]) + ** 2 + ), + freq_larmor=physics_variables.freq_plasma_larmor_toroidal_deuteron_profile + * (2 * np.pi), + ) + + # Since isotropic (v⟂)² = 2(v)² for a Maxwellian distribution, + # we can use the total velocity to calculate the Larmor radius for an isotropic profile + physics_variables.radius_plasma_triton_toroidal_larmor_isotropic_profile = ( + self.calculate_larmor_radius( + vel_perp=np.sqrt( + 2 + * np.concatenate([ + physics_variables.vel_plasma_triton_profile[::-1], + physics_variables.vel_plasma_triton_profile, + ]) + ** 2 + ), + freq_larmor=physics_variables.freq_plasma_larmor_toroidal_triton_profile + * (2 * np.pi), + ) ) # ============================ @@ -10438,6 +10459,17 @@ def output_detailed_physics(self): i ], ) + for i in range( + len(physics_variables.radius_plasma_triton_toroidal_larmor_isotropic_profile) + ): + po.ovarre( + self.mfile, + f"Plasma triton isotropic Larmor radius at point {i}", + f"(radius_plasma_triton_toroidal_larmor_isotropic_profile{i})", + physics_variables.radius_plasma_triton_toroidal_larmor_isotropic_profile[ + i + ], + ) po.osubhd(self.outfile, "Velocities:") From 4cae410f2828827097460eac994028ab46dc8b2e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Feb 2026 11:50:20 +0000 Subject: [PATCH 15/20] Add secondary x-axis for radius in electron frequency profile plot --- process/io/plot_proc.py | 35 +++++++++++++++++++++++++++-------- 1 file changed, 27 insertions(+), 8 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index d818b80d3..ed5dea66a 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12759,6 +12759,7 @@ def plot_electron_frequency_profile(axis, mfile_data, scan): linestyle="-", label=r"$f_{Larmor,toroidal,e}$", ) + x = np.linspace(0, 1, len(freq_plasma_electron_profile)) y = np.array(freq_plasma_electron_profile) / 1e9 # original curve @@ -12767,9 +12768,31 @@ def plot_electron_frequency_profile(axis, mfile_data, scan): axis.plot(-x, y, color="blue", linestyle="-", label="_nolegend_") axis.set_xlim(-1.025, 1.025) + axis.set_xlabel("$\\rho$ [r/a]") axis.set_ylabel("Frequency [GHz]") axis.grid(True, which="both", linestyle="--", alpha=0.5) + # Add secondary x-axis showing radius in metres below the primary axis + ax2 = axis.twiny() + rmajor = mfile_data.get("rmajor", scan=scan) + rminor = mfile_data.get("rminor", scan=scan) + + # Convert normalized radius to actual radius + # rho ranges from -1 to 1, which corresponds to r = rmajor - rminor to rmajor + rminor + rho_ticks = np.array([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1]) + r_ticks = rmajor + rho_ticks * rminor + + ax2.set_xticks(rho_ticks) + ax2.set_xticklabels([f"{r:.2f}" for r in r_ticks]) + ax2.set_xlabel("Radius [m]") + ax2.minorticks_on() + ax2.set_xlim(axis.get_xlim()) + + # Move secondary axis to the bottom + ax2.xaxis.set_ticks_position("bottom") + ax2.xaxis.set_label_position("bottom") + ax2.spines["bottom"].set_position(("outward", 40)) + axis.legend() @@ -13352,17 +13375,13 @@ def main_plot( plot_debye_length_profile(figs[15].add_subplot(232), m_file, scan) plot_velocity_profile(figs[15].add_subplot(233), m_file, scan) + plot_plasma_coloumb_logarithms(figs[15].add_subplot(231), m_file, scan) - ax_ion = figs[15].add_subplot(414) - ax_electron = figs[15].add_subplot(413, sharex=ax_ion) - - plot_electron_frequency_profile(ax_electron, m_file, scan) - - plot_ion_frequency_profile(ax_ion, m_file, scan) + plot_electron_frequency_profile(figs[15].add_subplot(212), m_file, scan) - plot_plasma_coloumb_logarithms(figs[15].add_subplot(231), m_file, scan) + plot_ion_frequency_profile(figs[16].add_subplot(411), m_file, scan) - plot_larmor_radius_profile(figs[16].add_subplot(311), m_file, scan) + plot_larmor_radius_profile(figs[16].add_subplot(313), m_file, scan) # Plot poloidal cross-section poloidal_cross_section( From 360cbb216a77a56a971bca642f4b84431db045b7 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Feb 2026 11:57:40 +0000 Subject: [PATCH 16/20] Add upper hybrid frequency calculation to DetailedPhysics class --- process/models/physics/physics.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 3c557da71..c627b8bf8 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9974,6 +9974,12 @@ def run(self): z_particle=1.0, ) ) + + # ============================ + # Upper hybrid frequencies + # ============================ + + # ============================ # Larmor radii @@ -10368,6 +10374,22 @@ def calculate_larmor_frequency( 2 * np.pi * m_particle ) + @staticmethod + def calculate_upper_hybrid_frequency( + freq_plasma: float | np.ndarray, freq_larmor: float | np.ndarray + ) -> float | np.ndarray: + """ + Calculate the upper hybrid frequency for a particle species. + + :param freq_plasma: Plasma frequency of the particle species (Hz). + :type freq_plasma: float | np.ndarray + :param freq_larmor: Larmor frequency of the particle species (Hz). + :type freq_larmor: float | np.ndarray + :returns: Upper hybrid frequency in Hz. + :rtype: float | np.ndarray + """ + return np.sqrt(freq_plasma**2 + freq_larmor**2) + @staticmethod def calculate_larmor_radius( vel_perp: float | np.ndarray, From 3fc9e5d412a0f3b76f3f38aec58d1219daf09061 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Feb 2026 11:58:37 +0000 Subject: [PATCH 17/20] Add upper hybrid frequency profile to physics variables --- process/data_structure/physics_variables.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 29c1e8916..f84d8e94b 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1362,6 +1362,9 @@ freq_plasma_larmor_toroidal_triton_profile: list[float] = None """Profile of triton Larmor frequency in plasma due to toroidal magnetic field (Hz)""" +freq_plasma_upper_hybrid_profile: list[float] = None +"""Profile of upper hybrid frequency in plasma (Hz)""" + def init_physics_module(): """Initialise the physics module""" @@ -1697,7 +1700,8 @@ def init_physics_variables(): freq_plasma_deuteron_profile, \ freq_plasma_larmor_toroidal_electron_profile, \ freq_plasma_larmor_toroidal_deuteron_profile, \ - freq_plasma_larmor_toroidal_triton_profile + freq_plasma_larmor_toroidal_triton_profile, \ + freq_plasma_upper_hybrid_profile m_beam_amu = 0.0 m_fuel_amu = 0.0 @@ -1974,3 +1978,4 @@ def init_physics_variables(): freq_plasma_larmor_toroidal_electron_profile = [] freq_plasma_larmor_toroidal_deuteron_profile = [] freq_plasma_larmor_toroidal_triton_profile = [] + freq_plasma_upper_hybrid_profile = [] From 2253daeebf4b4c848a09de5aad8d4fef4ad81f8b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Feb 2026 12:06:23 +0000 Subject: [PATCH 18/20] Add upper hybrid frequency profile calculation to DetailedPhysics class and plot --- process/io/plot_proc.py | 16 ++++++++++++++++ process/models/physics/physics.py | 20 +++++++++++++++++--- 2 files changed, 33 insertions(+), 3 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index ed5dea66a..11911ab01 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12752,6 +12752,13 @@ def plot_electron_frequency_profile(axis, mfile_data, scan): ) ] + freq_plasma_upper_hybrid_electron_profile = [ + mfile_data.data[f"freq_plasma_upper_hybrid_profile{i}"].get_scan(scan) + for i in range( + 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + axis.plot( np.linspace(-1, 1, len(freq_plasma_larmor_toroidal_electron_profile)), np.array(freq_plasma_larmor_toroidal_electron_profile) / 1e9, @@ -12766,6 +12773,15 @@ def plot_electron_frequency_profile(axis, mfile_data, scan): axis.plot(x, y, color="blue", linestyle="-", label=r"$\omega_{p,e}$") # mirrored across the y-axis (drawn at negative rho) axis.plot(-x, y, color="blue", linestyle="-", label="_nolegend_") + + axis.plot( + np.linspace(-1, 1, len(freq_plasma_upper_hybrid_electron_profile)), + np.array(freq_plasma_upper_hybrid_electron_profile) / 1e9, + color="purple", + linestyle="-", + label=r"$\omega_{UH,e}$", + ) + axis.set_xlim(-1.025, 1.025) axis.set_xlabel("$\\rho$ [r/a]") diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index c627b8bf8..d48afc783 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9974,12 +9974,18 @@ def run(self): z_particle=1.0, ) ) - + # ============================ # Upper hybrid frequencies # ============================ - + physics_variables.freq_plasma_upper_hybrid_profile = self.calculate_upper_hybrid_frequency( + freq_plasma=np.concatenate([ + physics_variables.freq_plasma_electron_profile[::-1], + physics_variables.freq_plasma_electron_profile, + ]), + freq_larmor=physics_variables.freq_plasma_larmor_toroidal_electron_profile, + ) # ============================ # Larmor radii @@ -10389,7 +10395,7 @@ def calculate_upper_hybrid_frequency( :rtype: float | np.ndarray """ return np.sqrt(freq_plasma**2 + freq_larmor**2) - + @staticmethod def calculate_larmor_radius( vel_perp: float | np.ndarray, @@ -10567,6 +10573,14 @@ def output_detailed_physics(self): physics_variables.freq_plasma_larmor_toroidal_triton_profile[i], ) + for i in range(len(physics_variables.freq_plasma_upper_hybrid_profile)): + po.ovarre( + self.mfile, + f"Plasma upper hybrid frequency at point {i}", + f"(freq_plasma_upper_hybrid_profile{i})", + physics_variables.freq_plasma_upper_hybrid_profile[i], + ) + po.osubhd(self.outfile, "Coulomb Logarithms:") for i in range( From 06edc986e4579b068349c7aee7cc9e8b728b723e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Feb 2026 13:26:15 +0000 Subject: [PATCH 19/20] Refactor plot colors and labels in velocity and frequency profiles for clarity --- process/io/plot_proc.py | 37 ++++++++++++++++++++++++++++--------- 1 file changed, 28 insertions(+), 9 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 11911ab01..288e63405 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -12692,7 +12692,7 @@ def plot_velocity_profile(axis, mfile_data, scan): axis.plot( np.linspace(0, 1, len(vel_plasma_deuteron_profile)), vel_plasma_deuteron_profile, - color="red", + color="pink", linestyle="-", label=r"$v_{D}$", ) @@ -12706,7 +12706,7 @@ def plot_velocity_profile(axis, mfile_data, scan): axis.plot( np.linspace(0, 1, len(vel_plasma_alpha_thermal_profile)), vel_plasma_alpha_thermal_profile, - color="orange", + color="red", linestyle="-", label=r"$v_{\alpha,thermal}$", ) @@ -12764,13 +12764,31 @@ def plot_electron_frequency_profile(axis, mfile_data, scan): np.array(freq_plasma_larmor_toroidal_electron_profile) / 1e9, color="red", linestyle="-", - label=r"$f_{Larmor,toroidal,e}$", + label=r"$f_{Larmor,toroidal,e}$ | Fundamental", + ) + + axis.plot( + np.linspace(-1, 1, len(freq_plasma_larmor_toroidal_electron_profile)), + 2 * np.array(freq_plasma_larmor_toroidal_electron_profile) / 1e9, + color="red", + linestyle="--", + label=r"$f_{Larmor,toroidal,e}$ | 2nd harmonic", + ) + + axis.plot( + np.linspace(-1, 1, len(freq_plasma_larmor_toroidal_electron_profile)), + 3 * np.array(freq_plasma_larmor_toroidal_electron_profile) / 1e9, + color="red", + linestyle=":", + label=r"$f_{Larmor,toroidal,e}$ | 3rd harmonic", ) x = np.linspace(0, 1, len(freq_plasma_electron_profile)) y = np.array(freq_plasma_electron_profile) / 1e9 # original curve - axis.plot(x, y, color="blue", linestyle="-", label=r"$\omega_{p,e}$") + axis.plot( + x, y, color="blue", linestyle="-", label=r"$\omega_{p,e}$ | Plasma Frequency" + ) # mirrored across the y-axis (drawn at negative rho) axis.plot(-x, y, color="blue", linestyle="-", label="_nolegend_") @@ -12779,10 +12797,11 @@ def plot_electron_frequency_profile(axis, mfile_data, scan): np.array(freq_plasma_upper_hybrid_electron_profile) / 1e9, color="purple", linestyle="-", - label=r"$\omega_{UH,e}$", + label=r"$\omega_{UH,e}$ | Upper Hybrid", ) axis.set_xlim(-1.025, 1.025) + axis.set_ylim(None, max(freq_plasma_larmor_toroidal_electron_profile) / 1e9 * 1.6) axis.set_xlabel("$\\rho$ [r/a]") axis.set_ylabel("Frequency [GHz]") @@ -12807,7 +12826,7 @@ def plot_electron_frequency_profile(axis, mfile_data, scan): # Move secondary axis to the bottom ax2.xaxis.set_ticks_position("bottom") ax2.xaxis.set_label_position("bottom") - ax2.spines["bottom"].set_position(("outward", 40)) + ax2.spines["bottom"].set_position(("outward", 30)) axis.legend() @@ -12905,7 +12924,7 @@ def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): axis.plot( np.linspace(0, 1, len(plasma_coulomb_log_electron_deuteron_profile)), plasma_coulomb_log_electron_deuteron_profile, - color="red", + color="pink", linestyle="-", label=r"$ln \Lambda_{e-D}$", ) @@ -12929,7 +12948,7 @@ def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): axis.plot( np.linspace(0, 1, len(plasma_coulomb_log_electron_alpha_thermal_profile)), plasma_coulomb_log_electron_alpha_thermal_profile, - color="purple", + color="red", linestyle="-", label=r"$ln \Lambda_{e-\alpha,thermal}$", ) @@ -13395,7 +13414,7 @@ def main_plot( plot_electron_frequency_profile(figs[15].add_subplot(212), m_file, scan) - plot_ion_frequency_profile(figs[16].add_subplot(411), m_file, scan) + plot_ion_frequency_profile(figs[16].add_subplot(311), m_file, scan) plot_larmor_radius_profile(figs[16].add_subplot(313), m_file, scan) From 02682d19751628b8130436d73f5aaf0221d22768 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sat, 21 Feb 2026 15:59:00 +0000 Subject: [PATCH 20/20] Fix toroidal field profile size in detailed physics test for accurate Larmor frequency calculation --- tests/unit/test_physics.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 93273518b..b09401ccb 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -3617,7 +3617,7 @@ def test_detailed_physics_run_computes_profiles(): ) # toroidal field profile for larmor frequency calc physics_variables.b_plasma_toroidal_profile = ( - np.ones_like(plasma.neprofile.profile_y) * 5.0 + np.ones(2 * len(plasma.neprofile.profile_y)) * 5.0 ) dp = DetailedPhysics(plasma) @@ -3637,7 +3637,8 @@ def test_detailed_physics_run_computes_profiles(): assert np.all(np.isfinite(physics_variables.freq_plasma_electron_profile)) assert ( - np.shape(physics_variables.freq_plasma_larmor_toroidal_electron_profile)[0] == n + np.shape(physics_variables.freq_plasma_larmor_toroidal_electron_profile)[0] + == n * 2 ) assert np.all( np.isfinite(physics_variables.freq_plasma_larmor_toroidal_electron_profile)