Bug in make_hexagonal_grid()

[matdander]

I am working on building up a simulator for the AO systems at the CHARA Array. Our WFS uses a hexagonal lenslet array, so I had to build one from a hexagonal grid using the make_hexagonal_grid() function.

Our lenslet arrays have a pitch of 300 microns. So I converted this to circum_diameter as noted by the documentation. After a confusing afternoon, I found that this results in an incorrect hexagonal pattern. If I directly give the lens pitch (2*apothem) I get the expected pattern. At first I thought that there was a vocabulary error in the documentation, but I see in the code that the function does use the circum_diameter. The culprit in the code doesn't jump out at me, I'll keep poking at it to see if I can find the problem.

The following code exemplifies the issue. I compute the separation of two grid points, once using a vertical pair and again using a diagonal pair. Both result in the same value on the same grid, as expected. However, the result using the circum_diameter does not produce the expected lenslet pitch, but using lenslet_pitch does.

from hcipy import *

nrings=3; #Number of rings in the hexagonal wfs grid.
lenslet_pitch=300.0e-6;#300 micron center-to-center
circum_diameter=2*lenslet_pitch/np.sqrt(3); #Long diagonal of the hexagon. Different than pitch. pitch = 2 * apothem

print("Make grid using circum_diameter: ")
wfs_grid=make_hexagonal_grid(circum_diameter,nrings,pointy_top=False);

i=9
x1=wfs_grid.points[i,0]
x2=wfs_grid.points[i+1,0]
y1=wfs_grid.points[i,1]
y2=wfs_grid.points[i+1,1]
print("Lenslet Separation: ",np.around(np.sqrt((y2-y1)**2+(x2-x1)**2),6))

i=13
x1=wfs_grid.points[i,0]
x2=wfs_grid.points[i+1,0]
y1=wfs_grid.points[i,1]
y2=wfs_grid.points[i+1,1]
print("Lenslet Separation: ",np.around(np.sqrt((y2-y1)**2+(x2-x1)**2),6))

print("Make grid using lenslet_pitch: ")
wfs_grid=make_hexagonal_grid(lenslet_pitch,nrings,pointy_top=False);

i=9
x1=wfs_grid.points[i,0]
x2=wfs_grid.points[i+1,0]
y1=wfs_grid.points[i,1]
y2=wfs_grid.points[i+1,1]
print("Lenslet Separation: ",np.around(np.sqrt((y2-y1)**2+(x2-x1)**2),6))

i=13
x1=wfs_grid.points[i,0]
x2=wfs_grid.points[i+1,0]
y1=wfs_grid.points[i,1]
y2=wfs_grid.points[i+1,1]
print("Lenslet Separation: ",np.around(np.sqrt((y2-y1)**2+(x2-x1)**2),6))

Results in: Make grid using circum_diameter: Lenslet Separation: 0.000346 Lenslet Separation: 0.000346 Make grid using lenslet_pitch: Lenslet Separation: 0.0003 Lenslet Separation: 0.0003

Feel free to point out if I've made a silly mistake!

[syhaffert]

This is an actual bug. I encountered it a couple times before but just assumed it as a quirk of hexagons and never bothered to really check the source code. The mistake is found here: https://github.com/ehpor/hcipy/blob/master/hcipy/field/util.py#L213

apothem = circum_diameter * np.sqrt(3) / 4

This part of the code transforms the circum_diameter parameter in the wrong way. It should actually be: apothem = circum_diameter 2 / np.sqrt(3) / 2 The first 2 / np.sqrt(3) factor is the transformation from long side to short side and the additional / 2 is to go from diameter to radius.

For now I suggest to use the lenslet pitch, as you found out. To fix this part of the code will require us to update a couple other things too. It will probably also create some effects for #155 .

[matdander]

I'm perfectly happy going with the lenslet pitch. I just wanted to make sure you all were aware as it might save someone a headache in the future.

Thanks.

[syhaffert]

I created a pull request to fix this issue in #168

[syhaffert]

I had an offline discussion with @ehpor about this issue. The main point here is the difference between the definition of a hexagonal grid and the grid on which hexagonal mirrors are defined. The hexagonal grid is defined in such a way the the points lie on the vertices of the circumscribed hexagon. While for a segmented aperture like Keck, the hexagonal grid has a size such that each point lies in the center of each hexagon. HCIpy uses the first definition as it is the actual definition of a hexagonal grid. Therefore, there is no bug in this code. It is a documentation issue.