Forward / backward resolved scattering

01/2021: updated to pyGDM v1.1+

In this example, we try to reproduce the directional visible light scattering from silicon spheres, reported by Fu et al. [1].

[1]: Fu, Y. H. et al.: Directional visible light scattering by silicon nanoparticles. Nat Commun 4, 1527 (2013) (https://doi.org/10.1038/ncomms2538)

from pyGDM2 import structures
from pyGDM2 import materials
from pyGDM2 import fields

from pyGDM2 import core
from pyGDM2 import propagators
from pyGDM2 import tools
from pyGDM2 import linear
from pyGDM2 import visu

import numpy as np
import matplotlib.pyplot as plt

Setting up the simulation

## --- Setup incident field
field_generator = fields.planewave
wavelengths = np.exp(np.linspace(np.log(300), np.log(1000), 30))
kwargs = dict(theta = [0.0])
efield = fields.efield(field_generator, wavelengths=wavelengths, kwargs=kwargs)


## --- Setup geometry (sphere D=150nm in vacuum)
scale_factor = 1.4
step = 18.75/scale_factor
radius = 4.*scale_factor
geometry = structures.sphere(step, R=radius, mesh='hex', ORIENTATION=2)
material = materials.silicon()
struct = structures.struct(step, geometry, material)

dyads = propagators.DyadsQuasistatic123(n1=1, n2=1)

sim = core.simulation(struct, efield, dyads)

structure initialization - automatic mesh detection: hex
structure initialization - consistency check: 1159/1159 dipoles valid

Run the simulation, get FW/BW scattering spectra

At first we run the main simulation core.scatter, then we calculate the scattering to the farfield separately for the upper and lower hemi-sphere:

## main simulation
E = core.scatter(sim, method='lu', verbose=True)

## FW and BW scattering spectrum
field_kwargs = tools.get_possible_field_params_spectra(sim)[0]
wl, scat_fw = tools.calculate_spectrum(sim, field_kwargs, linear.farfield,
                                      tetamin=np.pi/2., tetamax=np.pi,
                                      return_value='int_Es')
wl, scat_bw = tools.calculate_spectrum(sim, field_kwargs, linear.farfield,
                                      tetamin=0, tetamax=np.pi/2.,
                                      return_value='int_Es')

/home/hans/.local/lib/python3.8/site-packages/numba/core/dispatcher.py:237: UserWarning: Numba extension module 'numba_scipy' failed to load due to 'ValueError(No function '__pyx_fuse_0pdtr' found in __pyx_capi__ of 'scipy.special.cython_special')'.
  entrypoints.init_all()

timing for wl=300.00nm - setup: EE 7740.3ms, inv.: 4652.8ms, repropa.: 2363.3ms (1 field configs), tot: 14756.9ms
timing for wl=312.72nm - setup: EE 737.4ms, inv.: 1884.7ms, repropa.: 13.5ms (1 field configs), tot: 2636.4ms
timing for wl=325.97nm - setup: EE 261.1ms, inv.: 1791.0ms, repropa.: 13.4ms (1 field configs), tot: 2069.3ms
timing for wl=339.79nm - setup: EE 221.0ms, inv.: 1662.2ms, repropa.: 13.4ms (1 field configs), tot: 1897.6ms
timing for wl=354.20nm - setup: EE 268.1ms, inv.: 1899.1ms, repropa.: 13.8ms (1 field configs), tot: 2182.0ms
timing for wl=369.21nm - setup: EE 234.2ms, inv.: 3208.1ms, repropa.: 25.6ms (1 field configs), tot: 3469.7ms
timing for wl=384.86nm - setup: EE 659.2ms, inv.: 4652.4ms, repropa.: 27.2ms (1 field configs), tot: 5340.3ms
timing for wl=401.17nm - setup: EE 767.4ms, inv.: 4333.1ms, repropa.: 34.7ms (1 field configs), tot: 5136.6ms
timing for wl=418.18nm - setup: EE 761.8ms, inv.: 3195.4ms, repropa.: 13.9ms (1 field configs), tot: 3972.0ms
timing for wl=435.91nm - setup: EE 221.4ms, inv.: 1597.6ms, repropa.: 14.2ms (1 field configs), tot: 1833.8ms
timing for wl=454.39nm - setup: EE 221.5ms, inv.: 1758.3ms, repropa.: 14.4ms (1 field configs), tot: 1995.0ms
timing for wl=473.65nm - setup: EE 221.0ms, inv.: 1986.2ms, repropa.: 17.2ms (1 field configs), tot: 2225.2ms
timing for wl=493.72nm - setup: EE 230.5ms, inv.: 1729.8ms, repropa.: 13.4ms (1 field configs), tot: 1974.5ms
timing for wl=514.65nm - setup: EE 276.3ms, inv.: 4533.5ms, repropa.: 33.1ms (1 field configs), tot: 4844.3ms
timing for wl=536.47nm - setup: EE 863.4ms, inv.: 4576.2ms, repropa.: 33.7ms (1 field configs), tot: 5475.9ms
timing for wl=559.21nm - setup: EE 843.7ms, inv.: 4715.4ms, repropa.: 14.7ms (1 field configs), tot: 5575.4ms
timing for wl=582.92nm - setup: EE 225.0ms, inv.: 1852.6ms, repropa.: 15.9ms (1 field configs), tot: 2094.6ms
timing for wl=607.63nm - setup: EE 244.2ms, inv.: 1835.0ms, repropa.: 15.1ms (1 field configs), tot: 2095.2ms
timing for wl=633.38nm - setup: EE 224.6ms, inv.: 1677.1ms, repropa.: 13.4ms (1 field configs), tot: 1915.9ms
timing for wl=660.23nm - setup: EE 222.6ms, inv.: 1569.4ms, repropa.: 15.6ms (1 field configs), tot: 1808.2ms
timing for wl=688.22nm - setup: EE 219.9ms, inv.: 2657.2ms, repropa.: 35.0ms (1 field configs), tot: 2913.7ms
timing for wl=717.39nm - setup: EE 734.3ms, inv.: 5717.8ms, repropa.: 31.1ms (1 field configs), tot: 6485.1ms
timing for wl=747.80nm - setup: EE 808.6ms, inv.: 5602.3ms, repropa.: 28.8ms (1 field configs), tot: 6441.6ms
timing for wl=779.50nm - setup: EE 906.3ms, inv.: 2080.0ms, repropa.: 17.3ms (1 field configs), tot: 3004.3ms
timing for wl=812.55nm - setup: EE 257.8ms, inv.: 1811.4ms, repropa.: 16.2ms (1 field configs), tot: 2087.1ms
timing for wl=846.99nm - setup: EE 251.7ms, inv.: 1773.4ms, repropa.: 12.9ms (1 field configs), tot: 2038.9ms
timing for wl=882.90nm - setup: EE 304.1ms, inv.: 1889.3ms, repropa.: 13.3ms (1 field configs), tot: 2207.5ms
timing for wl=920.32nm - setup: EE 243.2ms, inv.: 2729.5ms, repropa.: 25.4ms (1 field configs), tot: 2999.8ms
timing for wl=959.33nm - setup: EE 662.7ms, inv.: 4702.6ms, repropa.: 31.2ms (1 field configs), tot: 5398.3ms
timing for wl=1000.00nm - setup: EE 777.3ms, inv.: 5263.1ms, repropa.: 29.4ms (1 field configs), tot: 6071.4ms

Plot the FW/BW spectra

Let’s see what scattering we get in both directions:

plt.figure()
plt.title(r"D=150nm Si sphere")

## --- scattering spectra FW & BW
plt.plot(wl, scat_fw, 'r', label='FW.')
plt.plot(wl, scat_bw, 'b', label='BW.')

plt.xlabel("wavelength (nm)")
plt.ylabel(r"scat. cross-section (nm^2)")
plt.xlim(400,800)

plt.plot([0], [0], color='g', dashes=[2,1], label='FW/BW')  # for legend entry
plt.legend(loc='center', frameon=False, ncol=1, bbox_to_anchor=(0.83, 0.3))


## --- logscale FW/BW ratio on right y-axis
plt.twinx()
plt.plot(wl, scat_fw/scat_bw, color='g', dashes=[2,1])
plt.plot([400,800], [1,1], color='k', lw=0.5, dashes=[2,1], label='FW/BW=1')

plt.ylabel(r"FW/BW ratio", rotation=270, labelpad=8, color='g')
plt.yscale('log')
plt.xlim( [400, 800] )
plt.ylim( [0.1, 10] )
plt.legend(loc='center', frameon=False, ncol=1, bbox_to_anchor=(0.8, 0.55))


plt.tight_layout()
plt.show()

Comparing this with the paper [1], cited at the very top, this looks pretty much exactly like their results (let’s say like a pixelated version of it…).