Grid of spheres

The pair of metal spheres from the chain example forms the unit cell of a square array with a period of 300 nm in the xy plane. A plane wave hits the array at normal incidence, traveling along z with its electric field along -x. The field at a point outside the spheres comes from expanding the waves of all cells in regular dipole waves about that point. The map shows the real part of E_z; at x = ±75 nm it vanishes by symmetry, and the printed values of about 1e-17 are rounding errors.

treams-rstreams 0.4.5

"""A square array of sphere pairs in the xy plane: the near field of the array."""

import numpy as np

import treams_rs as tr

# Enlarge for a full map: treams uses a 31 x 31 grid.
GRID_POINTS = 5

k0 = 2 * np.pi / 700
material = tr.Material(-16.5 + 1j)
lmax = 3
radii = [75, 75]
positions = [[-75, 0, 30], [75, 0, -30]]
lattice = tr.Lattice.square(300)
kpar = [0, 0]

spheres = [
    tr.sphere_tmatrix(k0=k0, lmax=lmax, radius=r, material=material) for r in radii
]
# A two-dimensional lattice is periodic in the xy plane.
array = tr.solve_periodic(
    tr.Cluster(spheres, positions=positions), lattice=lattice, kpar=kpar
)

inc = tr.plane_wave(direction=[0, 0, 1], pol=[-1, 0, 0], k0=k0)
sca = array.scatter(inc)

x = np.linspace(-150, 150, GRID_POINTS)
z = np.linspace(-150, 150, GRID_POINTS)
xx, zz = np.meshgrid(x, z, indexing="ij")
grid = np.stack((xx, np.zeros_like(xx), zz), axis=-1)

# Points inside a sphere of the reference cell get NaN, as in treams.
distances = np.linalg.norm(grid[..., None, :] - np.array(positions), axis=-1)
valid = np.all(distances > np.array(radii), axis=-1)

# The waves of all cells, expanded in regular dipole waves about each point,
# give the scattered field at that point.
ez = np.full_like(xx, np.nan)
ez[valid] = [
    np.real(
        inc.efield(r)[2]
        + sca.in_basis(
            tr.SphericalBasis.default(1, positions=[r]), kind="regular"
        ).efield(r)[2]
    )
    for r in grid[valid]
]

if __name__ == "__main__":
    print("Real part of E_z at y = 0 (rows: x, columns: z)")
    for row in ez:
        print("".join(f"{v:>13.6g}" for v in row))

# Source: https://github.com/tfp-photonics/treams/blob/1f5d0d6ebb007288f28bc9e16f6d266e8b55dc39/docs/examples/grid.py
# Adapted: reduced sizes, plotting removed.
# MIT license, copyright Dominik Beutel; see LICENSE.treams.
import numpy as np
import treams

k0 = 2 * np.pi / 700
materials = [treams.Material(-16.5 + 1j), treams.Material()]
lmax = 3
radii = [75, 75]
positions = [[-75, 0, 30], [75, 0, -30]]
lattice = treams.Lattice.square(300)
kpar = [0, 0]

spheres = [treams.TMatrix.sphere(lmax, k0, r, materials) for r in radii]
array = treams.TMatrix.cluster(spheres, positions).latticeinteraction.solve(
    lattice, kpar
)

inc = treams.plane_wave([0, 0, k0], [-1, 0, 0], k0=array.k0, material=array.material)
sca = array @ inc.expand(array.basis)

x = np.linspace(-150, 150, 5)
y = 0
z = np.linspace(-150, 150, 5)
xx, zz = np.meshgrid(x, z, indexing="ij")
yy = np.full_like(xx, y)
grid = np.stack((xx, yy, zz), axis=-1)

ez = np.zeros_like(xx)
valid = array.valid_points(grid, radii)
vals = []
for r in grid[valid]:
    swb = treams.SphericalWaveBasis.default(1, positions=[r])
    field = sca.expandlattice(basis=swb).efield(r)
    vals.append(np.real(inc.efield(r)[2] + field[2]))
ez[valid] = vals
ez[~valid] = np.nan

Output

Real part of E_z at y = 0 (rows: x, columns: z)
    0.0413553     0.506245     -2.06494     0.497878    0.0389712
  1.35913e-17  4.09508e-18          nan          nan  2.81282e-18
   -0.0413553    -0.506245      2.06494    -0.497878   -0.0389712
 -4.45271e-19          nan          nan  3.20194e-18 -5.06185e-18
    0.0413553     0.506245     -2.06494     0.497878    0.0389712

Differences from the treams version

• solve_periodic(Cluster(...), lattice=Lattice.square(300), kpar=[0, 0]) returns a periodic response that keeps its lattice and Bloch vector; treams writes TMatrix.cluster(...).latticeinteraction.solve(lattice, kpar).
• in_basis(basis, kind="regular") expands the waves of all cells about a point; treams writes sca.expandlattice(basis=swb).
• The periodic response has no valid_points, so the script finds the points outside the spheres of the unit cell with NumPy.