3388 update the beta limit section of the docs to be up to date and i…

…nclude all models (#3389)

* Refactor culblm function to calculate_beta_limit with enhanced documentation and type hints

* Update navigation structure in mkdocs.yml for Beta Limit section. Move the fast alpha pressure contribution into the plasma beta docs section

* Rename betaft to beta_fast_alpha across the codebase for consistency and clarity

* Rename switch variable iculbl to i_beta_component for clarity in beta limit scaling

* Rename betalim to beta_limit_upper for clarity and consistency across the codebase

* Rename betalim_lower to beta_limit_lower for consistency and clarity across the codebase

* Initial input of beta type constraint equations to docs

* Rename betpmx to beta_poloidal_max for clarity and consistency across the codebase

* Rename fbetap to fbeta_poloidal for clarity and consistency across the codebase

* Rename betap to beta_poloidal for clarity and consistency across the codebase

* Rename betat to beta_toroidal and add it as a proper physics variable for clarity and consistency across the codebase

* Add comment for normalised total beta calculation and update output label for clarity

* Rename betath to beta_thermal and add it as a proper physics variable for consistency across the codebase

* Add beta_thermal_poloidal and beta_thermal_toroidal variables for thermal beta calculations

* Add norm_beta_thermal variable and update output for normalised thermal beta calculation

* Rename normalised_total_beta to norm_beta_total and add norm_beta_toroidal and norm_beta_poloidal to update related calculations and references for consistency

* Enhance plasma beta documentation and update fast alpha beta calculation notes for clarity

* Add more detail for the alpha beta second model and add graph

* Rename 'betap' to 'beta_poloidal' in BootstrapFraction parameter classes for clarity to account for new bootstrap scalings

* 📝 Update plasma beta documentation to clarify relationships and constraints

* Rename norm_beta_thermal to beta_norm_thermal in Python and Fortran files for consistency

* Rename norm_beta_toroidal to beta_norm_toroidal in Python and Fortran files for consistency

* Rename norm_beta_poloidal to beta_norm_poloidal in Python and Fortran files for consistency

* Rename norm_beta_total to beta_norm_total in Python and Fortran files for consistency

* 🔄 Rename dnbeta to beta_norm_limit_upper for consistency across data files and scripts

* 🔄 Rename gammaft to f_beta_alpha_beam_thermal for consistency in Fortran and data files

* 🔄 Rename epbetmax to beta_poloidal_eps_max for consistency across files and documentation

* 🔄 Add beta_poloidal_eps calculation and documentation for poloidal beta limit

* 🔄 Rename fbetatry to fbeta_max for consistency across data files and documentation

* 🔄 Rename ifalphap to i_beta_fast_alpha for consistency across data files and scripts

* 🔄 Rename beta_limit_upper to beta_max for consistency across codebase and documentation

* 🔄 Rename beta_norm_limit_upper to beta_norm_max for consistency across data files and scripts

* 🔄 Update obsolete_vars.py to rename variables for consistency across the codebase

* 🔄 Rename beta_limit_lower to beta_min for consistency across the codebase

* 🔄 Rename fbetatry_lower to fbeta_min for consistency across the codebase

* 🐛 Update minimum beta constraint to work on the total beta and not thermal

* 🔄 Rename fbeta to fbeta_poloidal_eps for consistency across the codebase

* 🎨 Update plasma_beta.md for clarity and consistency in beta limit descriptions

* 🐛 Fix Tholerus pressure peaking parameter and expand docs

* Add Menard normalized beta limit plotting script

* Add original normalized beta limit plotting script and fix spacing in Menard plot

* Add beta_poloidal_max and beta_min to the output file

* ❇️ Add calculation for plasma thermal energy derived from thermal beta and create new variable e_plasma_beta_thermal

* Rename total_plasma_internal_energy to e_plasma_beta and update references throughout the code

* Update output descriptions for plasma thermal energy and total beta energy calculations

* 🎨 Refactor plasma thermal energy calculation and remove duplicate code

* Update bokeh version requirement to 3.4.0 in development dependencies

documentation/proc-pages/development/standards.md

@@ -533,15 +533,15 @@ subroutine constraint_eqn_001(args)
   !! - \( T_i \) -- density weighted average ion temperature [keV]
   !! - \( B_{tot} \) -- total toroidal + poloidal field [T]
 
-  use physics_variables, only: betaft, beta_beam, dene, ten, dnitot, tin, btot, beta
+  use physics_variables, only: beta_fast_alpha, beta_beam, dene, ten, dnitot, tin, btot, beta
   use constants, only: electron_charge,rmu0
 
   implicit none
 
   type(constraint_args_type), intent(out) :: args
   !! constraint derived type
 
-    args%cc = 1.0D0 - (betaft + beta_beam + &
+    args%cc = 1.0D0 - (beta_fast_alpha + beta_beam + &
       2.0D3*rmu0*electron_charge * (dene*ten + dnitot*tin)/btot**2 )/beta
     args%con = beta * (1.0D0 - args%cc)
     args%err = beta * args%cc

documentation/proc-pages/fusion-devices/stellarator.md

@@ -123,7 +123,7 @@ Stellarators try to achieve zero plasma current in order to allow safe divertor
 ### Beta limit
 
 The stellarator version calculates the plasma beta based on the input parameter and it is thus not necessary to Differently to the tokamak version, 
-The beta limit is assumed to be 5%, based on 3-D MHD calculations[^7]. To apply the beta limit, constraint equation no. 24 should be turned on with iteration variable no. 36 (`fbetatry`).
+The beta limit is assumed to be 5%, based on 3-D MHD calculations[^7]. To apply the beta limit, constraint equation no. 24 should be turned on with iteration variable no. 36 (`fbeta_max`).
 
 ### Density limit
 

documentation/proc-pages/physics-models/error.txt

@@ -388,9 +388,9 @@ where \(B_0\) is the axial vacuum toroidal field, and \(\beta\) is
 defined with respect to the total equilibrium \(\mathbf{B}\)-field
 \footnote{T. C. Hender et al., `Physics Assessment for the European
   Reactor Study', AEA Fusion Report AEA FUS 172 (1992)}. The beta
-coefficient \(g\) is set using input parameter \texttt{dnbeta}. To apply
+coefficient \(g\) is set using input parameter \texttt{beta_norm_max}. To apply
 the beta limit, constraint equation no. 24 should be turned on with
-iteration variable no. 36 (\texttt{fbetatry}). The limit can be applied
+iteration variable no. 36 (\texttt{fbeta_max}). The limit can be applied
 to either the total plasma beta, in which case switch \texttt{iculbl}
 should be set to 0, to only the thermal component of the plasma beta, in
 which case \texttt{iculbl} should be set to 1, or to the thermal plus
@@ -401,7 +401,7 @@ neutral beam components, in which case \texttt{iculbl} should be set to
 coefficient}{Aspect ratio scaling of beta g coefficient}}\label{aspect-ratio-scaling-of-beta-g-coefficient}
 
 Switch \texttt{gtscale} determines whether the beta \(g\) coefficient
-\texttt{dnbeta} should scale with aspect ratio
+\texttt{beta_norm_max} should scale with aspect ratio
 (\(\mathtt{gtscale \not= 0}\)), or be fixed at the input value
 (\texttt{gtscale\ =\ 0\}}). Note that \texttt{gtscale} is over-ridden if
 \texttt{iprofile} = 1.
@@ -411,14 +411,14 @@ Switch \texttt{gtscale} determines whether the beta \(g\) coefficient
 
 To apply a limit to the value of \(\epsilon.\beta_p\), where
 \(\epsilon = a/R\) is the inverse aspect ratio, constraint equation no.
-6 should be turned on with iteration variable no. 8 (\texttt{fbeta}).
+6 should be turned on with iteration variable no. 8 (\texttt{fbeta_poloidal_eps}).
 The limiting value of \(\epsilon.\beta_p\) should be set using input
-parameter \texttt{epbetmax}.
+parameter \texttt{beta_poloidal_eps_max}.
 
 \subsection{Fast Alpha Pressure
 Contribution}\label{fast-alpha-pressure-contribution}
 
-Switch \texttt{ifalphap} may be used to select the model used to
+Switch \texttt{i_beta_fast_alpha} may be used to select the model used to
 calculate the pressure contribution from the fast alpha particles, there
 are two options 1\footnote{N.A. Uckan and ITER Physics Group, `ITER
   Physics Design Guidelines: 1989', ITER Documentation Series, No. 10,
@@ -428,10 +428,10 @@ are two options 1\footnote{N.A. Uckan and ITER Physics Group, `ITER
 \[\begin{aligned}
 \frac{\beta_{\alpha}}{\beta_{th}} & = 0.29 \, \left( \langle T_{10} \rangle -
   0.37 \right) \, \left( \frac{n_{DT}}{n_e} \right)^2
-\hspace{20mm} \mbox{ifalphap = 0} \\
+\hspace{20mm} \mbox{i_beta_fast_alpha = 0} \\
 \frac{\beta_{\alpha}}{\beta_{th}} & = 0.26 \, \left( \langle T_{10} \rangle -
   0.65 \right)^{0.5} \, \left( \frac{n_{DT}}{n_e} \right)^2
-\hspace{16mm} \mbox{ifalphap = 1 (default)}
+\hspace{16mm} \mbox{i_beta_fast_alpha = 1 (default)}
 \end{aligned}\]
 
 The latter model is a better estimate at higher temperatures.
@@ -568,7 +568,7 @@ profile peaking factor \texttt{alphaj} is consistent with the input
 values for the safety factor on axis and at the plasma edge (\texttt{q0}
 and \texttt{q}, respectively), the plasma internal inductance \(l_i\) is
 consistent with this \texttt{alphaj}, and the beta \(g\) coefficient
-\texttt{dnbeta} scales with \(l_i\).
+\texttt{beta_norm_max} scales with \(l_i\).
 
 It is recommended that current scaling law \texttt{icurr\ =\ 4} is used
 if \texttt{iprofile\ =\ 1}. Switch \texttt{gtscale} is over-ridden if

documentation/proc-pages/physics-models/plasma_current/plasma_current.md

@@ -556,7 +556,7 @@ $$
 
 A limited degree of self-consistency between the plasma current profile and other parameters can be 
 enforced by setting switch `iprofile = 1`. This sets the current 
-profile peaking factor $\alpha_J$ (`alphaj`),  the normalised internal inductance $l_i$ (`rli`) and beta limit $g$-factor (`dnbeta`) using the 
+profile peaking factor $\alpha_J$ (`alphaj`),  the normalised internal inductance $l_i$ (`rli`) and beta limit $g$-factor (`beta_norm_max`) using the 
 safety factor on axis `q0` and the cylindrical safety factor $q_*$ (`qstar`):   
 
 $$\begin{aligned}

documentation/proc-pages/scripts/plotting_scripts/menard_beta_norm_plot.py

@@ -0,0 +1,24 @@
+import numpy as np
+from bokeh.models import ColumnDataSource
+from bokeh.plotting import figure, output_file, save
+
+
+x = np.linspace(1.0, 5, 500)
+y = 3.12 + 3.5 * (1 / x) ** 1.7
+source = ColumnDataSource(data=dict(x=x, y=y))
+
+plot = figure(
+    x_range=(1, 5),
+    y_range=(2, 8),
+    width=400,
+    height=400,
+    title="Menard Normalized Beta Limit",
+)
+plot.xaxis.axis_label = r"Aspect ratio, \  $$[A]$$"
+plot.yaxis.axis_label = r"Normalized beta limit, \  $$[\beta_N]$$"
+
+plot.line("x", "y", source=source, line_width=3, line_alpha=0.6)
+
+# Save the plot as HTML
+output_file("menard_beta_norm.html", title="Menard Normalized Beta Limit")
+save(plot)

documentation/proc-pages/scripts/plotting_scripts/original_beta_norm_plot.py

@@ -0,0 +1,24 @@
+import numpy as np
+from bokeh.models import ColumnDataSource
+from bokeh.plotting import figure, output_file, save
+
+
+x = np.linspace(1.0, 5, 500)
+y = 2.7 * (1 + 5 * (1 / x) ** 3.5)
+source = ColumnDataSource(data=dict(x=x, y=y))
+
+plot = figure(
+    x_range=(1, 5),
+    y_range=(2, 15),
+    width=400,
+    height=400,
+    title="Original Normalized Beta Limit",
+)
+plot.xaxis.axis_label = r"Aspect ratio, \  $$[A]$$"
+plot.yaxis.axis_label = r"Normalized beta limit, \  $$[\beta_N]$$"
+
+plot.line("x", "y", source=source, line_width=3, line_alpha=0.6)
+
+# Save the plot as HTML
+output_file("original_beta_norm.html", title="Original Normalized Beta Limit")
+save(plot)

examples/data/csv_output_large_tokamak_MFILE.DAT

@@ -343,23 +343,23 @@
  Safety_factor_at_95%_flux_surface_______________________________________ (q95)_________________________      3.4251E+00 ITV
  Cylindrical_safety_factor_(qcyl)________________________________________ (qstar)_______________________      2.8387E+00    
  Total_plasma_beta_______________________________________________________ (beta)________________________      3.6055E-02 ITV
- Total_poloidal_beta_____________________________________________________ (betap)_______________________      1.2531E+00 OP 
+ Total_poloidal_beta_____________________________________________________ (beta_poloidal)_______________________      1.2531E+00 OP 
  Total_toroidal_beta_____________________________________________________ ______________________________      3.7123E-02 OP 
  Fast_alpha_beta_________________________________________________________ (betaft)______________________      4.4814E-03 OP 
  Beam_ion_beta___________________________________________________________ (betanb)______________________      0.0000E+00 OP 
- (Fast_alpha_+_beam_beta)/(thermal_beta)_________________________________ (gammaft)_____________________      1.4194E-01 OP 
+ (Fast_alpha_+_beam_beta)/(thermal_beta)_________________________________ (f_beta_alpha_beam_thermal)_____________________      1.4194E-01 OP 
  Thermal_beta____________________________________________________________ ______________________________      3.1574E-02 OP 
  Thermal_poloidal_beta___________________________________________________ ______________________________      1.0974E+00 OP 
  Thermal_toroidal_beta_(=_beta-exp)______________________________________ ______________________________      3.2509E-02 OP 
- 2nd_stability_beta_:_beta_p_/_(R/a)_____________________________________ (eps*betap)___________________      4.1770E-01    
- 2nd_stability_beta_upper_limit__________________________________________ (epbetmax)____________________      1.3800E+00    
- Beta_g_coefficient______________________________________________________ (dnbeta)______________________      4.7593E+00    
+ 2nd_stability_beta_:_beta_p_/_(R/a)_____________________________________ (eps*beta_poloidal)___________________      4.1770E-01    
+ 2nd_stability_beta_upper_limit__________________________________________ (beta_poloidal_eps_max)____________________      1.3800E+00    
+ Beta_g_coefficient______________________________________________________ (beta_norm_max)______________________      4.7593E+00    
  Normalised_thermal_beta_________________________________________________ ______________________________      2.5527E+00    
  Normalised_total_beta___________________________________________________ ______________________________      2.9150E+00    
  Normalised_toroidal_beta________________________________________________ (normalised_toroidal_beta)____      3.0013E+00    
  Limit_on_thermal_beta___________________________________________________ (betalim)_____________________      5.8867E-02    
  Plasma_thermal_energy_(J)_______________________________________________ ______________________________      9.3900E+08 OP 
- Total_plasma_internal_energy_(J)________________________________________ (total_plasma_internal_energy)      1.0723E+09 OP 
+ Total_plasma_internal_energy_(J)________________________________________ (e_plasma_beta)      1.0723E+09 OP 
  Electron_temperature_(keV)______________________________________________ (te)__________________________      1.2238E+01 ITV
  Electron_temperature_on_axis_(keV)______________________________________ (te0)_________________________      2.5063E+01    
  Ion_temperature_(keV)___________________________________________________ (ti)__________________________      1.2238E+01    
@@ -1559,7 +1559,7 @@ alphan = 1.00
 alphat = 1.45 
 
 * (troyon-like) coefficient for beta scaling
-dnbeta = 3.0
+beta_norm_max = 3.0
 
 * Zohm elongation scaling adjustment factor (ishape=2; 3)
 fkzohm = 1.02
@@ -1571,7 +1571,7 @@ gamma = 0.3
 i_bootstrap_current = 4
 
 * Switch for beta limit scaling
-iculbl = 1
+i_beta_component = 1
 
 * Switch for plasma current scaling
 i_plasma_current    = 4
@@ -1580,7 +1580,7 @@ i_plasma_current    = 4
 i_density_limit = 7
 
 * Switch for fast alpha pressure calculation
-ifalphap = 1
+i_beta_fast_alpha = 1
 
 * Switch for inverse quadrature in l-mode scaling laws 5 and 9
 iinvqd = 1