.. image:: https://mybinder.org/badge_logo.svg :target: https://mybinder.org/v2/gh/ehpor/hcipy/HEAD?labpath=doc/tutorial_notebooks/CoreConceptsGridsFields/CoreConceptsGridsFields.ipynb :alt: Open in Binder Core Concepts: Grids and Fields =============================== In this notebook we’ll explore grids and fields, the two fundamental data structures in HCIPy. Grids define the spatial coordinates where we evaluate optical elements, and fields store physical quantities like electric field amplitude or phase at those points. .. code:: ipython3 import hcipy import numpy as np import matplotlib.pyplot as plt We’ll start with a uniform Cartesian grid, the most common grid type for general-purpose simulations. A ``CartesianGrid`` stores points on a rectangular lattice, making it ideal for Fourier transforms and standard image processing. .. code:: ipython3 # Uniform Cartesian grid grid = hcipy.make_uniform_grid([32, 32], [64, 32]) print(f"Grid dimensions: {grid.dims}") print(f"Grid spacing: {grid.delta}") print(f"Total extent: {grid.delta[0] * grid.dims[0]:.4f} x {grid.delta[1] * grid.dims[1]:.4f}") plt.figure(figsize=(10, 6)) plt.plot(grid.x, grid.y, '+', markersize=8, alpha=0.6) plt.xlabel('x') plt.ylabel('y') plt.title('Cartesian Grid Sampling Pattern') plt.grid(True, alpha=0.3) plt.gca().set_aspect('equal') plt.show() .. parsed-literal:: Grid dimensions: (32, 32) Grid spacing: [2. 1.] Total extent: 64.0000 x 32.0000 .. image:: output_3_1.png The ``.dims`` property tells us the number of points in each direction, ``.delta`` gives the spacing between adjacent points, and the total extent is simply their product. Notice that our grid has the same number of points in each direction (32 by 32) but the physical extent is wider in x (64 units) than in y (32 units) because the x-spacing is larger. These properties are the building blocks for defining the spatial scale of any simulation. Now let’s try a polar grid. Polar grids are useful for problems with circular symmetry — we construct one using ``SeparatedCoords`` and ``PolarGrid``. .. code:: ipython3 # Polar grid with logarithmic radial spacing r = np.logspace(-1, 1, 11) theta = np.linspace(0, 2 * np.pi, 10, endpoint=False) coords = hcipy.SeparatedCoords((r, theta)) polar_grid = hcipy.PolarGrid(coords) # Convert to Cartesian for visualization cart_grid = polar_grid.as_('cartesian') plt.figure(figsize=(8, 8)) plt.plot(cart_grid.x, cart_grid.y, '+', markersize=10) plt.axis('equal') plt.xlabel('x') plt.ylabel('y') plt.title('Polar Grid (Logarithmic Radial Spacing)') plt.grid(True, alpha=0.3) plt.show() print(f"Polar grid size: {polar_grid.size}") print(f"Radial range: {polar_grid.r.min():.2f} to {polar_grid.r.max():.2f}") .. image:: output_6_0.png .. parsed-literal:: Polar grid size: 110 Radial range: 0.10 to 10.00 With logarithmic radial spacing the points cluster near the center, which is useful when the optical field varies fastest near the origin. We converted the polar grid to Cartesian coordinates with ``.as_('cartesian')``. Grids are just coordinate containers, though — they don’t hold any physical values on their own. To attach data to grid points we need a field. A ``Field`` is a NumPy-like array that carries a reference to the grid it’s sampled on. This association means operations on the field preserve the grid, and visualization tools like ``imshow_field`` automatically know how to display the data. Let’s create a real-valued field with a Gaussian intensity pattern. .. code:: ipython3 # Pupil grid for field examples pupil_grid = hcipy.make_pupil_grid(128, diameter=1.5) # Real field: Gaussian intensity pattern x = pupil_grid.x y = pupil_grid.y gaussian_values = np.exp(-(x**2 + y**2) / 0.3) field = hcipy.Field(gaussian_values, pupil_grid) The ``Field`` constructor takes an array of values and the grid they correspond to — now every operation on ``field`` knows about the ``pupil_grid``. We can visualize the full 2D field and also use ``.shaped`` to reshape the flat data into the grid’s natural multi-dimensional layout for slicing. .. code:: ipython3 plt.figure(figsize=(10, 4)) plt.subplot(121) hcipy.imshow_field(field, cmap='inferno') plt.colorbar(label='Intensity') plt.title('Gaussian Field (2D)') plt.subplot(122) center_slice = field.shaped[field.shaped.shape[0]//2, :] plt.plot(pupil_grid.x.reshape(pupil_grid.dims)[0, :], center_slice) plt.xlabel('x') plt.ylabel('Intensity') plt.title('Cross-Section') plt.grid(True, alpha=0.3) plt.tight_layout() plt.show() print(f"Field size: {field.size} points") .. image:: output_11_0.png .. parsed-literal:: Field size: 16384 points The ``.shaped`` attribute lets us index into the field as if it were a 2D array — here we grab the center row to plot a cross-section of the Gaussian. Fields also support standard arithmetic: addition, scalar multiplication, and element-wise operations all work as expected, preserving the grid link. We can see the Gaussian profile is smooth and symmetric, as expected. Now let’s move to complex-valued fields, which are essential in wave optics. The electric field is typically written as E = A exp(i phi), where the amplitude A and phase phi encode the wavefront. HCIPy’s ``imshow_field`` can display both at once using a two-color visualization: brightness for amplitude and hue for phase. .. code:: ipython3 # Complex field: Gaussian amplitude with linear phase tilt grid = hcipy.make_uniform_grid([128, 128], [5, 5]) x, y = grid.coords r = grid.as_('polar').r z = np.exp(-r**2) * np.exp(-1j*(5 * x + 2 * y)) z = hcipy.Field(z, grid) This field combines a Gaussian amplitude profile with a linear phase tilt — the tilt comes from the exp(-i(5x + 2y)) factor and causes the wavefront to propagate at an angle. Let’s extract the amplitude and phase separately to see each component. .. code:: ipython3 plt.figure(figsize=(12, 5)) plt.subplot(131) hcipy.imshow_field(z) plt.title('Complex Field (Amplitude + Phase)') plt.subplot(132) hcipy.imshow_field(np.abs(z), cmap='inferno') plt.title('Amplitude Only') plt.subplot(133) hcipy.imshow_field(np.angle(z), vmin=-np.pi, vmax=np.pi, cmap='hsv') plt.title('Phase Only') plt.tight_layout() plt.show() print(f"Field is complex: {np.iscomplexobj(z)}") print(f"Phase range: {np.min(np.angle(z)):.2f} to {np.max(np.angle(z)):.2f} radians") .. image:: output_16_0.png .. parsed-literal:: Field is complex: True Phase range: -3.14 to 3.14 radians The left panel shows the combined view: brightness fades from center to edge (the Gaussian amplitude) while colors cycle through the rainbow (the linear phase tilt). The phase-only panel confirms the tilt is a smooth ramp, wrapping at +/- pi which creates the color discontinuities. Notice that ``np.abs`` and ``np.angle`` work directly on a ``Field``, returning new fields on the same grid. An important pattern in HCIPy is separating an optical element’s definition from how it’s sampled on a grid. Aperture generators are functions that take a grid and return a field — we can define the geometry once and evaluate it on any grid. .. code:: ipython3 # Aperture generators circular_aperture = hcipy.make_circular_aperture(diameter=1.0) rectangular_aperture = hcipy.make_rectangular_aperture([1.0, 0.5]) hexagonal_aperture = hcipy.make_hexagonal_aperture(1.0) We created three aperture generators — circular, rectangular, and hexagonal — and can now evaluate all three on the same pupil grid to compare them side-by-side. .. code:: ipython3 # Evaluate on the pupil grid circ_pupil = circular_aperture(pupil_grid) rect_pupil = rectangular_aperture(pupil_grid) hex_pupil = hexagonal_aperture(pupil_grid) fig, axes = plt.subplots(1, 3, figsize=(15, 4)) hcipy.imshow_field(circ_pupil, ax=axes[0], cmap='gray') axes[0].set_title('Circular Aperture') axes[0].axis('off') hcipy.imshow_field(rect_pupil, ax=axes[1], cmap='gray') axes[1].set_title('Rectangular Aperture') axes[1].axis('off') hcipy.imshow_field(hex_pupil, ax=axes[2], cmap='gray') axes[2].set_title('Hexagonal Aperture') axes[2].axis('off') plt.tight_layout() plt.show() .. image:: output_20_0.png Each generator correctly samples its geometry onto the grid at whatever resolution we provide. Even with a moderate pupil grid the three shapes are clearly distinguishable. When apertures have sharp edges — especially thin features like telescope spider arms — direct sampling can produce aliasing artifacts. The result is jagged or broken edges that don’t represent the true optical element. Supersampling fixes this by evaluating the aperture at higher resolution and averaging down, effectively integrating each pixel over the fine structure. .. code:: ipython3 # Magellan telescope aperture magellan = hcipy.make_magellan_aperture(normalized=True) # Direct sampling pupil_direct = magellan(pupil_grid) The Magellan telescope aperture has thin spider arms that are especially prone to aliasing. We can see what the direct evaluation looks like — then compare with a supersampled version. .. code:: ipython3 # 8x supersampling pupil_supersampled = hcipy.evaluate_supersampled(magellan, pupil_grid, 8) ``evaluate_supersampled`` takes three arguments: the aperture function, the target grid, and a supersampling factor. Here we’re evaluating at 8x the native resolution and averaging back down to the original grid. Now we can compare the two approaches directly. .. code:: ipython3 fig, axes = plt.subplots(1, 2, figsize=(14, 6)) hcipy.imshow_field(pupil_direct, ax=axes[0], cmap='gray') axes[0].set_title('Direct Sampling (Aliasing on spider arms)') hcipy.imshow_field(pupil_supersampled, ax=axes[1], cmap='gray') axes[1].set_title('8x Supersampling (Smooth edges)') plt.tight_layout() plt.show() # Zoomed comparison print("Zoomed view of spider arm region:") fig, axes = plt.subplots(1, 2, figsize=(14, 6)) hcipy.imshow_field(pupil_direct, ax=axes[0], cmap='gray') axes[0].set_title('Direct Sampling (zoomed)') axes[0].set_xlim(-0.05, 0.25) axes[0].set_ylim(-0.05, 0.25) hcipy.imshow_field(pupil_supersampled, ax=axes[1], cmap='gray') axes[1].set_title('8x Supersampling (zoomed)') axes[1].set_xlim(-0.05, 0.25) axes[1].set_ylim(-0.05, 0.25) plt.tight_layout() plt.show() .. image:: output_26_0.png .. parsed-literal:: Zoomed view of spider arm region: .. image:: output_26_2.png The difference is stark in the zoomed view: the direct-sampled spider arm is broken into disconnected pixels, while the supersampled version is smooth and continuous. For quantitative work — like computing contrast curves or coupling efficiencies — this accuracy matters. Finally, let’s look at grid weights. Each grid point has an associated weight representing the area it covers. For uniform Cartesian grids all weights are identical; for non-uniform grids they vary. These weights enable numerical integration over the grid. .. code:: ipython3 # Grid weights print(f"Individual point weight: {pupil_grid.weights:.6f} square meters") print(f"Total grid area: {(pupil_grid.ones() * pupil_grid.weights).sum():.4f} square meters") print(f"Grid dimensions: {pupil_grid.dims}") # Compute aperture area via weighted integration aperture = hcipy.make_circular_aperture(diameter=1.0)(pupil_grid) computed_area = (aperture * pupil_grid.weights).sum() expected_area = np.pi * (0.5)**2 print(f"\nCircular Aperture Area Calculation:") print(f" Computed area: {computed_area:.6f}") print(f" Expected area (πr²): {expected_area:.6f}") print(f" Relative error: {abs(computed_area - expected_area)/expected_area * 100:.3f}%") .. parsed-literal:: Individual point weight: 0.000137 square meters Total grid area: 2.2500 square meters Grid dimensions: (128, 128) Circular Aperture Area Calculation: Computed area: 0.785522 Expected area (πr²): 0.785398 Relative error: 0.016% The weighted sum reproduces the true area of a unit-diameter circle to within a fraction of a percent — even at this moderate resolution. Grid weights are used throughout HCIPy for everything from energy conservation checks to computing Zernike coefficient overlaps, so it’s worth knowing they’re there.