Skip to content

Grating

The chain of spheres and the cylinder from the chain and cylinder example repeat every 300 nm along x and form a grating. Here the chain uses the three axial wavenumbers kz + 2πn / 300 nm with |n| ≤ 1. A plane wave with kz = 0.005 nm⁻¹ travels almost along y. The field at a point outside the cylinders comes from expanding the waves of all cells in regular cylindrical waves about that point. The map shows the real part of E_x in the plane x = 0; the rows y = ±50 nm lie inside the cylinders.

"""A grating of sphere chains and cylinders: the near field in one period."""

import numpy as np

import treams_rs as tr

# Enlarge for a full map: treams uses a 16 x 16 grid.
GRID_POINTS = 4

k0 = 2 * np.pi / 700
material = tr.Material(-16.5 + 1j)
lmax = mmax = 3
radii = [65, 55]
positions = [[-40, -50, 0], [40, 50, 0]]
lattice_z = tr.Lattice(300)
lattice_x = tr.Lattice(300, "x")
kz = 0.005

# A chain of spheres along z in cylindrical waves, as in the cylindrical
# cluster example, with the axial orders |n| <= 1.1.
sphere = tr.sphere_tmatrix(k0=k0, lmax=lmax, radius=radii[0], material=material)
chain = tr.solve_periodic(sphere, lattice=lattice_z, kpar=kz)
bmax = 1.1 * lattice_z.reciprocal
basis = tr.CylindricalBasis.diffr_orders(kz, mmax, lattice_z, bmax)
kzs = np.unique(basis.kz)
chain_tmc = chain.to_cylindrical(basis)
cylinder = tr.cylinder_tmatrix(
    k0=k0, kz=kzs, mmax=mmax, radius=radii[1], material=material
)

# The chain and the cylinder repeat every 300 nm along x.
grating = tr.solve_periodic(
    tr.Cluster([chain_tmc, cylinder], positions=positions), lattice=lattice_x, kpar=0
)
inc = tr.plane_wave(
    direction=[0, np.sqrt(k0**2 - kz**2), kz],
    pol=[np.sqrt(0.5), -np.sqrt(0.5)],
    k0=k0,
)
sca = grating.scatter(inc)

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

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

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

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

# Source: https://github.com/tfp-photonics/treams/blob/1f5d0d6ebb007288f28bc9e16f6d266e8b55dc39/docs/examples/grating_tmatrixc.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 = mmax = 3
radii = [65, 55]
positions = [[-40, -50, 0], [40, 50, 0]]
lattice_z = treams.Lattice(300)
lattice_x = treams.Lattice(300, "x")
kz = 0.005

sphere = treams.TMatrix.sphere(lmax, k0, radii[0], materials)
chain = sphere.latticeinteraction.solve(lattice_z, kz)
bmax = 1.1 * lattice_z.reciprocal
cwb = treams.CylindricalWaveBasis.diffr_orders(kz, mmax, lattice_z, bmax)
kzs = np.unique(cwb.kz)
chain_tmc = treams.TMatrixC.from_array(chain, cwb)
cylinder = treams.TMatrixC.cylinder(kzs, mmax, k0, radii[1], materials)

cluster = treams.TMatrixC.cluster(
    [chain_tmc, cylinder], positions
).latticeinteraction.solve(lattice_x, 0)
inc = treams.plane_wave(
    [0, np.sqrt(k0**2 - kz**2), kz],
    [np.sqrt(0.5), -np.sqrt(0.5)],
    k0=chain.k0,
    material=chain.material,
)
sca = cluster @ inc.expand(cluster.basis)

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

ex = np.zeros_like(xx)
valid = cluster.valid_points(grid, radii)
vals = []
for r in grid[valid]:
    cwb = treams.CylindricalWaveBasis.default(kzs, 1, positions=[r])
    field = sca.expandlattice(basis=cwb).efield(r)
    vals.append(np.real(inc.efield(r)[0] + field[0]))
ex[valid] = vals
ex[~valid] = np.nan

Output

Real part of E_x at x = 0 (rows: y, columns: z)
    -0.705071    -0.944121    -0.988146    -0.761953
          nan          nan          nan          nan
          nan          nan          nan          nan
     0.698716     0.818856      0.71492     0.450867

Differences from the treams version

  • solve_periodic(Cluster(...), lattice=Lattice(300, "x"), kpar=0) returns a periodic response that keeps its lattice and Bloch vector; treams writes TMatrixC.cluster(...).latticeinteraction.solve(lattice_x, 0).
  • in_basis(basis, kind="regular") expands the waves of all cells about a point; treams writes sca.expandlattice(basis=cwb).
  • The periodic response has no valid_points, so the script finds the points outside the cylinders of the unit cell with NumPy.
  • The results agree with treams to a relative 1.3e-8. treams sums the lattice of cylinders with an automatic Ewald split (the parameter that divides the sum into a real-space and a reciprocal-space part) that leaves errors of about 5e-9. With a larger explicit split, treams agrees with treams-rs to 1e-12. One explicit split cannot serve this example, because its axial orders include both propagating and decaying waves; see differences.