Add robust Zernike evaluation

Author: jmeyers314

galsim/zernike.py

@@ -721,17 +721,40 @@ def gradY(self):
         newCoef /= self.R_outer
         return Zernike(newCoef, R_outer=self.R_outer, R_inner=self.R_inner)
 
-    def __call__(self, x, y):
+    def __call__(self, x, y, robust=False):
         """Evaluate this Zernike polynomial series at Cartesian coordinates x and y.
         Synonym for `evalCartesian`.
 
+        Parameters:
+            x:       x-coordinate of evaluation points.  Can be list-like.
+            y:       y-coordinate of evaluation points.  Can be list-like.
+            robust:  If True, use a more robust method for evaluating the polynomial.
+                     This can sometimes be slower, but is usually more accurate,
+                     especially for large Noll indices.  [default: False]
+        Returns:
+            Series evaluations as numpy array.
+        """
+        if robust:
+            return self.evalCartesianRobust(x, y)
+        return self.evalCartesian(x, y)
+
+    def evalCartesianRobust(self, x, y):
+        """Evaluate this Zernike polynomial series at Cartesian coordinates x and y using a more
+        robust method than the default `evalCartesian`.
+
         Parameters:
             x:    x-coordinate of evaluation points.  Can be list-like.
             y:    y-coordinate of evaluation points.  Can be list-like.
+
         Returns:
             Series evaluations as numpy array.
         """
-        return self.evalCartesian(x, y)
+        x = np.asarray(x)
+        y = np.asarray(y)
+        rho = (x + 1j*y) / self.R_outer
+        rhosq = np.abs(rho)**2
+        ar = _noll_coef_array(len(self.coef)-1, self.R_inner/self.R_outer).dot(self.coef[1:])
+        return horner2d(rhosq, rho, ar).real
 
     def evalCartesian(self, x, y):
         """Evaluate this Zernike polynomial series at Cartesian coordinates x and y.
@@ -743,7 +766,8 @@ def evalCartesian(self, x, y):
         Returns:
             Series evaluations as numpy array.
         """
-        return horner2d(x, y, self._coef_array_xy, dtype=float)
+        ar = self._coef_array_xy
+        return horner2d(x, y, ar, dtype=float)
 
     def evalPolar(self, rho, theta):
         """Evaluate this Zernike polynomial series at polar coordinates rho and theta.

tests/test_zernike.py

@@ -1583,25 +1583,25 @@ def test_large_j(run_slow):
     eps = R_inner/R_outer
 
     test_vals = [
-        (10, 1e-12),  # Z66
-        (20, 1e-12),  # Z231
-        (40, 1e-9),   # Z861
+        (10, 1e-12, 1e-12),  # Z66
+        (20, 1e-12, 1e-12),  # Z231
+        (40, 1e-9, 1e-12),   # Z861
     ]
     if run_slow:
         test_vals += [
-            (60, 1e-6),   # Z1891
-            (80, 1e-3),  # Z3321
-            # (100, 10),  # Z5151  # This one is catastrophic failure!
+            (60, 1e-6, 1e-12),  # Z1891
+            (80, 1e-3, 1e-12),  # Z3321
+            (100, None, 1e-11),   # Z5151
+            (200, None, 1e-11),   # Z20301
         ]
 
     print()
-    for n, tol in test_vals:
+    for n, tol, tol2 in test_vals:
         j = (n+1)*(n+2)//2
         _, m = galsim.zernike.noll_to_zern(j)
         print(f"Z{j} => (n, m) = ({n}, {m})")
         assert n == abs(m)
-        coefs = np.zeros(j+1)
-        coefs[j] = 1.0
+        coefs = [0]*j+[1]
         zk = Zernike(coefs, R_outer=R_outer, R_inner=R_inner)
 
         def analytic_zk(x, y):
@@ -1613,10 +1613,23 @@ def analytic_zk(x, y):
             else:
                 return r**n * np.sin(n*theta) * factor
 
+        analytic_vals = analytic_zk(x, y)
+        if n < 100:
+            np.testing.assert_allclose(
+                zk(x, y),
+                analytic_vals,
+                atol=tol, rtol=tol
+            )
+
+        robust_vals = zk.evalCartesianRobust(x, y)
         np.testing.assert_allclose(
-            zk(x, y),
-            analytic_zk(x, y),
-            atol=tol, rtol=tol
+            robust_vals,
+            analytic_vals,
+            atol=tol2, rtol=tol2
+        )
+        np.testing.assert_equal(
+            robust_vals,
+            zk(x, y, robust=True)
         )