No need for series expansion and eliminating divisions with new basis.

[ericagol]

@rodluger I think I may have figured out a stable recursion formula for the cases in which we found that we had series expansion for v_{max} and then iterate downwards. If the recursion is stable, there may be no need for a series expansion! The basic idea is simple: substitute \sin^2{\xi} = 1-\cos^2{\xi}, so that integrals are over \cos^{2v}\xi rather than \sin^2{\xi}. It looks like the reduction formulae for these are additive, so we don't get the cancellation which leads to numerical instability. I feel a bit embarrassed for not having figured this out for the starry paper, but such is the nature of research I guess.

I also may have figured out a simpler basis set whose coefficients are merely binomial coefficients, which can be computed via addition using Pascal's triangle, which will get rid of the need for lots of divisions. This may also be more numerically stable since I think the terms in the series all have the same sign. Instead of f(z) = z^n, the new basis is f(z) = (n+2)^{-1} \sum_{i=1}^{n/2} z^{2i} for even n, which looks more complicated, but actually may be easier to compute. For odd n, it becomes f(z) = \frac{z^3}{n+2} \sum_{i=1}^{(n-3)/2} z^{2i}. These have \tilde g terms of z^n for even n, and z^n-3/(n+2)z for odd n. So, this should give very little algebra for the conversion between these coefficients and the u_n or even less for Gimenez' a_n.

I'm wondering whether I should try to derive these with your SymPy notebooks so that it will make it easy to add these to paper.

[rodluger]

@ericagol Awesome! So this applies to both starry and limbdark? That will speed up the code quite a bit I imagine. And yes, SymPy derivations would be awesome. Let me know if you need help with those.

[ericagol]

@rodluger Yes, I think it should apply to starry as well. I'm still not yet sure if this will improve things, but if it does, I suspect it will be faster and more accurate.