Create plasma current class

process/core/io/plot_proc.py

@@ -62,6 +62,7 @@
 )
 from process.models.physics.current_drive import ElectronBernstein, ElectronCyclotron
 from process.models.physics.impurity_radiation import read_impurity_file
+from process.models.physics.plasma_current import PlasmaCurrent
 from process.models.tfcoil.superconducting import SUPERCONDUCTING_TF_TYPES
 
 mpl.rcParams["figure.max_open_warning"] = 40
@@ -13406,23 +13407,24 @@ def main_plot(
     plot_ebw_ecrh_coupling_graph(figs[12].add_subplot(111), m_file, scan)
 
     plot_bootstrap_comparison(figs[13].add_subplot(221), m_file, scan)
-    plot_h_threshold_comparison(figs[13].add_subplot(224), m_file, scan)
+    PlasmaCurrent.plot_plasma_current_comparison(figs[13].add_subplot(224), m_file, scan)
+    plot_h_threshold_comparison(figs[14].add_subplot(224), m_file, scan)
     plot_density_limit_comparison(figs[14].add_subplot(221), m_file, scan)
-    plot_confinement_time_comparison(figs[14].add_subplot(224), m_file, scan)
+    plot_confinement_time_comparison(figs[15].add_subplot(224), m_file, scan)
 
-    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)
+    plot_debye_length_profile(figs[16].add_subplot(232), m_file, scan)
+    plot_velocity_profile(figs[16].add_subplot(233), m_file, scan)
+    plot_plasma_coloumb_logarithms(figs[16].add_subplot(231), m_file, scan)
 
-    plot_electron_frequency_profile(figs[15].add_subplot(212), m_file, scan)
+    plot_electron_frequency_profile(figs[16].add_subplot(212), m_file, scan)
 
-    plot_ion_frequency_profile(figs[16].add_subplot(311), m_file, scan)
+    plot_ion_frequency_profile(figs[17].add_subplot(311), m_file, scan)
 
-    plot_larmor_radius_profile(figs[16].add_subplot(313), m_file, scan)
+    plot_larmor_radius_profile(figs[17].add_subplot(313), m_file, scan)
 
     # Plot poloidal cross-section
     poloidal_cross_section(
-        figs[17].add_subplot(121, aspect="equal"),
+        figs[18].add_subplot(121, aspect="equal"),
         m_file,
         scan,
         demo_ranges,
@@ -13432,80 +13434,80 @@ def main_plot(
 
     # Plot toroidal cross-section
     toroidal_cross_section(
-        figs[17].add_subplot(122, aspect="equal"),
+        figs[18].add_subplot(122, aspect="equal"),
         m_file,
         scan,
         demo_ranges,
         colour_scheme,
     )
 
     # Plot color key
-    ax17 = figs[17].add_subplot(222)
+    ax17 = figs[18].add_subplot(222)
     ax17.set_position([0.5, 0.5, 0.5, 0.5])
     color_key(ax17, m_file, scan, colour_scheme)
 
-    ax18 = figs[18].add_subplot(211)
+    ax18 = figs[19].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[19].add_subplot(211)
+    ax185 = figs[20].add_subplot(211)
     ax185.set_position([0.1, 0.61, 0.8, 0.32])
 
-    ax18b = figs[19].add_subplot(212)
+    ax18b = figs[20].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)
 
     # 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[20].add_subplot(221, aspect="equal")
+        ax19 = figs[21].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[20])
+        plot_superconducting_tf_wp(ax19, m_file, scan, figs[21])
 
         # TF coil turn structure
-        ax20 = figs[21].add_subplot(325, aspect="equal")
+        ax20 = figs[22].add_subplot(325, aspect="equal")
         ax20.set_position([0.025, 0.5, 0.4, 0.4])
-        plot_tf_cable_in_conduit_turn(ax20, figs[21], m_file, scan)
-        plot_205 = figs[21].add_subplot(223, aspect="equal")
+        plot_tf_cable_in_conduit_turn(ax20, figs[22], m_file, scan)
+        plot_205 = figs[22].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[21], m_file, scan)
+        plot_cable_in_conduit_cable(plot_205, figs[22], m_file, scan)
     else:
-        ax19 = figs[20].add_subplot(211, aspect="equal")
+        ax19 = figs[21].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[20])
-        plot_resistive_tf_info(ax19, m_file, scan, figs[20])
+        plot_resistive_tf_wp(ax19, m_file, scan, figs[21])
+        plot_resistive_tf_info(ax19, m_file, scan, figs[21])
     plot_tf_coil_structure(
-        figs[22].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme
+        figs[23].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme
     )
 
-    plot_plasma_outboard_toroidal_ripple_map(figs[23], m_file, scan)
+    plot_plasma_outboard_toroidal_ripple_map(figs[24], m_file, scan)
 
-    plot_tf_stress(figs[24].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file)
+    plot_tf_stress(figs[25].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file)
 
