Consider ways to improve accuracy/performance of 3D axial derivative interpolation

[leon-vv]

Ideas to improve the method:

• $A^{\nu+1}_m$ should be the second derivative (in z direction) of coefficient $A^{\nu}_m$. We could use this to improve the accuracy of the cubic spline interpolation along the z axis.
• The current theta interpolation uses values computed using numerical integration in Python using 64 bit floating point numbers. Therefore they can be expected to be accurate to at most a small multiple of the 64 bit float machine epsilon. Since for the pre-computation we can expend an arbitrary amount of computing power, we need to make sure these values are correct down to the max precision achievable with 64 bit floating point numbers. We can achieve this using arbitrary precision floating point number numerical integration.
• The current theta interpolation uses a cubic spline interpolation. We can consider other interpolation types, also we can make the cubic spline interpolation more accurate by explicitly providing the derivatives, which we can compute to arbitrary precision ourselves.

[leon-vv]

Regarding point two, I made an implementation of computing the coefficients using arbitrary precision floating point numbers (Python's mpmath ). It was extremely slow, taking hours for the computation of a thousand theta values. Anyway, I checked that the current fast computations using Numba give results accurate to better than 100 epsilon, usually down to a few epsilon. The gain in using arbitrary precision math is therefore likely to be small..

[leon-vv]

See also #24

[leon-vv]

Closing for now, as #24 already tracks the implementation of improved axial interpolation of the coefficients, while the 64 bit floating point arithmetic has been shown not to be a great source of error (compared with arbitrary precision floating point math)