We assume $\psi = \rho^2$, so that $\rho$ is proportional to the minor radius. This is necessary for the Zernike polynomials to properly resolve the magnetic axis. However, it is not necessary far from the axis and might come at the cost of accuracy near the LCFS.
We should generalize the radial coordinate. One option could be $\psi = \rho^2 + \rho^3 - \rho^4$, so that it still has $\rho^2$ scaling near the axis but also $\psi \propto \rho$ near the boundary for better accuracy there.
