jmeyers314
commented
3 years ago To start, Mike recently wrote on another issue:
Oh, one more comment about the conventions in Zernike. I also think some of those might have the 1j in the wrong place, which might explain some of the convention mismatches you were seeing. In particular in _xy_contribution the 1j multiplies the j term, which I think is backwards. I didn't try switching it, but I'd guess it would break some existing unit tests, so you might try to think through whether that should be changed or not.
I looked into this a bit more this afternoon. I think our Zernikes are actually okay. Here's a short function to print out a human-readable polynomial from the _coef_array_xy (which derives from _xy_contribution) and some other details for a given Zernike polynomial:
def Zinfo(j, flipxy=False):
print(f"j = {j}")
Z = galsim.zernike.Zernike([0]*j+[1])
n, m = galsim.zernike.noll_to_zern(j)
print(f"n, m = {n}, {m}")
var = (1 if m==0 else 2)*(n+1)
arr = Z._coef_array_xy/np.sqrt(var)
print("\n_coef_array_xy")
print(arr)
first = True
out = f"\nsqrt({var}) * ("
for (i, k), val in np.ndenumerate(arr):
if val != 0:
if not first:
out += " + "
first = False
ival = int(np.round(val))
if ival != 1:
out += str(int(np.round(val)))
if flipxy:
i, k = k, i
if i >= 1:
out += "x"
if i >= 2:
out += "^"+str(i)
if k >= 1:
out += "y"
if k >= 2:
out += "^"+str(k)
out += ")"
# Clean up some special cases
out = out.replace("(-1x", "(-x")
out = out.replace("(-1y", "(-y")
out = out.replace("+ -", "- ")
if out == "sqrt(1) * ()":
out = "sqrt(1) * (1)"
print(out)
print(f"Z_{j}(0, 1) = sqrt({var}) * {int(np.round(Z(0, 1)/np.sqrt(var)))}")
print(f"Z_{j}(1, 0) = sqrt({var}) * {int(np.round(Z(1, 0)/np.sqrt(var)))}")
print(f"Z_{j}(1, 1) = sqrt({var}) * {int(np.round(Z(1, 1)/np.sqrt(var)))}")A few points of comparison: Zemax says Z_2 is 4^(1/2) (p) * COS (A). So assuming that p is sqrt(x^2 + y^2) and A is arctan2(y, x), this is Z_2 = 2 x. The above function, when called with j=2, prints out:
j = 2
n, m = 1, 1
_coef_array_xy
[[0. 0.]
[1. 0.]]
sqrt(4) * (x)
Z_2(0, 1) = sqrt(4) * 0
Z_2(1, 0) = sqrt(4) * 1
Z_2(1, 1) = sqrt(4) * 1I'll note that there does seem to be a bit of disagreement in the wild on whether A should be arctan2(y, x) or arctan2(x, y). E.g., slides 6-8 here. The flipxy param above will switch the human readable output from one convention to the other. The paper here prints out a number of Zernikes in the cartesian basis, but adopts this alternate definition for A. Using that and the Zinfo fn, I've verified we're generating the right Zernike functions up to j=36.
Along the way I also checked that horner2d is behaving like I think it should; namely that arr[i, j] is the coefficient of x^i y^j.
So anyway, I think our Zernikes are okay.
aaronroodman
rmjarvis
cwwalter
mperrin
Aaron Roodman has warned us for years that our
OpticalPSFclass doesn't match what Zemax predicts. This issue is to take a deep dive and figure out exactly what coordinate system conventions we and Zemax are using, and then either document them, or if merited, adopt Zemax's conventions. Potential points of confusion include: