WIP:Update the heating and current drive docs to include all models

documentation/proc-pages/eng-models/heating_and_current_drive/NBI/culham_nb.md

@@ -0,0 +1,223 @@
+# Culham Neutral Beam Model | `culnbi()`
+
+- `iefrf/iefrffix` = 8 
+
+
+
+This routine calculates Neutral Beam current drive parameters
+using the corrections outlined in AEA FUS 172 to the ITER method.
+The result cannot be guaranteed for devices with aspect ratios far
+from that of ITER (approx. 2.8).
+
+| Output | Description |
+|----------|-------------|
+| $\mathtt{effnbss}$  | Neutral beam current drive efficiency in Amperes per Watt |
+| $\mathtt{fpion}$    | Fraction of NB power given to ions |
+| $\mathtt{fshine}$   | Shine-through fraction of the beam |
+
+$$
+\mathtt{frbeam} = \frac{R_{\text{tan}}}{R_0}
+$$
+
+Where $R_{\text{tan}}$ is major radius at which the centre-line of the beam is tangential to the toroidal direction. This can be user defined
+
+$$
+\left(1+ \frac{1}{A}\right) < \mathtt{frbeam}
+$$
+
+A quick sanity check is done to make sure no negative roots are formed when calculating $\mathtt{dpath}$ this prevents setups where the NBI beam would miss the plasma
+
+
+$$
+\mathtt{dpath} = R_0 \sqrt{\left(1+\frac{1}{A}\right)^2-\mathtt{frbeam}^2}
+$$
+
+Beams topping cross section is calculated via $\mathtt{sigbeam}$ found [here](../NBI/nbi_overview.md/#beam-stopping-cross-section-sigbeam). This produces $\mathtt{sigstop}$
+
+Calculate number of decay lengths to centre
+
+$$
+\mathtt{taubeam} = \mathtt{dpath} \times n_{\text{e,0}} \times \mathtt{sigstop}
+$$
+
+Calculate the shine through fraction of the beam
+
+$$
+\mathtt{fshine} = e^{\left(-2 \times \mathtt{dpath} \times n_{\text{e,0}} \times \mathtt{sigstop}\right)}
+$$
+
+Deuterium and tritium beam densities
+
+$$
+\mathtt{dend} = n_{\text{ion}} \times (1-\mathtt{ftritbm})
+$$
+
+$$
+\mathtt{dent} = n_{\text{ion}} \times \mathtt{ftritbm}
+$$
+
+Power split to the ions and electrons is clauclated with the $\mathtt{cfnbi()}$ method found [here](../NBI/nbi_overview.md/#ion-coupled-power-cfnbi) and outputs $\mathtt{fpion}$
+
+## Current drive efficiency | `etanb2()`
+
+This routine calculates the current drive efficiency in A/W of
+a neutral beam system, based on the 1990 ITER model,
+plus correction terms outlined in Culham Report AEA FUS 172.
+
+| Input       | Description                          |
+| :---------- | :----------------------------------- |
+| $\mathtt{abeam}$      | beam ion mass (amu)   |
+| $\mathtt{alphan}$, $\alpha_n$       | density profile factor   |
+| $\mathtt{alphat}$, $\alpha_T$       | temperature profile factor  |
+|  $\mathtt{aspect}$, $A$      |   aspect ratio                            |
+|  $\mathtt{dene}$, $n_{\text{e}}$     |    volume averaged electron density $(\text{m}^{-3})$                           |
+|  $\mathtt{dnla}$, $n_{\text{e,0}}$      |    line averaged electron density $(\text{m}^{-3})$                           |
+|  $\mathtt{enbeam}$      |  neutral beam energy $(\text{keV})$                             |
+|  $\mathtt{frbeam}$      |   R_tangent / R_major for neutral beam injection                            |
+|  $\mathtt{fshine}$      |  shine-through fraction of beam                             |
+|  $\mathtt{rmajor}$, $R$      |  plasma major radius $(\text{m})$                              |
+|  $\mathtt{rminor}$, $a$      |  plasma minor radius $(\text{m})$                             |
+|  $\mathtt{ten}$      |    density weighted average electron temperature $(\text{keV})$                             |
+|  $\mathtt{zeff}$, $Z_{\text{eff}}$     |   plasma effective charge                            |
+
+
+Charge of beam ions
+$$
+\mathtt{zbeam} = 1.0
+$$
+
+Fitting factor (IPDG89)
+
+$$
+\mathtt{bbd} = 1.0
+$$
+
+Volume averaged electron density ($10^{20} \text{m}^{-3}$)
+
+$$
+\mathtt{dene20} = n_{\text{e,20}}
+$$
+
+Line averaged electron density ($10^{20} \text{m}^{-3}$)
+
+$$    
+\mathtt{dnla20} = n_{\text{(e,0) 20}} 
+$$
+
+Critical energy ($\text{MeV}$) (power to electrons = power to ions) (IPDG89)
+N.B. ten is in keV
+
+$$
+\mathtt{ecrit} = 0.01 \times \mathtt{abeam} \times \mathtt{ten}
+$$
+
+Beam energy in MeV
+
+$$
+\mathtt{ebmev} = \frac{\mathtt{enbeam}}{10^3}
+$$
+
+x and y coefficients of function J0(x,y) (IPDG89)
+
+$$    
+\mathtt{xjs} = \frac{\mathtt{ebmev}}{\mathtt{bbd}\times \mathtt{ecrit}}  
+$$
+
+$$
+\mathtt{xj} = \sqrt{\mathtt{xjs}}
+$$
+
+$$
+\mathtt{yj} = \frac{0.8 \times Z_{\text{eff}}}{\mathtt{abeam}}
+$$
+
+Fitting function J0(x,y)
+
+$$
+\mathtt{j0} = \frac{xjs}{(4.0 + 3.0 \times \mathtt{yj} + \mathtt{xjs} \times (\mathtt{xj} + 1.39 + 0.61 \times yj^{0.7}))}
+$$
+
+Effective inverse aspect ratio, with a limit on its maximum value
+
+$$    
+ \mathtt{epseff} = \text{min}(0.2, (0.5 / A))
+$$   
+
+Reduction in the reverse electron current
+due to neoclassical effects
+
+$$
+\mathtt{gfac} = (1.55 + 0.85 / Z_{\text{eff}}) \times \sqrt{\mathtt{epseff}} - (0.2 + 1.55 / Z_{\text{eff}}) \times \mathtt{epseff}
+$$
+
+Reduction in the net beam driven current
+due to the reverse electron current
+
+$$
+\mathtt{ffac} = 1.0 - \frac{\mathtt{zbeam}}{Z_{\text{eff}}} \times (1.0 - \mathtt{gfac})
+$$
+
+Normalisation to allow results to be valid for
+non-ITER plasma size and density:
+
+Line averaged electron density ($10^{20} \text{m}^{-3}$) normalised to ITER
+
+$$    
+\mathtt{nnorm} = 1.0
+$$
+
+Distance along beam to plasma centre
+
+$$
+\mathtt{r} = \text{max}(R, R \times \mathtt{frbeam})
+$$
+
+$$
+\mathtt{eps1} = a / \mathtt{r}
+$$
+
+
+$$
+\mathtt{d} = R \times \sqrt{((1.0 + \mathtt{eps1})^2 - \mathtt{frbeam}^2)}
+$$
+
+Distance along beam to plasma centre for ITER
+assuming a tangency radius equal to the major radius
+
+$$
+\mathtt{epsitr} = 2.15 / 6.0
+$$  
+
+$$  
+\mathtt{dnorm} = 6.0 \times \sqrt{(2.0 \times \mathtt{epsitr} + \mathtt{epsitr}^2)}
+$$
+
+Normalisation to beam energy (assumes a simplified formula for
+the beam stopping cross-section)
+
+$$
+    \mathtt{ebnorm} = \mathtt{ebmev} \times ((\mathtt{nnorm} \times \mathtt{dnorm}) / (n_{\text{e,0}} \times \mathtt{d})) ^{1.0 / 0.78)}
+$$
+
+A_bd fitting coefficient, after normalisation with ebnorm
+
+$$   
+\mathtt{abd} = (
+    0.107
+    \times (1.0 - 0.35 \times \alpha_n + 0.14 \times \alpha_n^2)
+    \times (1.0 - 0.21 \times \alpha_T)
+    \times (1.0 - 0.2 \times \mathtt{ebnorm} + 0.09 \times \mathtt{ebnorm}^2)
+    )
+$$
+
+Normalised current drive efficiency ($\text{A/W} \text{m}^{2}$) (IPDG89)
+
+$$
+\mathtt{gamnb} = 5.0 \times \mathtt{abd} \times 0.1 \times \mathtt{ten} \times (1.0 - \mathtt{fshine}) \times \mathtt{frbeam} \times \frac{\mathtt{j0}}{0.2} \times \mathtt{ffac}
+$$
+
+Current drive efficiency (A/W)
+
+$$
+\text{Current drive efficiency [A/W]} = \frac{\mathtt{gamnb}}{\mathtt{dene20}\times R}
+$$

documentation/proc-pages/eng-models/heating_and_current_drive/NBI/iter_nb.md

@@ -0,0 +1,128 @@
+# ITER Neutral Beam Model | `iternb()`
+
+- `iefrf/iefrffix` = 5
+
+| Output | Description |
+|----------|-------------|
+| $\mathtt{effnbss}$  | Neutral beam current drive efficiency in $\text{A/W}$ |
+| $\mathtt{fpion}$    | Fraction of NB power given to ions |
+| $\mathtt{fshine}$   | Shine-through fraction of the beam |
+
+This model calculates the current drive parameters for a neutral beam system, based on the 1990 ITER model.[^1]
+
+
+Firstly the beam access is checked for such that
+$$
+\bigg(1+ \frac{1}{A}\bigg) > (R_{\text{tangential}}/R_0)
+$$
+
+The beam path length to centre is calculated:
+
+$$
+\underbrace{\mathtt{dpath}}_{\text{Path length to centre}} = R_0 \sqrt{\left(\left(1+\frac{1}{A}\right)^2-\mathtt{frbeam}^2\right)}
+$$
+
+
+Beam stopping cross-section ($\sigma_{\text{beam}}$) is calculated using the `sigbeam` method described [here](nbi_overview.md) :
+
+
+Calculate number of electron decay lengths to centre
+
+$$
+\tau_{\text{beam}} = \mathtt{dpath}\times n_e \sigma_{\text{beam}}
+$$
+
+Shine-through fraction of beam:
+$$
+f_{\text{shine}} = e^{(-2.0 \times  \mathtt{dpath} \times  n_e  \sigma_{\text{beam}})} \\
+$$
+
+Deuterium and tritium beam densities:
+$$
+n_D = n_i * (1.0 - \mathtt{ftritbm}) 
+$$
+
+$$
+n_T = n_i * \mathtt{ftritbm}
+$$
+
+Power split to ions / electrons is calculated via the the `cfnbi` method described [here](nbi_overview.md)
+
+
+## Current drive efficiency | `etanb()`
+
+This routine calculates the current drive efficiency of
+a neutral beam system, based on the 1990 ITER model.
+AEA FUS 251: A User's Guide to the PROCESS Systems Code
+ITER Physics Design Guidelines: 1989 IPDG89, N. A. Uckan et al,
+ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990
+
+| Input   | Description                                               |
+|---------|-----------------------------------------------------------|
+| $\mathtt{abeam}$, $m_{\text{u,ion}}$   | Beam ion mass ($\text{amu}$)                                       |
+| $\mathtt{alphan}$  | Density profile factor                                    |
+| $\mathtt{alphat}$  | Temperature profile factor                                |
+| $\mathtt{aspect}$, $A$  | Aspect ratio                                              |
+| $\mathtt{dene20}$, $n_{\text{e,20}}$    | Volume averaged electron density ($10^{20} \text{m}^{-3}$)                  |
+| $\mathtt{ebeam}$   | Neutral beam energy ($\text{keV}$)                                 |
+| $\mathtt{rmajor}$, R  | Plasma major radius ($\text{m}$)                                   |
+| $\mathtt{ten}$     | Density weighted average electron temperature ($\text{keV}$)       |
+| $\mathtt{zeff}$, $Z_{\text{eff}}$    | Plasma effective charge                                   |
+
+| Output  | Description                                               |
+|---------|-----------------------------------------------------------|
+| $\mathtt{etanb}$   | Neutral beam current drive efficiency in $\text{A/W}$ |
+
+
+
+
+$$
+\mathtt{zbeam} = 1.0
+$$
+
+$$
+\mathtt{bbd} = 1.0
+$$
+
+
+Ratio of E_beam/E_crit
+
+$$
+\mathtt{xjs} = \frac{\mathtt{ebeam}}{10 \  m_{\text{u,ion}} \ T_e} 
+$$
+
+$$
+\mathtt{xj} = \sqrt{\mathtt{xjs}}
+$$
+
+$$
+\mathtt{yj} = 0.8 \frac{Z_{\text{eff}}}{m_{\text{u,ion}}}
+$$
+
+$$
+\mathtt{rjfunc} = \frac{\mathtt{xjs}}{((4.0 + 3.0\mathtt{yj} + \mathtt{xjs}(\mathtt{xj} + 1.39 + 0.61\mathtt{yj}^{0.7})))}
+$$
+
+$$
+\mathtt{epseff} = \frac{0.5}{A}
+$$
+
+$$
+\mathtt{gfac} = \left(1.55 + \frac{0.85}{Z_{\text{eff}}}\right)\left(\sqrt{\mathtt{epseff}}-\left(0.2+\frac{1.55}{Z_{\text{eff}}}\right)\mathtt{epseff}\right)
+$$
+
+$$
+\mathtt{ffac} = \frac{1}{\mathtt{zbeam}} - \frac{(1-\mathtt{gfac})}{Z_{\text{eff}}}
+$$
+
+$$
+\mathtt{abd} = 0.107 (1-0.35 \ \mathtt{alphan}+0.14 \ \mathtt{alphan}^2)(1-0.21 \ \mathtt{alphat})(1-0.2\times 10^{-3}\mathtt{ebeam}+0.09\times 10^{-6} \ \mathtt{ebeam}^2)
+$$
+
+$$
+\text{Current drive efficiency [A/W]} = \mathtt{abd} \times\frac{5}{R_0} \times0.1\frac{\mathtt{ten}}{n_{\text{e},20}} \times \frac{\mathtt{rjfunc}}{0.2}\mathtt{ffac}
+$$
+
+
+
+[^1]: N. A. Uckan and ITER Physics Group, *"ITER Physics Design Guidelines: 1989"*, ITER Documentation Series, No. 10, IAEA/ITER/DS/10 (1990)