Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Divide by zero in make_zernike_basis with grids that include 0 as grid point. #91

Open
syhaffert opened this issue Jan 4, 2021 · 1 comment
Labels
bug Something isn't working

Comments

@syhaffert
Copy link
Collaborator

The make_zernike_basis has a divide by zero error when used with a grid that has 0 as grid point.
The specific line that is causing the issue is line 139 in zernike.py in the function zernike_radial.

h3 = -4 * (q - 2) * (q - 3) / float((p + q - 2) * (p - q + 4))
h2 = h3 * (p + q) * (p - q + 2) / float(4 * (q - 1)) + (q - 2)
h1 = q * (q - 1) / 2.0 - q * h2 + h3 * (p + q + 2) * (p - q) / 8.0

r2 = zernike_radial(2, 2, r, cache)
res = h1 * zernike_radial(p, q, r, cache) + (h2 + h3 / r2) * zernike_radial(n, q - 2, r, cache)

The line with h3 / r2 is creates a NaN if zero is included as grid point.
This does not happen analytically because the limit r->0 converges, just like the sinc function.
We evaluate all polynomials numerically so there are no analytical formulas for the zero point.
I have not read the paper yet that the documentation quotes. Maybe they propose a solution for this.

A quick and dirty fix is to shift the grid by a small amount. I used a 1 um shift for a grid of 10m in diameter.

@syhaffert syhaffert added the bug Something isn't working label Jan 4, 2021
@syhaffert
Copy link
Collaborator Author

To circumvent this issue we can switch to the recurrence relation of T. Andersen 2018 (https://www.osapublishing.org/oe/fulltext.cfm?uri=oe-26-15-18878&id=395240)

The new recurrence relation does not include 1/r^2 terms because it starts at lowest order and works towards the higher order modes. I will implement the new recurrence relation and compare the two.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

1 participant