search

ehpor / hcipy

A framework for performing optical propagation simulations, meant for high contrast imaging, in Python.
https://hcipy.org
MIT License
90 stars 30 forks source link

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

Open syhaffert opened 3 years ago

syhaffert commented 3 years ago

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 commented 3 years ago

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.