OceanNuclear
commented
1 week ago After a detailed check of the maths, I can confirm that this is indeed a mistake.
I have plotted
However, the plot below
, generated using the data @timothy-nunn gave me, unfortunately happens to use data whose final term (the term mistakenly raised to $N+1+1$) happen to have coefficient
sum(args)$=0$, so the error is invisible in this case.
Regardless, I am pretty confident that @timothy-nunn is right, so I have implemented the change.
Describe the bug
Documentation issue
The equation (14) in
equilibria.rstprovides the Lao polynomial parameterisation for the plasma profiles.The equation, however, disagrees with Lao et al. [1] and the docstring in
bluemira/equilibria/profiles/laopoly. I believe that the equation should instead read: $g\big(\overline{\psi},\boldsymbol{\alpha}\big) = \sum_{n=0}^{N}\alpha_{n+1}\overline{\psi}^{~n}-\overline{\psi}^{~N+1}\sum_{n=0}^{N} \alpha_{n+1}$.That is, the first normalised $\psi$ should be to the power of $n$, not $n+1$.
However, rewriting it in keeping with Lao et al. may reduce confusion further: $
g\big(\overline{\psi},\boldsymbol{\alpha}\big) = \sum_{n=0}^{N}\alpha_{n}\overline{\psi}^{~n}-\overline{\psi}^{~N+1}\sum_{n=0}^{N} \alpha_{n}$. Where $N$ is the number of coefficients -1: $N = |\alpha| -1$.e.g. for $\alpha=[1, 2]$, $N=2-1=1$ and the equation becomes $g(x) = (1x^0 + 2x^1) - x^2(1+2)$
Implementation issue
The final term of
bluemira/equilibria/profiles/laopolyraises the normalised equilibrium toN+2(see the above definition) but I believe it should be toN+1.Because N is defined as the number of constants minus one, this final term of the equation raises the normalised psi to the power of the number of alpha constants. However, the implementation does
x ** (len(args) + 1)which is essentially raising it to $N+1 +1$.[1] Lao, L. L., et al. "Reconstruction of current profile parameters and plasma shapes in tokamaks." Nuclear fusion 25.11 (1985): 1611.