Skip to content

Single sphere

A sphere with permittivity 16 + 0.5i and radius 75 nm has electric and magnetic multipole resonances between 430 and 1000 THz. The spectrum shows the extinction and scattering efficiencies, the cross sections averaged over all directions and polarizations of the incident light and divided by the geometric cross section πr². Keeping only the dipole modes (l = 1) shows which resonances higher multipoles cause. A plane wave at 1000 THz then gives the intensity of the total field around the sphere.

"""Extinction and scattering spectrum of a sphere 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 / 700, 1 / 300, FREQUENCIES)
material = tr.Material(16 + 0.5j)
lmax = 4
radius = 75
spheres = [
    tr.sphere_tmatrix(k0=k0, lmax=lmax, radius=radius, material=material) for k0 in k0s
]

# Efficiencies are cross sections divided by the geometric cross section.
area = np.pi * radius**2
averages = [tm.average_cross_sections for tm in spheres]
xs_sca = np.array([avg.scattering for avg in averages]) / area
xs_ext = np.array([avg.extinction for avg in averages]) / area

# Keep only the dipole modes (l = 1) of the same T-matrices.
dipoles = tr.SphericalBasis.default(1)
averages_lmax1 = [tm.select(dipoles).average_cross_sections for tm in spheres]
xs_sca_lmax1 = np.array([avg.scattering for avg in averages_lmax1]) / area
xs_ext_lmax1 = np.array([avg.extinction for avg in averages_lmax1]) / area

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

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

# Points inside the sphere 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 section / geometric cross section)")
    headers = ("frequency/THz", "extinction", "scattering", "ext. l=1", "sca. l=1")
    print("".join(f"{h:>15}" for h in headers))
    frequencies = k0s * 299792.458 / (2 * np.pi)
    columns = (frequencies, xs_ext, xs_sca, xs_ext_lmax1, xs_sca_lmax1)
    for row in zip(*columns, strict=True):
        print("".join(f"{v:>15.6g}" for v in row))
    print("\nIntensity at y = 0 (rows: x, columns: z)")
    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/sphere.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 / 700, 1 / 300, 6)
materials = [treams.Material(16 + 0.5j), treams.Material()]
lmax = 4
radius = 75
spheres = [treams.TMatrix.sphere(lmax, k0, radius, materials) for k0 in k0s]

xs_sca = np.array([tm.xs_sca_avg for tm in spheres]) / (np.pi * radius**2)
xs_ext = np.array([tm.xs_ext_avg for tm in spheres]) / (np.pi * radius**2)

swb_lmax1 = treams.SphericalWaveBasis.default(1)
spheres_lmax1 = [tm[swb_lmax1] for tm in spheres]
xs_sca_lmax1 = np.array([tm.xs_sca_avg for tm in spheres_lmax1]) / (np.pi * radius**2)
xs_ext_lmax1 = np.array([tm.xs_ext_avg for tm in spheres_lmax1]) / (np.pi * radius**2)

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

x = np.linspace(-100, 100, 7)
y = 0
z = np.linspace(-100, 100, 7)
xx, zz = np.meshgrid(x, z, indexing="ij")
yy = np.full_like(xx, y)
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 section / geometric cross section)
  frequency/THz     extinction     scattering       ext. l=1       sca. l=1
        428.275        1.33296       0.987177        1.33118       0.986812
        542.482        3.51292        3.10639        3.50198        3.10339
        656.688        5.39143         4.3844        5.15451        4.32194
        770.895        1.02102        0.60478       0.766772       0.467907
        885.102        3.04969        1.59416        1.23458        1.02832
        999.308        3.23746        2.55864        3.02429        2.47346

Intensity at y = 0 (rows: x, columns: z)
   0.532304   0.332909   0.305677   0.437686   0.456637   0.408559   0.393593
   0.394684   0.209014        nan        nan        nan   0.684083   0.539167
   0.256366        nan        nan        nan        nan        nan   0.594802
   0.190624        nan        nan        nan        nan        nan    0.54269
   0.256366        nan        nan        nan        nan        nan   0.594802
   0.394684   0.209014        nan        nan        nan   0.684083   0.539167
   0.532304   0.332909   0.305677   0.437686   0.456637   0.408559   0.393593

Differences from the treams version

  • sphere_tmatrix takes keywords and the sphere's material; the exterior medium is vacuum unless given. treams takes one list of materials, ending with the exterior.
  • average_cross_sections returns named scattering and extinction values; treams reads xs_sca_avg and xs_ext_avg.
  • select(SphericalBasis.default(1)) keeps the dipole modes; treams indexes with tm[swb_lmax1].
  • The script gives plane_wave a unit direction, which does not depend on k0; treams passes the wave vector [0, 0, k0], which treams-rs also accepts.
  • scatter returns the scattered wave; treams writes tm @ inc.expand(tm.basis).
  • The intensity map starts as NaN, so points inside the sphere stay NaN.