Dielectric nano-sphere, constant index

01/2021: updated to pyGDM v1.1+

Comparing pyGDM to Mie theory for a constant ref.-index dielectric sphere (n=2, D=300nm).

At first we load the modules:

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
import copy

## --- load pre-calculated Mie-data
wl_mie, qext_mie, qsca_mie = np.loadtxt("scat_mie_n2_D300nm.txt").T

Simulation setup

We now create several simulations which calculate scattering spectra for nano-spheres using different discretization density and meshes. Then we compare the results to Mie theory.

#==============================================================================
# pyGDM setup
#==============================================================================
## ---------- Setup scale-factors for differntly fine meshes
## --- cubic mesh structures
factors_cube = [0.8, 1.0, 1.25]

## --- hexagonal mesh structures
factors_hex = [0.75, 1.0, 1.2]


#==============================================================================
# Setup incident field
#==============================================================================
field_generator = fields.planewave
## log-interval spectrum (denser at low lambda):
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 environment (vacuum) and geometries (spheres D=300nm, n=2)
#==============================================================================
n1, n2 = 1.0, 1.0     # vacuum environment
dyads = propagators.DyadsQuasistatic123(n1, n2)

material = materials.dummy(2.0)

## --- cubic mesh
sim_cube = []
for fact in factors_cube:
    step = 37.5/fact
    radius = 4.*fact     # 37.5*4 = 150nm  -->  D=300nm
    geometry = structures.sphere(step, R=radius, mesh='cube')

    struct = structures.struct(step, geometry, material)

    sim = core.simulation(struct, efield, dyads)
    sim_cube.append( copy.deepcopy(sim) )

## --- hexagonal mesh
sim_hex = []
for fact in factors_hex:
    step = 37.5/fact
    radius = 4.*fact
    geometry = structures.sphere(step, R=radius, mesh='hex', ORIENTATION=2)

    struct = structures.struct(step, geometry, material)

    sim = core.simulation(struct, efield, dyads)
    sim_hex.append( copy.deepcopy(sim) )

structure initialization - automatic mesh detection: cube
structure initialization - consistency check: 171/171 dipoles valid
structure initialization - automatic mesh detection: cube
structure initialization - consistency check: 305/305 dipoles valid
structure initialization - automatic mesh detection: cube
structure initialization - consistency check: 619/619 dipoles valid
structure initialization - automatic mesh detection: hex
structure initialization - consistency check: 195/195 dipoles valid
structure initialization - automatic mesh detection: hex
structure initialization - consistency check: 425/425 dipoles valid
structure initialization - automatic mesh detection: hex
structure initialization - consistency check: 763/763 dipoles valid

Run the simulations

## --- cubic mesh
a_ext_cube = []
labels_cube = []
for sim in sim_cube:
    plt.figure()
    plt.title('(cube) N_dipoles = {}'.format(len(sim.struct.geometry)))
    visu.structure(sim)
    print('(cube) ----- N_dipoles =', len(sim.struct.geometry), end='')

    sim.scatter(verbose=False)

    ## extinction spectrum
    field_kwargs = tools.get_possible_field_params_spectra(sim)[0]
    wl, spec = tools.calculate_spectrum(sim, field_kwargs, linear.extinct)
    a_ext = spec.T[0]
    a_geo = tools.get_geometric_cross_section(sim)

    a_ext_cube.append(a_ext/a_geo)
    labels_cube.append(len(sim.struct.geometry))
    print('done.')


## --- hexagonal mesh
a_ext_hex = []
labels_hex = []
for sim in sim_hex:
    plt.figure()
    plt.title('(hex) N_dipoles = {}'.format(len(sim.struct.geometry)))
    visu.structure(sim)
    print('(hex) ----- N_dipoles =', len(sim.struct.geometry), end='')

    sim.scatter(verbose=False)

    ## extinction spectrum
    field_kwargs = tools.get_possible_field_params_spectra(sim)[0]
    wl, spec = tools.calculate_spectrum(sim, field_kwargs, linear.extinct)
    a_ext = spec.T[0]
    a_geo = tools.get_geometric_cross_section(sim)

    a_ext_hex.append(a_ext/a_geo)
    labels_hex.append(len(sim.struct.geometry))
    print('done.')

/home/hans/.local/lib/python3.8/site-packages/pyGDM2/visu.py:49: UserWarning: 3D data. Falling back to XY projection...
  warnings.warn("3D data. Falling back to XY projection...")

(cube) ----- N_dipoles = 171

/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()

done.

(cube) ----- N_dipoles = 305done.

(cube) ----- N_dipoles = 619done.

(hex) ----- N_dipoles = 195done.

(hex) ----- N_dipoles = 425done.

(hex) ----- N_dipoles = 763done.

Plot the spectra

plt.figure(figsize=(11,4))


## --- cubic
plt.subplot(121)
plt.plot(wl_mie, qext_mie, 'r--', dashes=[2,1],label='Mie')

colors = plt.cm.Blues( np.linspace(0.5, 1.0, len(a_ext_cube)) )
for i, (ae, lab) in enumerate(zip(a_ext_cube, labels_cube)):
    c = colors[i]
    lab='N={}'.format(lab)
    plt.plot(wl, ae, color=c, lw=0.5)
    plt.scatter(wl, ae, marker='x', linewidth=1.5, color=c, label=lab)

plt.legend(loc='best', fontsize=12, title='cubic mesh')
plt.xlabel("wavelength (nm)")
plt.ylabel("extinction efficiency")
plt.xlim( [wl.min(), wl.max()] )


## --- hexagonal
plt.subplot(122)
plt.plot(wl_mie, qext_mie, 'r--', dashes=[2,1],label='Mie')

colors = plt.cm.Blues( np.linspace(0.5, 1.0, len(a_ext_hex)) )
for i, (ae, lab) in enumerate(zip(a_ext_hex, labels_hex)):
    c = colors[i]
    lab='N={}'.format(lab)
    plt.plot(wl, ae, color=c, lw=0.5)
    plt.scatter(wl, ae, marker='x', linewidth=1.5, color=c, label=lab)

plt.legend(loc='best', fontsize=12, title='haxagonal mesh')
plt.xlabel("wavelength (nm)")
plt.ylabel("extinction efficiency")
plt.xlim( [wl.min(), wl.max()] )


plt.tight_layout()
plt.show()

Looks as if already quite a few meshpoints are sufficient for an OK agreement with Mie theory.