Skip to content

Single cylinder

An infinite cylinder along z with permittivity 16 + 0.5i and radius 75 nm scatters light that travels in the xy plane, so the axial wavenumber kz is 0. The spectrum shows the extinction and scattering efficiencies: the cross widths averaged over all directions in the xy plane and both polarizations, divided by the diameter. A cross width is the cross section per unit length of the cylinder. Keeping only the modes of order m = 0 shows which resonances the higher orders cause. A plane wave at 1000 THz then gives the intensity of the total field around the cylinder.

"""Cross widths of an infinite cylinder and the intensity around it."""

import numpy as np

import treams_rs as tr

# Enlarge these for full sweeps: treams uses 200 frequencies and a 101 x 101 grid.
FREQUENCIES = 6
GRID_POINTS = 7

k0s = 2 * np.pi * np.linspace(1 / 6000, 1 / 300, FREQUENCIES)
material = tr.Material(16 + 0.5j)
mmax = 4
radius = 75
# Light that travels in the xy plane has the axial wavenumber kz = 0.
kzs = [0]
cylinders = [
    tr.cylinder_tmatrix(k0=k0, kz=kzs, mmax=mmax, radius=radius, material=material)
    for k0 in k0s
]

# Efficiencies are cross widths divided by the diameter.
averages = [tm.average_cross_widths for tm in cylinders]
xw_sca = np.array([avg.scattering for avg in averages]) / (2 * radius)
xw_ext = np.array([avg.extinction for avg in averages]) / (2 * radius)

# Keep only the modes of order m = 0 of the same T-matrices.
order_zero = tr.CylindricalBasis.default(kzs, 0)
averages_mmax0 = [tm.select(order_zero).average_cross_widths for tm in cylinders]
xw_sca_mmax0 = np.array([avg.scattering for avg in averages_mmax0]) / (2 * radius)
xw_ext_mmax0 = np.array([avg.extinction for avg in averages_mmax0]) / (2 * radius)

# Total field of a plane wave at the highest frequency, in the plane z = 0.
tm = cylinders[-1]
inc = tr.plane_wave(direction=[1, 0, 0], pol="positive_helicity", k0=tm.k0)
sca = tm.scatter(inc)

x = np.linspace(-100, 100, GRID_POINTS)
y = np.linspace(-100, 100, GRID_POINTS)
xx, yy = np.meshgrid(x, y, indexing="ij")
grid = np.stack((xx, yy, np.zeros_like(xx)), axis=-1)

# Points inside the cylinder get NaN: the expansion holds only outside it.
intensity = np.full_like(xx, np.nan)
valid = tm.valid_points(grid, [radius])
total = inc.efield(grid[valid]) + sca.efield(grid[valid])
intensity[valid] = 0.5 * np.sum(np.abs(total) ** 2, axis=-1)

if __name__ == "__main__":
    print("Efficiencies (cross width / diameter)")
    headers = ("frequency/THz", "extinction", "scattering", "ext. m=0", "sca. m=0")
    print("".join(f"{h:>15}" for h in headers))
    frequencies = k0s * 299792.458 / (2 * np.pi)
    columns = (frequencies, xw_ext, xw_sca, xw_ext_mmax0, xw_sca_mmax0)
    for row in zip(*columns, strict=True):
        print("".join(f"{v:>15.6g}" for v in row))
    print("\nIntensity at z = 0 (rows: x, columns: y)")
    for row in intensity:
        print("".join(f"{v:>11.6g}" for v in row))

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

k0s = 2 * np.pi * np.linspace(1 / 6000, 1 / 300, 6)
materials = [treams.Material(16 + 0.5j), treams.Material()]
mmax = 4
radius = 75
kzs = [0]
cylinders = [treams.TMatrixC.cylinder(kzs, mmax, k0, radius, materials) for k0 in k0s]

xw_sca = np.array([tm.xw_sca_avg for tm in cylinders]) / (2 * radius)
xw_ext = np.array([tm.xw_ext_avg for tm in cylinders]) / (2 * radius)

cwb_mmax0 = treams.CylindricalWaveBasis.default(kzs, 0)
cylinders_mmax0 = [tm[cwb_mmax0] for tm in cylinders]
xw_sca_mmax0 = np.array([tm.xw_sca_avg for tm in cylinders_mmax0]) / (2 * radius)
xw_ext_mmax0 = np.array([tm.xw_ext_avg for tm in cylinders_mmax0]) / (2 * radius)

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

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

intensity = np.zeros_like(xx)
valid = tm.valid_points(grid, [radius])
intensity[~valid] = np.nan
intensity[valid] = 0.5 * np.sum(
    np.abs(inc.efield(grid[valid]) + sca.efield(grid[valid])) ** 2, -1
)

Output

Efficiencies (cross width / diameter)
  frequency/THz     extinction     scattering       ext. m=0       sca. m=0
        49.9654       0.129913      0.0888491       0.128959      0.0883775
        239.834        2.55707         2.4827        2.47201        2.41391
        429.703        3.37011        3.14936         1.8746        1.79484
        619.571        1.65808        1.10053        0.46281       0.300058
         809.44        1.48166       0.903129       0.680195       0.598727
        999.308        2.68786        2.33452        1.02044       0.971437

Intensity at z = 0 (rows: x, columns: y)
   0.900581   0.792748   0.609413   0.494534   0.609413   0.792748   0.900581
    0.57268   0.486247        nan        nan        nan   0.486247    0.57268
   0.302176        nan        nan        nan        nan        nan   0.302176
   0.316976        nan        nan        nan        nan        nan   0.316976
   0.418043        nan        nan        nan        nan        nan   0.418043
   0.287753   0.250069        nan        nan        nan   0.250069   0.287753
   0.174415    0.15122   0.177147   0.200564   0.177147    0.15122   0.174415

Differences from the treams version

  • cylinder_tmatrix takes keywords and the cylinder's material; the exterior medium is vacuum unless given. treams takes one list of materials, ending with the exterior.
  • average_cross_widths returns named scattering and extinction values; treams reads xw_sca_avg and xw_ext_avg.
  • select(CylindricalBasis.default(kzs, 0)) keeps the modes of order 0; treams indexes with tm[cwb_mmax0].
  • scatter returns the scattered wave; treams writes tm @ inc.expand(tm.basis).
  • The intensity map starts as NaN, so points inside the cylinder stay NaN.