-    plot_current_profiles_over_time(figs[25].add_subplot(111), m_file, scan)
+    plot_current_profiles_over_time(figs[26].add_subplot(111), m_file, scan)
 
     plot_cs_coil_structure(
-        figs[26].add_subplot(121, aspect="equal"), figs[26], m_file, scan
+        figs[27].add_subplot(121, aspect="equal"), figs[27], m_file, scan
     )
     plot_cs_turn_structure(
-        figs[26].add_subplot(224, aspect="equal"), figs[26], m_file, scan
+        figs[27].add_subplot(224, aspect="equal"), figs[27], m_file, scan
     )
 
     plot_first_wall_top_down_cross_section(
-        figs[27].add_subplot(221, aspect="equal"), m_file, scan
+        figs[28].add_subplot(221, aspect="equal"), 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_first_wall_poloidal_cross_section(figs[28].add_subplot(122), m_file, scan)
+    plot_fw_90_deg_pipe_bend(figs[28].add_subplot(337), m_file, scan)
 
-    plot_blkt_pipe_bends(figs[28], m_file, scan)
-    ax_blanket = figs[28].add_subplot(122, aspect="equal")
+    plot_blkt_pipe_bends(figs[29], m_file, scan)
+    ax_blanket = figs[29].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]")
@@ -13548,13 +13550,13 @@ def main_plot(
     )
 
     plot_main_power_flow(
-        figs[29].add_subplot(111, aspect="equal"), m_file, scan, figs[29]
+        figs[30].add_subplot(111, aspect="equal"), m_file, scan, figs[30]
     )
 
-    ax24 = figs[30].add_subplot(111)
+    ax24 = figs[31].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[30])
+    plot_system_power_profiles_over_time(ax24, m_file, scan, figs[31])
 
 
 def create_thickness_builds(m_file, scan: int):
@@ -13631,7 +13633,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(31)]
+    pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(32)]
 
     # run main_plot
     main_plot(

process/data_structure/physics_variables.py

@@ -961,6 +961,31 @@
 plasma_current: float = None
 """plasma current (A)"""
 
+c_plasma_peng_analytic: float = None
+"""Peng analytic plasma current (A)"""
+
+c_plasma_peng_double_null: float = None
+"""Peng double null divertor plasma current (A)"""
+
+c_plasma_cyclindrical: float = None
+"""Cylindrical plasma current (A)"""
+
+c_plasma_ipdg89: float = None
+"""ITER IPDG89 plasma current (A)"""
+
+c_plasma_todd_empirical_i: float = None
+"""Todd empirical plasma current I (A)"""
+
+c_plasma_todd_empirical_ii: float = None
+"""Todd empirical plasma current II (A)"""
+c_plasma_connor_hastie: float = None
+"""Connor-Hastie plasma current (A)"""
+
+c_plasma_sauter: float = None
+"""Sauter plasma current (A)"""
+
+c_plasma_fiesta_st: float = None
+"""FIESTA ST plasma current (A)"""
 
 p_plasma_neutron_mw: float = None
 """Neutron fusion power from just the plasma [MW]"""
@@ -1611,6 +1636,15 @@ def init_physics_variables():
         pflux_fw_rad_mw, \
         pden_ion_electron_equilibration_mw, \
         plasma_current, \
+        c_plasma_peng_analytic, \
+        c_plasma_peng_double_null, \
+        c_plasma_cyclindrical, \
+        c_plasma_ipdg89, \
+        c_plasma_todd_empirical_i, \
+        c_plasma_todd_empirical_ii, \
+        c_plasma_connor_hastie, \
+        c_plasma_sauter, \
+        c_plasma_fiesta_st, \
         p_plasma_neutron_mw, \
         p_neutron_total_mw, \
         pden_neutron_total_mw, \
@@ -1893,6 +1927,15 @@ def init_physics_variables():
     pflux_fw_rad_mw = 0.0
     pden_ion_electron_equilibration_mw = 0.0
     plasma_current = 0.0
+    c_plasma_peng_analytic = 0.0
+    c_plasma_peng_double_null = 0.0
+    c_plasma_cyclindrical = 0.0
+    c_plasma_ipdg89 = 0.0
+    c_plasma_todd_empirical_i = 0.0
+    c_plasma_todd_empirical_ii = 0.0
+    c_plasma_connor_hastie = 0.0
+    c_plasma_sauter = 0.0
+    c_plasma_fiesta_st = 0.0
     p_plasma_neutron_mw = 0.0
     p_neutron_total_mw = 0.0
     pden_neutron_total_mw = 0.0

process/main.py

@@ -102,6 +102,7 @@
     PlasmaBeta,
     PlasmaInductance,
 )
+from process.models.physics.plasma_current import PlasmaCurrent
 from process.models.physics.plasma_geometry import PlasmaGeom
 from process.models.physics.plasma_profiles import PlasmaProfile
 from process.models.power import Power
@@ -706,6 +707,7 @@ def __init__(self):
         self.plasma_bootstrap_current = PlasmaBootstrapCurrent(
             plasma_profile=self.plasma_profile
         )
+        self.plasma_current = PlasmaCurrent()
         self.physics = Physics(
             plasma_profile=self.plasma_profile,
             current_drive=self.current_drive,
@@ -714,6 +716,7 @@ def __init__(self):
             plasma_density_limit=self.plasma_density_limit,
             plasma_exhaust=self.plasma_exhaust,
             plasma_bootstrap_current=self.plasma_bootstrap_current,
+            plasma_current=self.plasma_current,
         )
         self.physics_detailed = DetailedPhysics(
             plasma_profile=self.plasma_profile,