diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 79bc79024..edd06e8b0 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -70,6 +70,7 @@ ) from process.models.physics.impurity_radiation import read_impurity_file from process.models.physics.l_h_transition import PlasmaConfinementTransitionModel +from process.models.physics.physics import DetailedPhysics from process.models.physics.plasma_current import PlasmaCurrent, PlasmaCurrentModel from process.models.tfcoil.base import TFCoilShapeModel from process.models.tfcoil.superconducting import SUPERCONDUCTING_TF_TYPES @@ -12769,406 +12770,6 @@ 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_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, - color="red", - linestyle="-", - 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" 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() - 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. - - Parameters - ---------- - axis : - - mfile_data : - - scan : - - """ - len_plasma_debye_electron_profile = [ - mfile_data.data[f"len_plasma_debye_electron_profile{i}"].get_scan(scan) - for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) - ] - - # Convert to micrometres (1e-6 m) - len_plasma_debye_electron_profile_um = [ - length * 1e6 for length in len_plasma_debye_electron_profile - ] - - axis.plot( - np.linspace(0, 1, len(len_plasma_debye_electron_profile_um)), - len_plasma_debye_electron_profile_um, - color="blue", - linestyle="-", - label=r"$\lambda_{Debye,e}$", - ) - - axis.set_ylabel(r"Debye Length [$\mu$m]") - - axis.set_xlabel("$\\rho \\ [r/a]$") - axis.grid(True, which="both", linestyle="--", alpha=0.5) - axis.set_xlim([0, 1.025]) - axis.minorticks_on() - axis.legend() - - -def plot_velocity_profile(axis, mfile_data, scan): - """Plot the electron thermal velocity profile on the given axis. - - Parameters - ---------- - axis : - - mfile_data : - - scan : - - """ - vel_plasma_electron_profile = [ - mfile_data.data[f"vel_plasma_electron_profile{i}"].get_scan(scan) - for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) - ] - vel_plasma_deuteron_profile = [ - mfile_data.data[f"vel_plasma_deuteron_profile{i}"].get_scan(scan) - for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) - ] - vel_plasma_triton_profile = [ - 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))) - ] - - 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, - color="blue", - linestyle="-", - label=r"$v_{e}$", - ) - axis.plot( - np.linspace(0, 1, len(vel_plasma_deuteron_profile)), - vel_plasma_deuteron_profile, - color="pink", - linestyle="-", - label=r"$v_{D}$", - ) - axis.plot( - np.linspace(0, 1, len(vel_plasma_triton_profile)), - vel_plasma_triton_profile, - color="green", - linestyle="-", - label=r"$v_{T}$", - ) - axis.plot( - np.linspace(0, 1, len(vel_plasma_alpha_thermal_profile)), - vel_plasma_alpha_thermal_profile, - color="red", - 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]") - axis.set_xlabel("$\\rho \\ [r/a]$") - axis.grid(True, which="both", linestyle="--", alpha=0.5) - axis.set_xlim([0, 1.025]) - axis.minorticks_on() - axis.legend() - - -def plot_electron_frequency_profile(axis, mfile_data, scan): - """Plot the electron thermal frequency profile on the given axis. - - Parameters - ---------- - axis : - - mfile_data : - - scan : - - """ - freq_plasma_electron_profile = [ - mfile_data.data[f"freq_plasma_electron_profile{i}"].get_scan(scan) - for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) - ] - freq_plasma_larmor_toroidal_electron_profile = [ - mfile_data.data[f"freq_plasma_larmor_toroidal_electron_profile{i}"].get_scan( - scan - ) - for i in range( - 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(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, - color="red", - linestyle="-", - 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}$ | Plasma Frequency" - ) - # 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}$ | 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]") - 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", 30)) - - axis.legend() - - -def plot_ion_frequency_profile(axis, mfile_data, scan): - freq_plasma_larmor_toroidal_deuteron_profile = [ - mfile_data.data[f"freq_plasma_larmor_toroidal_deuteron_profile{i}"].get_scan( - scan - ) - for i in range( - 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) - ) - ] - - freq_plasma_larmor_toroidal_triton_profile = [ - mfile_data.data[f"freq_plasma_larmor_toroidal_triton_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_deuteron_profile)), - np.array(freq_plasma_larmor_toroidal_deuteron_profile) / 1e6, - color="red", - linestyle="-", - label=r"$f_{Larmor,toroidal,D}$", - ) - axis.plot( - np.linspace(-1, 1, len(freq_plasma_larmor_toroidal_triton_profile)), - np.array(freq_plasma_larmor_toroidal_triton_profile) / 1e6, - color="green", - linestyle="-", - label=r"$f_{Larmor,toroidal,T}$", - ) - - axis.set_ylabel("Frequency [MHz]") - axis.set_xlabel("$\\rho \\ [r/a]$") - axis.grid(True, which="both", linestyle="--", alpha=0.5) - axis.minorticks_on() - axis.legend() - - -def plot_plasma_coloumb_logarithms(axis, mfile_data, scan): - """Plot the plasma coloumb logarithms on the given axis. - - Parameters - ---------- - axis : - - mfile_data : - - scan : - - """ - plasma_coulomb_log_electron_electron_profile = [ - mfile_data.data[f"plasma_coulomb_log_electron_electron_profile{i}"].get_scan( - scan - ) - for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) - ] - - plasma_coulomb_log_electron_deuteron_profile = [ - mfile_data.data[f"plasma_coulomb_log_electron_deuteron_profile{i}"].get_scan( - scan - ) - for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan))) - ] - - plasma_coulomb_log_electron_triton_profile = [ - mfile_data.data[f"plasma_coulomb_log_electron_triton_profile{i}"].get_scan(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))) - ] - - 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, - color="blue", - linestyle="-", - label=r"$ln \Lambda_{e-e}$", - ) - - axis.plot( - np.linspace(0, 1, len(plasma_coulomb_log_electron_deuteron_profile)), - plasma_coulomb_log_electron_deuteron_profile, - color="pink", - linestyle="-", - label=r"$ln \Lambda_{e-D}$", - ) - - axis.plot( - np.linspace(0, 1, len(plasma_coulomb_log_electron_triton_profile)), - plasma_coulomb_log_electron_triton_profile, - color="green", - linestyle="-", - 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.plot( - np.linspace(0, 1, len(plasma_coulomb_log_electron_alpha_thermal_profile)), - plasma_coulomb_log_electron_alpha_thermal_profile, - color="red", - 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) - axis.minorticks_on() - axis.legend() - - def plot_equality_constraint_equations(axis: plt.Axes, m_file_data: mf.MFile, scan: int): """Plot the equality constraints for a solution and their normalised residuals @@ -13746,19 +13347,37 @@ def main_plot( plot_density_limit_comparison(figs[15].add_subplot(221), m_file, scan) plot_confinement_time_comparison(figs[16].add_subplot(224), m_file, scan) - plot_debye_length_profile(figs[17].add_subplot(232), m_file, scan) - plot_velocity_profile(figs[17].add_subplot(233), m_file, scan) - plot_plasma_coloumb_logarithms(figs[17].add_subplot(231), m_file, scan) + DetailedPhysics.plot_debye_length_profile(figs[17].add_subplot(232), m_file, scan) + DetailedPhysics.plot_velocity_profile(figs[17].add_subplot(233), m_file, scan) + DetailedPhysics.plot_plasma_coloumb_logarithms( + figs[17].add_subplot(231), m_file, scan + ) + DetailedPhysics.plot_collision_time_profile(figs[17].add_subplot(234), m_file, scan) + DetailedPhysics.plot_collision_frequency_profile( + figs[17].add_subplot(235), m_file, scan + ) + DetailedPhysics.plot_mean_free_path_profile(figs[17].add_subplot(236), m_file, scan) + + DetailedPhysics.plot_ion_slowing_down_time_profile( + figs[18].add_subplot(231), m_file, scan + ) + + DetailedPhysics.plot_resistivity_profile(figs[18].add_subplot(232), m_file, scan) + + ax_electron_freq = figs[19].add_subplot(211) + DetailedPhysics.plot_electron_frequency_profile(ax_electron_freq, m_file, scan) - plot_electron_frequency_profile(figs[17].add_subplot(212), m_file, scan) + ax_ion_freq = figs[19].add_subplot(413, sharex=ax_electron_freq) + DetailedPhysics.plot_ion_frequency_profile(ax_ion_freq, m_file, scan) - plot_ion_frequency_profile(figs[18].add_subplot(311), m_file, scan) + ax_larmor = figs[19].add_subplot(414, sharex=ax_electron_freq) + DetailedPhysics.plot_larmor_radius_profile(ax_larmor, m_file, scan) - plot_larmor_radius_profile(figs[18].add_subplot(313), m_file, scan) + figs[19].subplots_adjust(hspace=0.5) # Plot poloidal cross-section poloidal_cross_section( - figs[19].add_subplot(121, aspect="equal"), + figs[20].add_subplot(121, aspect="equal"), m_file, scan, demo_ranges, @@ -13768,7 +13387,7 @@ def main_plot( # Plot toroidal cross-section toroidal_cross_section( - figs[19].add_subplot(122, aspect="equal"), + figs[20].add_subplot(122, aspect="equal"), m_file, scan, demo_ranges, @@ -13776,27 +13395,27 @@ def main_plot( ) # Plot color key - ax17 = figs[19].add_subplot(222) + ax17 = figs[20].add_subplot(222) ax17.set_position([0.5, 0.5, 0.5, 0.5]) color_key(ax17, m_file, scan, colour_scheme) plot_full_machine_poloidal_cross_section( - figs[20].add_subplot(111, aspect="equal"), + figs[21].add_subplot(111, aspect="equal"), m_file, scan, radial_build, colour_scheme, ) - ax18 = figs[21].add_subplot(211) + ax18 = figs[22].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[22].add_subplot(211) + ax185 = figs[23].add_subplot(211) ax185.set_position([0.1, 0.61, 0.8, 0.32]) - ax18b = figs[22].add_subplot(212) + ax18b = figs[23].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) @@ -13804,62 +13423,62 @@ 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[23].add_subplot(221, aspect="equal") + ax19 = figs[24].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[23]) + plot_superconducting_tf_wp(ax19, m_file, scan, figs[24]) # TF coil turn structure - ax20 = figs[24].add_subplot(325, aspect="equal") + ax20 = figs[25].add_subplot(325, aspect="equal") ax20.set_position([0.025, 0.5, 0.4, 0.4]) - plot_tf_cable_in_conduit_turn(ax20, figs[24], m_file, scan) - plot_205 = figs[24].add_subplot(223, aspect="equal") + plot_tf_cable_in_conduit_turn(ax20, figs[25], m_file, scan) + plot_205 = figs[25].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[24], m_file, scan) + plot_cable_in_conduit_cable(plot_205, figs[25], m_file, scan) else: - ax19 = figs[23].add_subplot(211, aspect="equal") + ax19 = figs[24].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[23]) - plot_resistive_tf_info(ax19, m_file, scan, figs[23]) + plot_resistive_tf_wp(ax19, m_file, scan, figs[24]) + plot_resistive_tf_info(ax19, m_file, scan, figs[24]) plot_tf_coil_structure( - figs[25].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme + figs[26].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme ) - plot_plasma_outboard_toroidal_ripple_map(figs[26], m_file, scan) + plot_plasma_outboard_toroidal_ripple_map(figs[27], m_file, scan) - plot_tf_stress(figs[27].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) + plot_tf_stress(figs[28].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) - plot_current_profiles_over_time(figs[28].add_subplot(111), m_file, scan) + plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan) plot_cs_coil_structure( - figs[29].add_subplot(121, aspect="equal"), figs[29], m_file, scan + figs[30].add_subplot(121, aspect="equal"), figs[30], m_file, scan ) plot_cs_turn_structure( - figs[29].add_subplot(224, aspect="equal"), figs[29], m_file, scan + figs[30].add_subplot(224, aspect="equal"), figs[30], m_file, scan ) plot_first_wall_top_down_cross_section( - figs[30].add_subplot(221, aspect="equal"), m_file, scan + figs[31].add_subplot(221, aspect="equal"), m_file, scan ) - plot_first_wall_poloidal_cross_section(figs[30].add_subplot(122), m_file, scan) - plot_fw_90_deg_pipe_bend(figs[30].add_subplot(337), m_file, scan) + plot_first_wall_poloidal_cross_section(figs[31].add_subplot(122), m_file, scan) + plot_fw_90_deg_pipe_bend(figs[31].add_subplot(337), m_file, scan) - plot_blkt_pipe_bends(figs[31], m_file, scan) - ax_blanket = figs[31].add_subplot(122, aspect="equal") - plot_blkt_structure(ax_blanket, figs[31], m_file, scan, radial_build, colour_scheme) + plot_blkt_pipe_bends(figs[32], m_file, scan) + ax_blanket = figs[32].add_subplot(122, aspect="equal") + plot_blkt_structure(ax_blanket, figs[32], m_file, scan, radial_build, colour_scheme) plot_main_power_flow( - figs[32].add_subplot(111, aspect="equal"), m_file, scan, figs[32] + figs[33].add_subplot(111, aspect="equal"), m_file, scan, figs[33] ) - ax24 = figs[33].add_subplot(111) + ax24 = figs[34].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[33]) + plot_system_power_profiles_over_time(ax24, m_file, scan, figs[34]) def create_thickness_builds(m_file, scan: int): @@ -13936,7 +13555,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(34)] + pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(35)] # run main_plot main_plot( diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index fccd7a829..8d818f237 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1148,6 +1148,9 @@ res_plasma: float = None """plasma resistance (ohm)""" +rho_plasma_spitzer_classical_profile: list[float] = None +"""Profile of plasma Spitzer classical resistivity (ohm m)""" + t_plasma_res_diffusion: float = None """plasma current resistive diffusion time (s)""" @@ -1327,6 +1330,9 @@ plasma_coulomb_log_electron_alpha_thermal_profile: list[float] = None """Profile of electron-alpha Coulomb logarithm in plasma""" +t_plasma_electron_alpha_spitzer_slow_profile: list[float] = None +"""Profile of electron-alpha Spitzer slowing down time in plasma (s)""" + freq_plasma_electron_profile: list[float] = None """Electron plasma frequency profile (Hz)""" @@ -1345,6 +1351,42 @@ freq_plasma_upper_hybrid_profile: list[float] = None """Profile of upper hybrid frequency in plasma (Hz)""" +t_plasma_electron_electron_collision_profile: list[float] = None +"""Profile of electron-electron collision time in plasma (s)""" + +t_plasma_electron_deuteron_collision_profile: list[float] = None +"""Profile of electron-deuteron collision time in plasma (s)""" + +t_plasma_electron_triton_collision_profile: list[float] = None +"""Profile of electron-triton collision time in plasma (s)""" + +t_plasma_electron_alpha_thermal_collision_profile: list[float] = None +"""Profile of electron-alpha collision time in plasma (s)""" + +freq_plasma_electron_electron_collision_profile: list[float] = None +"""Profile of electron-electron collision frequency in plasma (Hz)""" + +freq_plasma_electron_deuteron_collision_profile: list[float] = None +"""Profile of electron-deuteron collision frequency in plasma (Hz)""" + +freq_plasma_electron_triton_collision_profile: list[float] = None +"""Profile of electron-triton collision frequency in plasma (Hz)""" + +freq_plasma_electron_alpha_thermal_collision_profile: list[float] = None +"""Profile of electron-alpha collision frequency in plasma (Hz)""" + +len_plasma_electron_electron_mean_free_path_profile: list[float] = None +"""Profile of electron-electron mean free path in plasma (m)""" + +len_plasma_electron_deuteron_mean_free_path_profile: list[float] = None +"""Profile of electron-deuteron mean free path in plasma (m)""" + +len_plasma_electron_triton_mean_free_path_profile: list[float] = None +"""Profile of electron-triton mean free path in plasma (m)""" + +len_plasma_electron_alpha_thermal_mean_free_path_profile: list[float] = None +"""Profile of electron-alpha mean free path in plasma (m)""" + def init_physics_module(): """Initialise the physics module""" @@ -1639,6 +1681,7 @@ def init_physics_variables(): f_nd_plasma_oxygen_electron, \ f_res_plasma_neo, \ res_plasma, \ + rho_plasma_spitzer_classical_profile, \ t_plasma_res_diffusion, \ a_plasma_surface, \ a_plasma_surface_outboard, \ @@ -1688,12 +1731,25 @@ def init_physics_variables(): plasma_coulomb_log_electron_triton_profile, \ plasma_coulomb_log_deuteron_triton_profile, \ plasma_coulomb_log_electron_alpha_thermal_profile, \ + t_plasma_electron_alpha_spitzer_slow_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, \ - freq_plasma_upper_hybrid_profile + freq_plasma_upper_hybrid_profile, \ + t_plasma_electron_electron_collision_profile, \ + t_plasma_electron_deuteron_collision_profile, \ + t_plasma_electron_triton_collision_profile, \ + t_plasma_electron_alpha_thermal_collision_profile, \ + freq_plasma_electron_electron_collision_profile, \ + freq_plasma_electron_deuteron_collision_profile, \ + freq_plasma_electron_triton_collision_profile, \ + freq_plasma_electron_alpha_thermal_collision_profile, \ + len_plasma_electron_electron_mean_free_path_profile, \ + len_plasma_electron_deuteron_mean_free_path_profile, \ + len_plasma_electron_triton_mean_free_path_profile, \ + len_plasma_electron_alpha_thermal_mean_free_path_profile m_beam_amu = 0.0 m_fuel_amu = 0.0 @@ -1930,6 +1986,7 @@ def init_physics_variables(): f_nd_plasma_oxygen_electron = 0.0 f_res_plasma_neo = 0.0 res_plasma = 0.0 + rho_plasma_spitzer_classical_profile = [] t_plasma_res_diffusion = 0.0 a_plasma_surface = 0.0 a_plasma_surface_outboard = 0.0 @@ -1984,3 +2041,16 @@ def init_physics_variables(): freq_plasma_larmor_toroidal_deuteron_profile = [] freq_plasma_larmor_toroidal_triton_profile = [] freq_plasma_upper_hybrid_profile = [] + t_plasma_electron_electron_collision_profile = [] + t_plasma_electron_deuteron_collision_profile = [] + t_plasma_electron_triton_collision_profile = [] + t_plasma_electron_alpha_thermal_collision_profile = [] + t_plasma_electron_alpha_spitzer_slow_profile = [] + freq_plasma_electron_electron_collision_profile = [] + freq_plasma_electron_deuteron_collision_profile = [] + freq_plasma_electron_triton_collision_profile = [] + freq_plasma_electron_alpha_thermal_collision_profile = [] + len_plasma_electron_electron_mean_free_path_profile = [] + len_plasma_electron_deuteron_mean_free_path_profile = [] + len_plasma_electron_triton_mean_free_path_profile = [] + len_plasma_electron_alpha_thermal_mean_free_path_profile = [] diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index d313c02ec..d91f29f2f 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -2,9 +2,11 @@ import math from enum import IntEnum +import matplotlib.pyplot as plt import numba as nb import numpy as np +import process.core.io.mfile as mf import process.models.physics.fusion_reactions as reactions import process.models.physics.impurity_radiation as impurity_radiation import process.models.physics.radiation_power as physics_funcs @@ -5055,6 +5057,120 @@ def run(self): for i in range(len(physics_variables.len_plasma_debye_electron_profile)) ]) + # ============================ + # Collision times + # ============================ + + physics_variables.t_plasma_electron_electron_collision_profile = self.calculate_electron_electron_collision_time( + temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y, + nd_plasma_electrons=self.plasma_profile.neprofile.profile_y, + plasma_coulomb_log_electron_electron=physics_variables.plasma_coulomb_log_electron_electron_profile, + ) + + physics_variables.t_plasma_electron_deuteron_collision_profile = self.calculate_electron_ion_collision_time( + temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y, + nd_plasma_ions=( + self.plasma_profile.neprofile.profile_y + * physics_variables.f_plasma_fuel_deuterium + * ( + physics_variables.nd_plasma_fuel_ions_vol_avg + / physics_variables.nd_plasma_electrons_vol_avg + ) + ), + plasma_coulomb_log_electron_ion=physics_variables.plasma_coulomb_log_electron_deuteron_profile, + n_charge_ion=1, + ) + + physics_variables.t_plasma_electron_triton_collision_profile = self.calculate_electron_ion_collision_time( + temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y, + nd_plasma_ions=( + self.plasma_profile.neprofile.profile_y + * physics_variables.f_plasma_fuel_tritium + * ( + physics_variables.nd_plasma_fuel_ions_vol_avg + / physics_variables.nd_plasma_electrons_vol_avg + ) + ), + plasma_coulomb_log_electron_ion=physics_variables.plasma_coulomb_log_electron_triton_profile, + n_charge_ion=1, + ) + + physics_variables.t_plasma_electron_alpha_thermal_collision_profile = self.calculate_electron_ion_collision_time( + temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y, + nd_plasma_ions=( + self.plasma_profile.neprofile.profile_y + * ( + physics_variables.nd_plasma_alphas_vol_avg + / physics_variables.nd_plasma_electrons_vol_avg + ) + ), + plasma_coulomb_log_electron_ion=physics_variables.plasma_coulomb_log_electron_alpha_thermal_profile, + n_charge_ion=2, + ) + + # ============================ + # Collision frequencies + # ============================ + + physics_variables.freq_plasma_electron_electron_collision_profile = ( + 1 / physics_variables.t_plasma_electron_electron_collision_profile + ) + + physics_variables.freq_plasma_electron_deuteron_collision_profile = ( + 1 / physics_variables.t_plasma_electron_deuteron_collision_profile + ) + + physics_variables.freq_plasma_electron_triton_collision_profile = ( + 1 / physics_variables.t_plasma_electron_triton_collision_profile + ) + + physics_variables.freq_plasma_electron_alpha_thermal_collision_profile = ( + 1 / physics_variables.t_plasma_electron_alpha_thermal_collision_profile + ) + + # ============================ + # Mean free paths + # ============================ + + physics_variables.len_plasma_electron_electron_mean_free_path_profile = ( + physics_variables.t_plasma_electron_electron_collision_profile + * physics_variables.vel_plasma_electron_profile + ) + + physics_variables.len_plasma_electron_deuteron_mean_free_path_profile = ( + physics_variables.t_plasma_electron_deuteron_collision_profile + * physics_variables.vel_plasma_electron_profile + ) + + physics_variables.len_plasma_electron_triton_mean_free_path_profile = ( + physics_variables.t_plasma_electron_triton_collision_profile + * physics_variables.vel_plasma_electron_profile + ) + + physics_variables.len_plasma_electron_alpha_thermal_mean_free_path_profile = ( + physics_variables.t_plasma_electron_alpha_thermal_collision_profile + * physics_variables.vel_plasma_electron_profile + ) + + # ============================ + # Spitzer slow down time + # ============================ + + physics_variables.t_plasma_electron_alpha_spitzer_slow_profile = self.calculate_spitzer_ion_slowing_down_time( + m_ion=constants.ALPHA_MASS, + plasma_coulomb_log_electron_ion=physics_variables.plasma_coulomb_log_electron_alpha_thermal_profile, + n_charge_ion=2, + ) + + # ============================ + # Resistivites + # ============================ + + physics_variables.res_plasma_spitzer_profile = self.calculate_spitzer_resistivity( + n_charge=physics_variables.n_charge_plasma_effective_profile, + electron_ion_coulomb_log=physics_variables.plasma_coulomb_log_electron_deuteron_profile, + ) + @staticmethod @nb.njit(cache=True) def calculate_debye_length( @@ -5339,6 +5455,154 @@ def calculate_average_relative_velocity( """ return np.sqrt(velocity_1**2 + velocity_2**2) + @staticmethod + @nb.njit(cache=True) + def calculate_electron_electron_collision_time( + temp_plasma_electron_kev: float | np.ndarray, + nd_plasma_electrons: float | np.ndarray, + plasma_coulomb_log_electron_electron: float | np.ndarray, + ) -> float | np.ndarray: + """Calculate the electron-electron collision time for a plasma (τₑₑ). + + Parameters + ---------- + temp_plasma_electron_kev : float | np.ndarray + Electron temperature in keV. + nd_plasma_electrons : float | np.ndarray + Electron density (m^-3). + plasma_coulomb_log_electron_electron : float | np.ndarray + Coulomb logarithm for electron-electron collisions. + + Returns + ------- + float | np.ndarray + Electron-electron collision time (s). + """ + return ( + 12 + * np.sqrt(2) + * np.pi**1.5 + * constants.EPSILON0**2 + * constants.ELECTRON_MASS**0.5 + * (temp_plasma_electron_kev * constants.KILOELECTRON_VOLT) ** 1.5 + ) / ( + plasma_coulomb_log_electron_electron + * constants.ELECTRON_CHARGE**4 + * nd_plasma_electrons + ) + + @staticmethod + @nb.njit(cache=True) + def calculate_electron_ion_collision_time( + temp_plasma_electron_kev: float | np.ndarray, + nd_plasma_ions: float | np.ndarray, + plasma_coulomb_log_electron_ion: float | np.ndarray, + n_charge_ion: float = 1.0, + ) -> float | np.ndarray: + """Calculate the electron-ion collision time for a plasma (τₑᵢ). + + Parameters + ---------- + temp_plasma_electron_kev : float | np.ndarray + Electron temperature in keV. + nd_plasma_ions : float | np.ndarray + Ion density (m^-3). + plasma_coulomb_log_electron_ion : float | np.ndarray + Coulomb logarithm for electron-ion collisions. + n_charge_ion : float + Charge number (Z) of the ion. + + Returns + ------- + float | np.ndarray + Electron-ion collision time (s). + """ + return ( + 12 + * np.pi**1.5 + * constants.EPSILON0**2 + * constants.ELECTRON_MASS**0.5 + * (temp_plasma_electron_kev * constants.KILOELECTRON_VOLT) ** 1.5 + ) / ( + np.sqrt(2) + * nd_plasma_ions + * n_charge_ion**2 + * plasma_coulomb_log_electron_ion + * constants.ELECTRON_CHARGE**4 + ) + + def calculate_spitzer_ion_slowing_down_time( + self, + m_ion: float, + plasma_coulomb_log_electron_ion: float | np.ndarray, + n_charge_ion: float = 1.0, + ) -> float: + """ + Calculate the Spitzer slowing down time for ions in a plasma. + + Parameters + ---------- + m_ion : float + Mass of the ion (kg). + plasma_coulomb_log_electron_ion : float | np.ndarray + Coulomb logarithm for electron-ion collisions. + n_charge_ion : float + Charge number (Z) of the ion. + + Returns + ------- + float + Spitzer slowing down time (s). + """ + + return ( + (3 * (2 * np.pi) ** 1.5 * constants.EPSILON0**2) + * ( + m_ion + * (self.plasma_profile.teprofile.profile_y * constants.KILOELECTRON_VOLT) + ** 1.5 + ) + / ( + self.plasma_profile.neprofile.profile_y + * constants.ELECTRON_CHARGE**4 + * plasma_coulomb_log_electron_ion + * n_charge_ion**2 + * np.sqrt(constants.ELECTRON_MASS) + ) + ) + + def calculate_spitzer_resistivity( + self, n_charge: float | np.ndarray, electron_ion_coulomb_log: float | np.ndarray + ) -> float | np.ndarray: + """ + Calculate the classical Spitzer resistivity for a plasma. + + Parameters + ---------- + n_charge : float | np.ndarray + Charge number (Z) of the ion. + electron_ion_coulomb_log : float | np.ndarray + Coulomb logarithm for electron-ion collisions. + + Returns + ------- + float | np.ndarray + Spitzer resistivity (Ohm m). + """ + + return ( + (4 * np.sqrt(2 * np.pi) / 3) + * (1 / (4 * np.pi * constants.EPSILON0) ** 2) + * ( + n_charge + * constants.ELECTRON_CHARGE**2 + * constants.ELECTRON_MASS**0.5 + * electron_ion_coulomb_log + ) + / (self.plasma_profile.teprofile.profile_y * constants.KILOELECTRON_VOLT) + ** 1.5 + ) + def output_detailed_physics(self): """Outputs detailed physics variables to file.""" @@ -5521,3 +5785,954 @@ def output_detailed_physics(self): f"(plasma_coulomb_log_electron_alpha_thermal_profile{i})", physics_variables.plasma_coulomb_log_electron_alpha_thermal_profile[i], ) + + po.osubhd(self.outfile, "Collision Times:") + + for i in range( + len(physics_variables.t_plasma_electron_electron_collision_profile) + ): + po.ovarre( + self.mfile, + f"Electron-electron collision time at point {i}", + f"(t_plasma_electron_electron_collision_profile{i})", + physics_variables.t_plasma_electron_electron_collision_profile[i], + ) + + for i in range( + len(physics_variables.t_plasma_electron_deuteron_collision_profile) + ): + po.ovarre( + self.mfile, + f"Electron-deuteron collision time at point {i}", + f"(t_plasma_electron_deuteron_collision_profile{i})", + physics_variables.t_plasma_electron_deuteron_collision_profile[i], + ) + + for i in range( + len(physics_variables.t_plasma_electron_triton_collision_profile) + ): + po.ovarre( + self.mfile, + f"Electron-triton collision time at point {i}", + f"(t_plasma_electron_triton_collision_profile{i})", + physics_variables.t_plasma_electron_triton_collision_profile[i], + ) + + for i in range( + len(physics_variables.t_plasma_electron_alpha_thermal_collision_profile) + ): + po.ovarre( + self.mfile, + f"Electron-alpha thermal collision time at point {i}", + f"(t_plasma_electron_alpha_thermal_collision_profile{i})", + physics_variables.t_plasma_electron_alpha_thermal_collision_profile[i], + ) + + po.osubhd(self.outfile, "Collision Frequencies:") + + for i in range( + len(physics_variables.freq_plasma_electron_electron_collision_profile) + ): + po.ovarre( + self.mfile, + f"Electron-electron collision frequency at point {i}", + f"(freq_plasma_electron_electron_collision_profile{i})", + physics_variables.freq_plasma_electron_electron_collision_profile[i], + ) + for i in range( + len(physics_variables.freq_plasma_electron_deuteron_collision_profile) + ): + po.ovarre( + self.mfile, + f"Electron-deuteron collision frequency at point {i}", + f"(freq_plasma_electron_deuteron_collision_profile{i})", + physics_variables.freq_plasma_electron_deuteron_collision_profile[i], + ) + for i in range( + len(physics_variables.freq_plasma_electron_triton_collision_profile) + ): + po.ovarre( + self.mfile, + f"Electron-triton collision frequency at point {i}", + f"(freq_plasma_electron_triton_collision_profile{i})", + physics_variables.freq_plasma_electron_triton_collision_profile[i], + ) + for i in range( + len(physics_variables.freq_plasma_electron_alpha_thermal_collision_profile) + ): + po.ovarre( + self.mfile, + f"Electron-alpha thermal collision frequency at point {i}", + f"(freq_plasma_electron_alpha_thermal_collision_profile{i})", + physics_variables.freq_plasma_electron_alpha_thermal_collision_profile[ + i + ], + ) + + po.osubhd(self.outfile, "Mean Free Paths:") + + for i in range( + len(physics_variables.len_plasma_electron_electron_mean_free_path_profile) + ): + po.ovarre( + self.mfile, + f"Electron-electron mean free path at point {i}", + f"(len_plasma_electron_electron_mean_free_path_profile{i})", + physics_variables.len_plasma_electron_electron_mean_free_path_profile[i], + ) + for i in range( + len(physics_variables.len_plasma_electron_deuteron_mean_free_path_profile) + ): + po.ovarre( + self.mfile, + f"Electron-deuteron mean free path at point {i}", + f"(len_plasma_electron_deuteron_mean_free_path_profile{i})", + physics_variables.len_plasma_electron_deuteron_mean_free_path_profile[i], + ) + for i in range( + len(physics_variables.len_plasma_electron_triton_mean_free_path_profile) + ): + po.ovarre( + self.mfile, + f"Electron-triton mean free path at point {i}", + f"(len_plasma_electron_triton_mean_free_path_profile{i})", + physics_variables.len_plasma_electron_triton_mean_free_path_profile[i], + ) + for i in range( + len( + physics_variables.len_plasma_electron_alpha_thermal_mean_free_path_profile + ) + ): + po.ovarre( + self.mfile, + f"Electron-alpha thermal mean free path at point {i}", + f"(len_plasma_electron_alpha_thermal_mean_free_path_profile{i})", + physics_variables.len_plasma_electron_alpha_thermal_mean_free_path_profile[ + i + ], + ) + + po.osubhd(self.outfile, "Spitzer slowing down times:") + + for i in range( + len(physics_variables.t_plasma_electron_alpha_spitzer_slow_profile) + ): + po.ovarre( + self.mfile, + f"Electron-alpha Spitzer slowing down time at point {i}", + f"(t_plasma_electron_alpha_spitzer_slow_profile{i})", + physics_variables.t_plasma_electron_alpha_spitzer_slow_profile[i], + ) + + po.osubhd(self.outfile, "Resistivities:") + + for i in range(len(physics_variables.res_plasma_spitzer_profile)): + po.ovarre( + self.mfile, + f"Plasma Spitzer resistivity at point {i}", + f"(res_plasma_spitzer_profile{i})", + physics_variables.res_plasma_spitzer_profile[i], + ) + + @staticmethod + 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_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, + color="red", + linestyle="-", + 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" 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() + axis.legend() + + @staticmethod + def plot_debye_length_profile(axis: plt.Axes, mfile_data: mf.MFile, scan: int): + """Plot the Debye length profile on the given axis. + + Parameters + ---------- + axis : plt.Axes + Axis object to plot on + mfile_data : mf.MFile + MFILE data object + scan : int + Scan number to use + + """ + len_plasma_debye_electron_profile = [ + mfile_data.data[f"len_plasma_debye_electron_profile{i}"].get_scan(scan) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + # Convert to micrometres (1e-6 m) + len_plasma_debye_electron_profile_um = [ + length * 1e6 for length in len_plasma_debye_electron_profile + ] + + axis.plot( + np.linspace(0, 1, len(len_plasma_debye_electron_profile_um)), + len_plasma_debye_electron_profile_um, + color="blue", + linestyle="-", + label=r"$\lambda_{Debye,e}$", + ) + + axis.set_ylabel(r"Debye Length [$\mu$m]") + + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.set_xlim([0, 1.025]) + axis.minorticks_on() + axis.legend() + + @staticmethod + def plot_velocity_profile(axis: plt.Axes, mfile_data: mf.MFile, scan: int) -> None: + """Plot the electron thermal velocity profile on the given axis. + + Parameters + ---------- + axis : plt.Axes + Axis object to plot on + mfile_data : mf.MFile + MFILE data object + scan : int + Scan number to use + """ + + vel_plasma_electron_profile = [ + mfile_data.data[f"vel_plasma_electron_profile{i}"].get_scan(scan) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + vel_plasma_deuteron_profile = [ + mfile_data.data[f"vel_plasma_deuteron_profile{i}"].get_scan(scan) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + vel_plasma_triton_profile = [ + 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)) + ) + ] + + 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, + color="blue", + linestyle="-", + label=r"$v_{e}$", + ) + axis.plot( + np.linspace(0, 1, len(vel_plasma_deuteron_profile)), + vel_plasma_deuteron_profile, + color="pink", + linestyle="-", + label=r"$v_{D}$", + ) + axis.plot( + np.linspace(0, 1, len(vel_plasma_triton_profile)), + vel_plasma_triton_profile, + color="green", + linestyle="-", + label=r"$v_{T}$", + ) + axis.plot( + np.linspace(0, 1, len(vel_plasma_alpha_thermal_profile)), + vel_plasma_alpha_thermal_profile, + color="red", + 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]") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.set_xlim([0, 1.025]) + axis.minorticks_on() + axis.legend() + + @staticmethod + def plot_electron_frequency_profile( + axis: plt.Axes, mfile_data: mf.MFile, scan: int + ) -> None: + """Plot the electron thermal frequency profile on the given axis. + + Parameters + ---------- + axis : plt.Axes + Axis object to plot on + mfile_data : mf.MFile + MFILE data object + scan : int + Scan number to use + """ + + freq_plasma_electron_profile = [ + mfile_data.data[f"freq_plasma_electron_profile{i}"].get_scan(scan) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + freq_plasma_larmor_toroidal_electron_profile = [ + mfile_data.data[f"freq_plasma_larmor_toroidal_electron_profile{i}"].get_scan( + scan + ) + for i in range( + 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(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, + color="red", + linestyle="-", + 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}$ | Plasma Frequency" + ) + # 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}$ | 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]") + 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", 30)) + + axis.legend() + + @staticmethod + def plot_ion_frequency_profile( + axis: plt.Axes, mfile_data: mf.MFile, scan: int + ) -> None: + """Plot the ion thermal frequency profile on the given axis. + + Parameters + ---------- + axis : plt.Axes + Axis object to plot on + mfile_data : mf.MFile + MFILE data object + scan : int + Scan number to use + + """ + + freq_plasma_larmor_toroidal_deuteron_profile = [ + mfile_data.data[f"freq_plasma_larmor_toroidal_deuteron_profile{i}"].get_scan( + scan + ) + for i in range( + 2 * int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + freq_plasma_larmor_toroidal_triton_profile = [ + mfile_data.data[f"freq_plasma_larmor_toroidal_triton_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_deuteron_profile)), + np.array(freq_plasma_larmor_toroidal_deuteron_profile) / 1e6, + color="red", + linestyle="-", + label=r"$f_{Larmor,toroidal,D}$", + ) + axis.plot( + np.linspace(-1, 1, len(freq_plasma_larmor_toroidal_triton_profile)), + np.array(freq_plasma_larmor_toroidal_triton_profile) / 1e6, + color="green", + linestyle="-", + label=r"$f_{Larmor,toroidal,T}$", + ) + + axis.set_ylabel("Frequency [MHz]") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + + @staticmethod + def plot_plasma_coloumb_logarithms( + axis: plt.Axes, mfile_data: mf.MFile, scan: int + ) -> None: + """Plot the plasma coloumb logarithms on the given axis. + + Parameters + ---------- + axis : plt.Axes + Axis object to plot on + mfile_data : mf.MFile + MFILE data object + scan : int + Scan number to use + + """ + plasma_coulomb_log_electron_electron_profile = [ + mfile_data.data[f"plasma_coulomb_log_electron_electron_profile{i}"].get_scan( + scan + ) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + plasma_coulomb_log_electron_deuteron_profile = [ + mfile_data.data[f"plasma_coulomb_log_electron_deuteron_profile{i}"].get_scan( + scan + ) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + plasma_coulomb_log_electron_triton_profile = [ + mfile_data.data[f"plasma_coulomb_log_electron_triton_profile{i}"].get_scan( + 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)) + ) + ] + + 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, + color="blue", + linestyle="-", + label=r"$ln \Lambda_{e-e}$", + ) + + axis.plot( + np.linspace(0, 1, len(plasma_coulomb_log_electron_deuteron_profile)), + plasma_coulomb_log_electron_deuteron_profile, + color="pink", + linestyle="-", + label=r"$ln \Lambda_{e-D}$", + ) + + axis.plot( + np.linspace(0, 1, len(plasma_coulomb_log_electron_triton_profile)), + plasma_coulomb_log_electron_triton_profile, + color="green", + linestyle="-", + 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.plot( + np.linspace(0, 1, len(plasma_coulomb_log_electron_alpha_thermal_profile)), + plasma_coulomb_log_electron_alpha_thermal_profile, + color="red", + 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) + axis.minorticks_on() + axis.legend() + + @staticmethod + def plot_collision_time_profile( + axis: plt.Axes, mfile_data: mf.MFile, scan: int + ) -> None: + """Plot the plasma collision times on the given axis. + + Parameters + ---------- + axis : plt.Axes + Axis object to plot on + mfile_data : mf.MFile + MFILE data object + scan : int + Scan number to use + + """ + t_plasma_electron_electron_collision_profile = [ + mfile_data.data[f"t_plasma_electron_electron_collision_profile{i}"].get_scan( + scan + ) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + t_plasma_electron_deuteron_collision_profile = [ + mfile_data.data[f"t_plasma_electron_deuteron_collision_profile{i}"].get_scan( + scan + ) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + t_plasma_electron_triton_collision_profile = [ + mfile_data.data[f"t_plasma_electron_triton_collision_profile{i}"].get_scan( + scan + ) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + t_plasma_electron_alpha_thermal_collision_profile = [ + mfile_data.data[ + f"t_plasma_electron_alpha_thermal_collision_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(t_plasma_electron_electron_collision_profile)), + t_plasma_electron_electron_collision_profile, + color="blue", + linestyle="-", + label=r"$\tau_{e-e}$", + ) + + axis.plot( + np.linspace(0, 1, len(t_plasma_electron_deuteron_collision_profile)), + t_plasma_electron_deuteron_collision_profile, + color="pink", + linestyle="-", + label=r"$\tau_{e-D}$", + ) + + axis.plot( + np.linspace(0, 1, len(t_plasma_electron_triton_collision_profile)), + t_plasma_electron_triton_collision_profile, + color="green", + linestyle="-", + label=r"$\tau_{e-T}$", + ) + + axis.plot( + np.linspace(0, 1, len(t_plasma_electron_alpha_thermal_collision_profile)), + t_plasma_electron_alpha_thermal_collision_profile, + color="red", + linestyle="-", + label=r"$\tau_{e-\alpha,thermal}$", + ) + + axis.set_yscale("log") + axis.set_ylabel("Collision Time [s]") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + + @staticmethod + def plot_collision_frequency_profile( + axis: plt.Axes, mfile_data: mf.MFile, scan: int + ) -> None: + """Plot the plasma collision frequencies on the given axis. + + Parameters + ---------- + axis : plt.Axes + Axis object to plot on + mfile_data : mf.MFile + MFILE data object + scan : int + Scan number to use + + """ + freq_plasma_electron_electron_collision_profile = [ + mfile_data.data[ + f"freq_plasma_electron_electron_collision_profile{i}" + ].get_scan(scan) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + freq_plasma_electron_deuteron_collision_profile = [ + mfile_data.data[ + f"freq_plasma_electron_deuteron_collision_profile{i}" + ].get_scan(scan) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + freq_plasma_electron_triton_collision_profile = [ + mfile_data.data[ + f"freq_plasma_electron_triton_collision_profile{i}" + ].get_scan(scan) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + freq_plasma_electron_alpha_thermal_collision_profile = [ + mfile_data.data[ + f"freq_plasma_electron_alpha_thermal_collision_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(freq_plasma_electron_electron_collision_profile)), + freq_plasma_electron_electron_collision_profile, + color="blue", + linestyle="-", + label=r"$\nu_{e-e}$", + ) + + axis.plot( + np.linspace(0, 1, len(freq_plasma_electron_deuteron_collision_profile)), + freq_plasma_electron_deuteron_collision_profile, + color="pink", + linestyle="-", + label=r"$\nu_{e-D}$", + ) + + axis.plot( + np.linspace(0, 1, len(freq_plasma_electron_triton_collision_profile)), + freq_plasma_electron_triton_collision_profile, + color="green", + linestyle="-", + label=r"$\nu_{e-T}$", + ) + + axis.plot( + np.linspace(0, 1, len(freq_plasma_electron_alpha_thermal_collision_profile)), + freq_plasma_electron_alpha_thermal_collision_profile, + color="red", + linestyle="-", + label=r"$\nu_{e-\alpha,thermal}$", + ) + axis.set_yscale("log") + axis.set_ylabel("Collision Frequency [Hz]") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + + @staticmethod + def plot_mean_free_path_profile( + axis: plt.Axes, mfile_data: mf.MFile, scan: int + ) -> None: + """Plot the plasma mean free path on the given axis. + + Parameters + ---------- + axis : plt.Axes + Axis object to plot on + mfile_data : mf.MFile + MFILE data object + scan : int + Scan number to use + + """ + len_plasma_electron_electron_mean_free_path_profile = [ + mfile_data.data[ + f"len_plasma_electron_electron_mean_free_path_profile{i}" + ].get_scan(scan) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + len_plasma_electron_deuteron_mean_free_path_profile = [ + mfile_data.data[ + f"len_plasma_electron_deuteron_mean_free_path_profile{i}" + ].get_scan(scan) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + len_plasma_electron_triton_mean_free_path_profile = [ + mfile_data.data[ + f"len_plasma_electron_triton_mean_free_path_profile{i}" + ].get_scan(scan) + for i in range( + int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)) + ) + ] + + len_plasma_electron_alpha_thermal_mean_free_path_profile = [ + mfile_data.data[ + f"len_plasma_electron_alpha_thermal_mean_free_path_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(len_plasma_electron_electron_mean_free_path_profile)), + len_plasma_electron_electron_mean_free_path_profile, + color="blue", + linestyle="-", + label=r"$\lambda_{mfp,e-e}$", + ) + + axis.plot( + np.linspace(0, 1, len(len_plasma_electron_deuteron_mean_free_path_profile)), + len_plasma_electron_deuteron_mean_free_path_profile, + color="pink", + linestyle="-", + label=r"$\lambda_{mfp,e-D}$", + ) + + axis.plot( + np.linspace(0, 1, len(len_plasma_electron_triton_mean_free_path_profile)), + len_plasma_electron_triton_mean_free_path_profile, + color="green", + linestyle="-", + label=r"$\lambda_{mfp,e-T}$", + ) + axis.plot( + np.linspace( + 0, 1, len(len_plasma_electron_alpha_thermal_mean_free_path_profile) + ), + len_plasma_electron_alpha_thermal_mean_free_path_profile, + color="red", + linestyle="-", + label=r"$\lambda_{mfp,e-\alpha,thermal}$", + ) + axis.set_yscale("log") + axis.set_ylabel("Mean Free Path [m]") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + + @staticmethod + def plot_ion_slowing_down_time_profile( + axis: plt.Axes, mfile_data: mf.MFile, scan: int + ) -> None: + """Plot the plasma Spitzer slowing down time on the given axis. + + Parameters + ---------- + axis : plt.Axes + Axis object to plot on + mfile_data : mf.MFile + MFILE data object + scan : int + Scan number to use + + """ + t_plasma_electron_alpha_spitzer_slow_profile = [ + mfile_data.data[f"t_plasma_electron_alpha_spitzer_slow_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(t_plasma_electron_alpha_spitzer_slow_profile)), + t_plasma_electron_alpha_spitzer_slow_profile, + color="red", + linestyle="-", + label=r"$\tau_{e-\alpha,Spitzer}$", + ) + + axis.set_yscale("log") + axis.set_ylabel("Spitzer Slowing Down Time [s]") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + + @staticmethod + def plot_resistivity_profile( + axis: plt.Axes, mfile_data: mf.MFile, scan: int + ) -> None: + """Plot the plasma resistivity on the given axis. + + Parameters + ---------- + axis : plt.Axes + Axis object to plot on + mfile_data : mf.MFile + MFILE data object + scan : int + Scan number to use + + """ + res_plasma_spitzer_profile = [ + mfile_data.data[f"res_plasma_spitzer_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(res_plasma_spitzer_profile)), + res_plasma_spitzer_profile, + color="red", + linestyle="-", + label=r"$\eta_{Spitzer}$", + ) + + axis.set_yscale("log") + axis.set_ylabel("Resistivity [Ohm m]") + axis.set_xlabel("$\\rho \\ [r/a]$") + axis.grid(True, which="both", linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend()