From 1db3795dc9f76785b5ea557d219732707ff13055 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 1 Dec 2025 11:31:48 +0000 Subject: [PATCH 1/4] Add function to plot radiation density contours in radial profile --- process/io/plot_proc.py | 149 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 149 insertions(+) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 22612945c9..3a561b7e7a 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -4384,6 +4384,154 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: # --- +def plot_rad_contour(axis, mfile_data, scan, impp): + # read in the impurity data + imp_data = read_imprad_data(2, impp) + + # find impurity densities + imp_frac = np.array([ + mfile_data.data["f_nd_impurity_electrons(01)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(02)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(03)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(04)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(05)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(06)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(07)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(08)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(09)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(10)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(11)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(12)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(13)"].get_scan(scan), + mfile_data.data["f_nd_impurity_electrons(14)"].get_scan(scan), + ]) + + n_plasma_profile_elements = int( + mfile_data.data["n_plasma_profile_elements"].get_scan(scan) + ) + alphan = mfile_data.data["alphan"].get_scan(scan) + alphat = mfile_data.data["alphat"].get_scan(scan) + nd_plasma_electron_on_axis = mfile_data.data["nd_plasma_electron_on_axis"].get_scan( + scan + ) + temp_plasma_electron_on_axis_kev = mfile_data.data[ + "temp_plasma_electron_on_axis_kev" + ].get_scan(scan) + radius_plasma_pedestal_density_norm = mfile_data.data[ + "radius_plasma_pedestal_density_norm" + ].get_scan(scan) + radius_plasma_pedestal_temp_norm = mfile_data.data[ + "radius_plasma_pedestal_temp_norm" + ].get_scan(scan) + + if i_plasma_pedestal == 0: + # Intialise the radius + rho = np.linspace(0, 1.0, n_plasma_profile_elements) + + # The density profile + ne = nd_plasma_electron_on_axis * (1 - rho**2) ** alphan + + # The temperature profile + te = temp_plasma_electron_on_axis_kev * (1 - rho**2) ** alphat + + if i_plasma_pedestal == 1: + rho = np.linspace(0, 1, n_plasma_profile_elements) + + # The density and temperature profile + ne = np.zeros_like(rho) + te = np.zeros_like(rho) + for q in range(rho.shape[0]): + if rho[q] <= radius_plasma_pedestal_density_norm: + ne[q] = ( + nd_plasma_pedestal_electron + + (ne0 - nd_plasma_pedestal_electron) + * (1 - rho[q] ** 2 / radius_plasma_pedestal_density_norm**2) + ** alphan + ) + else: + ne[q] = nd_plasma_separatrix_electron + ( + nd_plasma_pedestal_electron - nd_plasma_separatrix_electron + ) * (1 - rho[q]) / ( + 1 - min(0.9999, radius_plasma_pedestal_density_norm) + ) + + if rho[q] <= radius_plasma_pedestal_temp_norm: + te[q] = ( + temp_plasma_pedestal_kev + + (te0 - temp_plasma_pedestal_kev) + * (1 - (rho[q] / radius_plasma_pedestal_temp_norm) ** tbeta) + ** alphat + ) + else: + te[q] = temp_plasma_separatrix_kev + ( + temp_plasma_pedestal_kev - temp_plasma_separatrix_kev + ) * (1 - rho[q]) / (1 - min(0.9999, radius_plasma_pedestal_temp_norm)) + + # Intailise the radiation profile arrays + pimpden = np.zeros([imp_data.shape[0], te.shape[0]]) + lz = np.zeros([imp_data.shape[0], te.shape[0]]) + prad = np.zeros(te.shape[0]) + + # Intailise the impurity radiation profile + for k in range(te.shape[0]): + for i in range(imp_data.shape[0]): + if te[k] <= imp_data[i][0][0]: + lz[i][k] = imp_data[i][0][1] + elif te[k] >= imp_data[i][imp_data.shape[1] - 1][0]: + lz[i][k] = imp_data[i][imp_data.shape[1] - 1][1] + else: + # Use np.interp for log-log interpolation + log_te_data = np.log([row[0] for row in imp_data[i]]) + log_lz_data = np.log([row[1] for row in imp_data[i]]) + lz[i][k] = np.exp(np.interp(np.log(te[k]), log_te_data, log_lz_data)) + pimpden[i][k] = imp_frac[i] * ne[k] * ne[k] * lz[i][k] + + for l in range(imp_data.shape[0]): # noqa: E741 + prad[k] = prad[k] + pimpden[l][k] * 1.0e-6 + + p_rad_grid, r_grid, z_grid = interp1d_profile(prad, mfile_data, scan) + + # Plot the upper half contour + p_rad_upper = axis.contourf( + r_grid, z_grid, p_rad_grid, levels=50, cmap="plasma", zorder=2 + ) + # Plot the lower half contour (mirror) + axis.contourf(r_grid, -z_grid, p_rad_grid, levels=50, cmap="plasma", zorder=2) + + axis.figure.colorbar( + p_rad_upper, + ax=axis, + label=r"$P_{\mathrm{rad}}$ $[\mathrm{MW.m}^{-3}]$", + location="left", + anchor=(-0.25, 0.5), + ) + + axis.set_xlabel("R [m]") + axis.set_xlim(rmajor - 1.2 * rminor, rmajor + 1.2 * rminor) + axis.set_ylim( + -1.2 * rminor * mfile_data.data["kappa"].get_scan(scan), + 1.2 * mfile_data.data["kappa"].get_scan(scan) * rminor, + ) + axis.set_ylabel("Z [m]") + axis.set_title("Line & Bremsstrahlung Radiation Density Contours") + axis.plot( + rmajor, + 0, + marker="o", + color="red", + markersize=6, + markeredgecolor="black", + zorder=100, + ) + # enable minor ticks and grid for clearer reading + axis.minorticks_on() + axis.grid(True, which="major", linestyle="--", linewidth=0.8, alpha=0.7, zorder=1) + + axis.grid(True, which="minor", linestyle=":", linewidth=0.4, alpha=0.5, zorder=1) + # make minor ticks visible on all sides and draw ticks inward for compact look + axis.tick_params(which="both", direction="in", top=True, right=True) + + def plot_vacuum_vessel_and_divertor(axis, mfile_data, scan, colour_scheme): """Function to plot vacuum vessel and divertor boxes @@ -12555,6 +12703,7 @@ def main_plot( plot_plasma_effective_charge_profile(fig5.add_subplot(221), m_file_data, scan) plot_ion_charge_profile(fig5.add_subplot(223), m_file_data, scan) + plot_rad_contour(fig5.add_subplot(122), m_file_data, scan, imp) plot_fusion_rate_profiles(fig6.add_subplot(122), fig6, m_file_data, scan) From aa2eafffe5c4c39718c10a595ce61a84b86a30f7 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 1 Dec 2025 11:40:34 +0000 Subject: [PATCH 2/4] Enhance plasma profile calculations in plot_radprofile function --- process/io/plot_proc.py | 56 +++++++++++++++++++++++++---------------- 1 file changed, 35 insertions(+), 21 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 3a561b7e7a..183f28be0f 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -4248,6 +4248,24 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: mfile_data.data["f_nd_impurity_electrons(14)"].get_scan(scan), ]) + n_plasma_profile_elements = int( + mfile_data.data["n_plasma_profile_elements"].get_scan(scan) + ) + alphan = mfile_data.data["alphan"].get_scan(scan) + alphat = mfile_data.data["alphat"].get_scan(scan) + nd_plasma_electron_on_axis = mfile_data.data["nd_plasma_electron_on_axis"].get_scan( + scan + ) + temp_plasma_electron_on_axis_kev = mfile_data.data[ + "temp_plasma_electron_on_axis_kev" + ].get_scan(scan) + radius_plasma_pedestal_density_norm = mfile_data.data[ + "radius_plasma_pedestal_density_norm" + ].get_scan(scan) + radius_plasma_pedestal_temp_norm = mfile_data.data[ + "radius_plasma_pedestal_temp_norm" + ].get_scan(scan) + if i_plasma_pedestal == 0: # Intialise the radius rho = np.linspace(0, 1.0) @@ -4258,22 +4276,22 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: # The temperature profile te = te0 * (1 - rho**2) ** alphat + if i_plasma_pedestal == 0: + # Intialise the radius + rho = np.linspace(0, 1.0, n_plasma_profile_elements) + + # The density profile + ne = nd_plasma_electron_on_axis * (1 - rho**2) ** alphan + + # The temperature profile + te = temp_plasma_electron_on_axis_kev * (1 - rho**2) ** alphat + if i_plasma_pedestal == 1: - # Intialise the normalised radius - rhoped = ( - radius_plasma_pedestal_density_norm + radius_plasma_pedestal_temp_norm - ) / 2.0 - rhocore1 = np.linspace(0, 0.95 * rhoped) - rhocore2 = np.linspace(0.95 * rhoped, rhoped) - rhocore = np.append(rhocore1, rhocore2) - rhosep = np.linspace(rhoped, 1) - rho = np.append(rhocore, rhosep) + rho = np.linspace(0, 1, n_plasma_profile_elements) # The density and temperature profile - # done in such away as to allow for plotting pedestals - # with different radius_plasma_pedestal_density_norm and radius_plasma_pedestal_temp_norm - ne = np.zeros(rho.shape[0]) - te = np.zeros(rho.shape[0]) + ne = np.zeros_like(rho) + te = np.zeros_like(rho) for q in range(rho.shape[0]): if rho[q] <= radius_plasma_pedestal_density_norm: ne[q] = ( @@ -4285,9 +4303,7 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: else: ne[q] = nd_plasma_separatrix_electron + ( nd_plasma_pedestal_electron - nd_plasma_separatrix_electron - ) * (1 - rho[q]) / ( - 1 - min(0.9999, radius_plasma_pedestal_density_norm) - ) + ) * (1 - rho[q]) / (1 - radius_plasma_pedestal_density_norm) if rho[q] <= radius_plasma_pedestal_temp_norm: te[q] = ( @@ -4299,7 +4315,7 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: else: te[q] = temp_plasma_separatrix_kev + ( temp_plasma_pedestal_kev - temp_plasma_separatrix_kev - ) * (1 - rho[q]) / (1 - min(0.9999, radius_plasma_pedestal_temp_norm)) + ) * (1 - rho[q]) / (1 - radius_plasma_pedestal_temp_norm) # Intailise the radiation profile arrays pimpden = np.zeros([imp_data.shape[0], te.shape[0]]) @@ -4451,9 +4467,7 @@ def plot_rad_contour(axis, mfile_data, scan, impp): else: ne[q] = nd_plasma_separatrix_electron + ( nd_plasma_pedestal_electron - nd_plasma_separatrix_electron - ) * (1 - rho[q]) / ( - 1 - min(0.9999, radius_plasma_pedestal_density_norm) - ) + ) * (1 - rho[q]) / (1 - radius_plasma_pedestal_density_norm) if rho[q] <= radius_plasma_pedestal_temp_norm: te[q] = ( @@ -4465,7 +4479,7 @@ def plot_rad_contour(axis, mfile_data, scan, impp): else: te[q] = temp_plasma_separatrix_kev + ( temp_plasma_pedestal_kev - temp_plasma_separatrix_kev - ) * (1 - rho[q]) / (1 - min(0.9999, radius_plasma_pedestal_temp_norm)) + ) * (1 - rho[q]) / (1 - radius_plasma_pedestal_temp_norm) # Intailise the radiation profile arrays pimpden = np.zeros([imp_data.shape[0], te.shape[0]]) From ed03f98425941aac0d87c70d8dacafd3b1bf8cce Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 1 Dec 2025 13:09:15 +0000 Subject: [PATCH 3/4] Add checks for closed plasma boundary in radiation contour plots --- process/io/plot_proc.py | 36 ++++++++++++++++++++++++------------ 1 file changed, 24 insertions(+), 12 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 183f28be0f..0f629c3fb0 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -4266,16 +4266,6 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: "radius_plasma_pedestal_temp_norm" ].get_scan(scan) - if i_plasma_pedestal == 0: - # Intialise the radius - rho = np.linspace(0, 1.0) - - # The density profile - ne = ne0 * (1 - rho**2) ** alphan - - # The temperature profile - te = te0 * (1 - rho**2) ** alphat - if i_plasma_pedestal == 0: # Intialise the radius rho = np.linspace(0, 1.0, n_plasma_profile_elements) @@ -4542,6 +4532,16 @@ def plot_rad_contour(axis, mfile_data, scan, impp): axis.grid(True, which="major", linestyle="--", linewidth=0.8, alpha=0.7, zorder=1) axis.grid(True, which="minor", linestyle=":", linewidth=0.4, alpha=0.5, zorder=1) + props_core_reduce = {"boxstyle": "round", "facecolor": "khaki", "alpha": 0.8} + axis.text( + 0.02, + 0.02, + rf"$f_{{\text{{core,reduce}}}}$ = {1.0}", + transform=axis.transAxes, + fontsize=8, + verticalalignment="bottom", + bbox=props_core_reduce, + ) # make minor ticks visible on all sides and draw ticks inward for compact look axis.tick_params(which="both", direction="in", top=True, right=True) @@ -12717,14 +12717,26 @@ def main_plot( plot_plasma_effective_charge_profile(fig5.add_subplot(221), m_file_data, scan) plot_ion_charge_profile(fig5.add_subplot(223), m_file_data, scan) - plot_rad_contour(fig5.add_subplot(122), m_file_data, scan, imp) + + if m_file_data.data["i_plasma_shape"].get_scan(scan) == 1: + plot_rad_contour(fig5.add_subplot(122), m_file_data, scan, imp) + + i_shape = int(m_file_data.data["i_plasma_shape"].get_scan(scan)) + if i_shape != 1: + msg = ( + "Radiation contour plots require a closed (Sauter) plasma boundary " + "(i_plasma_shape == 1). " + f"Current i_plasma_shape = {i_shape}. Contour plots are skipped; " + "see the 1D radiation plots for available information." + ) + # Add explanatory text to both figures reserved for contour outputs + fig5.text(0.75, 0.5, msg, ha="center", va="center", wrap=True, fontsize=12) plot_fusion_rate_profiles(fig6.add_subplot(122), fig6, m_file_data, scan) if m_file_data.data["i_plasma_shape"].get_scan(scan) == 1: plot_fusion_rate_contours(fig7, fig8, m_file_data, scan) - i_shape = int(m_file_data.data["i_plasma_shape"].get_scan(scan)) if i_shape != 1: msg = ( "Fusion-rate contour plots require a closed (Sauter) plasma boundary " From 68ec297acd05685088c118dfb23d429683a961dc Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 5 Dec 2025 17:00:34 +0000 Subject: [PATCH 4/4] Requested changes --- process/io/plot_proc.py | 140 ++++++++++++++++++++++++---------------- 1 file changed, 83 insertions(+), 57 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 0f629c3fb0..7972ff6a22 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -19,6 +19,7 @@ import textwrap from argparse import RawTextHelpFormatter from importlib import resources +from typing import Any import matplotlib as mpl import matplotlib.backends.backend_pdf as bpdf @@ -4228,24 +4229,12 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: prof.set_title("Raw Data: Line & Bremsstrahlung radiation profile") # read in the impurity data - imp_data = read_imprad_data(2, impp) + imp_data = read_imprad_data(_skiprows=2, data_path=impp) # find impurity densities imp_frac = np.array([ - mfile_data.data["f_nd_impurity_electrons(01)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(02)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(03)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(04)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(05)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(06)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(07)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(08)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(09)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(10)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(11)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(12)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(13)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(14)"].get_scan(scan), + mfile_data.data[f"f_nd_impurity_electrons({i:02d})"].get_scan(scan) + for i in range(1, 15) ]) n_plasma_profile_elements = int( @@ -4266,9 +4255,10 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: "radius_plasma_pedestal_temp_norm" ].get_scan(scan) + rho = np.linspace(0, 1.0, n_plasma_profile_elements) + if i_plasma_pedestal == 0: # Intialise the radius - rho = np.linspace(0, 1.0, n_plasma_profile_elements) # The density profile ne = nd_plasma_electron_on_axis * (1 - rho**2) ** alphan @@ -4277,8 +4267,6 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: te = temp_plasma_electron_on_axis_kev * (1 - rho**2) ** alphat if i_plasma_pedestal == 1: - rho = np.linspace(0, 1, n_plasma_profile_elements) - # The density and temperature profile ne = np.zeros_like(rho) te = np.zeros_like(rho) @@ -4390,28 +4378,49 @@ def plot_radprofile(prof, mfile_data, scan, impp, demo_ranges) -> float: # --- -def plot_rad_contour(axis, mfile_data, scan, impp): - # read in the impurity data +def plot_rad_contour( + axis: "mpl.axes.Axes", mfile_data: "Any", scan: int, impp: int +) -> None: + """ + Plots the contour of line and bremsstrahlung radiation density for a plasma cross-section. + + This function reads impurity and plasma profile data, computes the radiation density profile, + interpolates it onto a 2D grid, and plots the upper and lower half contours on the provided axis. + + Parameters + ---------- + axis : matplotlib.axes.Axes + The matplotlib axis object to plot the contours on. + mfile_data : Any + Data object containing plasma and impurity profile information. + scan : int + The scan index to extract profile data for plotting. + impp : int + The impurity index to select the relevant impurity radiation data. + + Returns + ------- + None + This function modifies the provided axis in-place and does not return any value. + + Notes + ----- + - The function assumes the existence of several global or previously defined variables and functions, + such as `read_imprad_data`, `interp1d_profile`, and plasma pedestal parameters. + - The plotted contours represent the radiation density in units of MW.m^-3. + - The function adds colorbar, axis labels, title, and core reduction annotation to the plot. + """ + # Read in the impurity data imp_data = read_imprad_data(2, impp) + # imp data is a 3D array with shape (num_impurities, num_temp_points, (temp, lz, zav)) - # find impurity densities + # Find the relative number density of each impurity imp_frac = np.array([ - mfile_data.data["f_nd_impurity_electrons(01)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(02)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(03)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(04)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(05)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(06)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(07)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(08)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(09)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(10)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(11)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(12)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(13)"].get_scan(scan), - mfile_data.data["f_nd_impurity_electrons(14)"].get_scan(scan), + mfile_data.data[f"f_nd_impurity_electrons({i:02d})"].get_scan(scan) + for i in range(1, 15) ]) + # Get necessary plasma profile parameters n_plasma_profile_elements = int( mfile_data.data["n_plasma_profile_elements"].get_scan(scan) ) @@ -4430,23 +4439,27 @@ def plot_rad_contour(axis, mfile_data, scan, impp): "radius_plasma_pedestal_temp_norm" ].get_scan(scan) - if i_plasma_pedestal == 0: - # Intialise the radius - rho = np.linspace(0, 1.0, n_plasma_profile_elements) + # Initialize the radius + rho = np.linspace(0, 1.0, n_plasma_profile_elements) + # Re-construct te and ne profiles based on pedestal settings + # Parabolic profiles if no pedestal + if i_plasma_pedestal == 0: # The density profile ne = nd_plasma_electron_on_axis * (1 - rho**2) ** alphan # The temperature profile te = temp_plasma_electron_on_axis_kev * (1 - rho**2) ** alphat + # Profiles with pedestal if i_plasma_pedestal == 1: - rho = np.linspace(0, 1, n_plasma_profile_elements) - # The density and temperature profile + # Initiliase empty normalised array with zeros ne = np.zeros_like(rho) te = np.zeros_like(rho) + # Reconstruct the temperature and density profiles with pedestal for q in range(rho.shape[0]): + # Core density region if rho[q] <= radius_plasma_pedestal_density_norm: ne[q] = ( nd_plasma_pedestal_electron @@ -4455,10 +4468,12 @@ def plot_rad_contour(axis, mfile_data, scan, impp): ** alphan ) else: + # Pedestal density region ne[q] = nd_plasma_separatrix_electron + ( nd_plasma_pedestal_electron - nd_plasma_separatrix_electron ) * (1 - rho[q]) / (1 - radius_plasma_pedestal_density_norm) + # Core temperature region if rho[q] <= radius_plasma_pedestal_temp_norm: te[q] = ( temp_plasma_pedestal_kev @@ -4467,6 +4482,7 @@ def plot_rad_contour(axis, mfile_data, scan, impp): ** alphat ) else: + # Pedestal temperature region te[q] = temp_plasma_separatrix_kev + ( temp_plasma_pedestal_kev - temp_plasma_separatrix_kev ) * (1 - rho[q]) / (1 - radius_plasma_pedestal_temp_norm) @@ -4477,21 +4493,32 @@ def plot_rad_contour(axis, mfile_data, scan, impp): prad = np.zeros(te.shape[0]) # Intailise the impurity radiation profile - for k in range(te.shape[0]): - for i in range(imp_data.shape[0]): - if te[k] <= imp_data[i][0][0]: - lz[i][k] = imp_data[i][0][1] - elif te[k] >= imp_data[i][imp_data.shape[1] - 1][0]: - lz[i][k] = imp_data[i][imp_data.shape[1] - 1][1] + for rho in range(te.shape[0]): + # imp data is a 3D array with shape (num_impurities, num_temp_points, (temp, lz, zav)) + for impurity in range(imp_data.shape[0]): + # Check if profile temperature is lower than dataset minimum. + # If so, use the minimum loss function value + if te[rho] <= imp_data[impurity][0][0]: + lz[impurity][rho] = imp_data[impurity][0][1] + + # Check if profile temperature is higher than dataset maximum. + # If so, use the maximum loss function value + elif te[rho] >= imp_data[impurity][imp_data.shape[1] - 1][0]: + lz[impurity][rho] = imp_data[impurity][imp_data.shape[1] - 1][1] else: - # Use np.interp for log-log interpolation - log_te_data = np.log([row[0] for row in imp_data[i]]) - log_lz_data = np.log([row[1] for row in imp_data[i]]) - lz[i][k] = np.exp(np.interp(np.log(te[k]), log_te_data, log_lz_data)) - pimpden[i][k] = imp_frac[i] * ne[k] * ne[k] * lz[i][k] + # If profile valie is within dataset range, use log-log interpolation to find value for loss function + log_te_data = np.log([row[0] for row in imp_data[impurity]]) + log_lz_data = np.log([row[1] for row in imp_data[impurity]]) + lz[impurity][rho] = np.exp( + np.interp(np.log(te[rho]), log_te_data, log_lz_data) + ) + # Find the power density for each impurity at each rho + pimpden[impurity][rho] = ( + imp_frac[impurity] * ne[rho] * ne[rho] * lz[impurity][rho] + ) - for l in range(imp_data.shape[0]): # noqa: E741 - prad[k] = prad[k] + pimpden[l][k] * 1.0e-6 + for impurity in range(imp_data.shape[0]): + prad[rho] = prad[rho] + pimpden[impurity][rho] * 1.0e-6 p_rad_grid, r_grid, z_grid = interp1d_profile(prad, mfile_data, scan) @@ -12655,7 +12682,7 @@ def main_plot( "\033[91m Warning : Impossible to recover impurity data, try running the macro in the main/utility folder" ) print(" -> No impurity plot done\033[0m") - + i_shape = int(m_file_data.data["i_plasma_shape"].get_scan(scan)) # Setup params for text plots plt.rcParams.update({"font.size": 8}) @@ -12718,10 +12745,9 @@ def main_plot( plot_plasma_effective_charge_profile(fig5.add_subplot(221), m_file_data, scan) plot_ion_charge_profile(fig5.add_subplot(223), m_file_data, scan) - if m_file_data.data["i_plasma_shape"].get_scan(scan) == 1: + if i_shape == 1: plot_rad_contour(fig5.add_subplot(122), m_file_data, scan, imp) - i_shape = int(m_file_data.data["i_plasma_shape"].get_scan(scan)) if i_shape != 1: msg = ( "Radiation contour plots require a closed (Sauter) plasma boundary " @@ -12751,7 +12777,7 @@ def main_plot( plot_plasma_pressure_profiles(fig9.add_subplot(222), m_file_data, scan) plot_plasma_pressure_gradient_profiles(fig9.add_subplot(224), m_file_data, scan) # Currently only works with Sauter geometry as plasma has a closed surface - i_shape = int(m_file_data.data["i_plasma_shape"].get_scan(scan)) + if i_shape == 1: plot_plasma_poloidal_pressure_contours( fig9.add_subplot(121, aspect="equal"), m_file_data, scan