Chain of spheres

Two metal spheres with permittivity -16.5 + 1i form the unit cell of a chain that repeats every 300 nm along z. A plane wave travels along x with its electric field along z, so the Bloch wave vector along the chain is zero. solve_periodic solves the interaction of all cells at once. The scattered 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.

treams-rstreams 0.4.5

"""A chain of sphere pairs, periodic along z: 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 = [[-30, 0, -75], [30, 0, 75]]
lattice = 300
kz = 0

spheres = [
    tr.sphere_tmatrix(k0=k0, lmax=lmax, radius=r, material=material) for r in radii
]
# A one-dimensional lattice is periodic along z.
chain = tr.solve_periodic(
    tr.Cluster(spheres, positions=positions), lattice=lattice, kpar=kz
)

inc = tr.plane_wave(direction=[1, 0, 0], pol=[0, 0, 1], k0=k0)
sca = chain.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/chain.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 = [[-30, 0, -75], [30, 0, 75]]
lattice = 300
kz = 0

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

inc = treams.plane_wave([k0, 0, 0], [0, 0, 1], k0=chain.k0, material=chain.material)
sca = chain @ inc.expand(chain.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 = chain.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.369221    -0.353395    -0.369221    -0.372175    -0.369221
    -0.068106          nan    -0.068106    -0.086169    -0.068106
     0.457191          nan     0.457191          nan     0.457191
   -0.0959318   -0.0939806   -0.0959318          nan   -0.0959318
    -0.369093    -0.372793    -0.369093    -0.349305    -0.369093

Differences from the treams version

• solve_periodic(Cluster(...), lattice=300, kpar=0) returns a periodic response that keeps its lattice and Bloch vector; treams writes TMatrix.cluster(...).latticeinteraction.solve(lattice, kz) and returns a T-matrix.
• scatter returns the waves of all cells. in_basis(basis, kind="regular") expands them 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, by the same rule as treams.