Skip to content

Chain and cylinder

A chain of metal spheres along z and an infinite cylinder stand side by side and couple. The chain has spheres of radius 65 nm every 300 nm, with its axis at (x, y) = (-40, -50) nm; the cylinder has radius 55 nm, with its axis at (40, 50) nm. Both have permittivity -16.5 + 1i. As in the chain example, to_cylindrical describes the chain by cylindrical waves. The cylinder uses the same seven axial wavenumbers, because waves with different kz do not couple. Cluster(...).solve() then gives the T-matrix of the pair. The map shows the real part of E_z for a plane wave with kz = 0.005 nm⁻¹; the row x = 0 lies inside the chain's cylinder.

"""A sphere chain next to a cylinder: the near field of the coupled pair."""

import numpy as np

import treams_rs as tr

# Enlarge for a full map: treams uses a 61 x 31 grid.
GRID_POINTS = (7, 5)

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 = tr.Lattice(300)
kz = 0.005

# A chain of spheres along z, expressed in cylindrical waves.
sphere = tr.sphere_tmatrix(k0=k0, lmax=lmax, radius=radii[0], material=material)
chain = tr.solve_periodic(sphere, lattice=lattice, kpar=kz)
bmax = 3.1 * lattice.reciprocal
basis = tr.CylindricalBasis.diffr_orders(kz, mmax, lattice, bmax)
chain_tmc = chain.to_cylindrical(basis)

# The cylinder needs the same axial wavenumbers as the chain to couple to it.
cylinder = tr.cylinder_tmatrix(
    k0=k0, kz=np.unique(basis.kz), mmax=mmax, radius=radii[1], material=material
)
cluster = tr.Cluster([chain_tmc, cylinder], positions=positions).solve()

inc = tr.plane_wave(
    direction=[np.sqrt(k0**2 - kz**2), 0, kz],
    pol=[np.sqrt(0.5), np.sqrt(0.5)],
    k0=k0,
)
sca = cluster.scatter(inc)

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

# Points inside either cylinder get NaN: the expansion holds only outside them.
ez = np.full_like(xx, np.nan, dtype=complex)
valid = cluster.valid_points(grid, radii)
ez[valid] = (inc.efield(grid[valid]) + sca.efield(grid[valid]))[..., 2]

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

# Source: https://github.com/tfp-photonics/treams/blob/1f5d0d6ebb007288f28bc9e16f6d266e8b55dc39/docs/examples/cluster_tmatrixc.py
# Adapted: reduced sizes, NaN inside the cylinders, plotting and the copies of
# neighboring cells 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 = treams.Lattice(300)
kz = 0.005

sphere = treams.TMatrix.sphere(lmax, k0, radii[0], materials)
chain = sphere.latticeinteraction.solve(lattice, kz)
bmax = 3.1 * lattice.reciprocal
cwb = treams.CylindricalWaveBasis.diffr_orders(kz, mmax, lattice, bmax)
chain_tmc = treams.TMatrixC.from_array(chain, cwb)

cylinder = treams.TMatrixC.cylinder(np.unique(cwb.kz), mmax, k0, radii[1], materials)

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

x = np.linspace(-300, 300, 7)
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, complex)
valid = cluster.valid_points(grid, radii)
ez[valid] = (inc.efield(grid[valid]) + sca.efield(grid[valid]))[..., 2]
ez[~valid] = np.nan

Output

Real part of E_z at y = 0 (rows: x, columns: z)
    -0.596941    -0.510939    -0.343232    -0.113167     0.127571
    -0.327837   -0.0471874     0.263634     0.643765     0.884022
     0.146854     0.547241     0.153631      1.48105      1.39048
          nan          nan          nan          nan          nan
     0.405381     0.379218     0.270026     0.148964    0.0154212
     0.381432     0.210826   0.00837416    -0.189323    -0.354667
    0.0227267    -0.200939    -0.396959    -0.536238     -0.59962

Differences from the treams version

  • Cluster([chain_tmc, cylinder], positions=positions) holds the uncoupled T-matrices and solve() couples them; treams writes TMatrixC.cluster(...).interaction.solve().
  • to_cylindrical(basis) converts the solved chain; treams writes TMatrixC.from_array(chain, cwb).
  • cylinder_tmatrix(kz=...) takes the axial wavenumbers as a keyword; treams takes them as the first argument of TMatrixC.cylinder.
  • Points inside a cylinder get NaN in both scripts. The treams version leaves them at zero and repeats the map for the neighboring cells; both scripts here stop at one cell.