Multipole Decomposition 1

New module ``multipole`` in v1.1.2

!!CAUTION!! : ``multipole`` is a new functionality. Please report possible problems and errors.

Here we demonstrate the exact multipole decomposition of the internal fields following reference [1]. The formalism allows calculating the extinction and scattering spectra of arbitrary shaped nanostructures.

In this example we reproduce Mie theory for nano-spheres of various lossless and lossy materials and different sizes.

[1] Alaee, R., Rockstuhl, C. and Fernandez-Corbaton, I. An electromagnetic multipole expansion beyond the long-wavelength approximation. Optics Communications 407, 17-21 (2018)

Note: This example requires the package PyMieScatt (https://pymiescatt.readthedocs.io/) for the calculation of the Mie spectra.

Simulation setup

import numpy as np
import matplotlib.pyplot as plt
import copy

from pyGDM2 import structures
from pyGDM2 import materials
from pyGDM2 import fields
from pyGDM2 import tools
from pyGDM2 import propagators
from pyGDM2 import linear
from pyGDM2 import core
from pyGDM2 import visu
from pyGDM2 import multipole

import PyMieScatt as ps


## some general config
n1 = n2 = n3 = 1.0
# method = 'lu'
method = 'cupy'
spacing = 10000
mesh = 'hex'


## specific simulations to run
cases = [
    dict(D=220, mat=materials.dummy(3.), step=15, name='lossless: n=3, D=220nm'),
    dict(D=175, mat=materials.germanium(), step=10, name='lossy: Germanium, D=175nm'),
    dict(D=50, mat=materials.gold(), step=4, name='plasmonic: Gold, D=50nm'),
    ]

Calculate Mie and run simulations

We calculate the Mie spectra and run the simulations + their multipole analysis. finally we plot the scattering and extinction for each case.

for case in cases:
    #%% --- Mie theory
    name = case['name']
    D = case['D']
    step = case['step']
    mat = case['mat']


    wavelengths = np.arange(400.0,900.0, 10)
    csca_Mie = []
    cext_Mie = []
    for wl in wavelengths :
        Que = ps.MieQ(mat.epsilon(wl)**0.5, wl, D, nMedium=n2, asCrossSection=True, asDict=True)
        csca_Mie.append(Que['Csca'])
        cext_Mie.append(Que['Cext'])


    #%% --- setup GDM sim
    ## setup structure and simulation
    dyads = propagators.DyadsQuasistatic123(n1, n2, n3, spacing)

    ## set up single sphere structure
    geo = structures.sphere(step, R=D/(2*step), mesh=mesh)
    struct_single = structures.struct(step, geo, mat)

    ## illumination
    field_generator = fields.plane_wave
    field_kwargs = dict(E_s=1, E_p=0, inc_angle=0, inc_plane='xz') # lin-pol X, normal incidence from below

    efield = fields.efield(field_generator, wavelengths=wavelengths, kwargs=copy.deepcopy(field_kwargs))
    sim = core.simulation(struct=struct_single, efield=efield, dyads=dyads)

    #visu.structure(sim, tit=name + f' - Ndp:{len(geo)}')
    print(name + f' - Ndp:{len(geo)}')


    #%% --- run sim and do multipole analysis
    sim.scatter(method=method, verbose=False)

    wl, ext_spec = tools.calculate_spectrum(sim, 0, linear.extinct)
    ex, sc, ab = ext_spec.T

    wl, spec_pm_ex = tools.calculate_spectrum(sim, 0, multipole.extinct)
    ext_p_ex, ext_m_ex, ext_Q_ex, ext_Qm_ex = spec_pm_ex.T

    wl, mpsc_spec = tools.calculate_spectrum(sim, 0, multipole.scs)
    sc_p, sc_m, sc_Qe, sc_Qm = mpsc_spec.T


    #%% --- plot
    plt.figure(figsize=(14,5))

    plt.subplot(121, title=name)
    plt.plot(wl, ex, lw=1, label='full ext.')
    plt.plot(wl, ext_p_ex, color='C1', label='ext-p')
    plt.plot(wl, ext_m_ex, color='C2', label='ext-m')
    plt.plot(wl, ext_Q_ex, color='C3', label='ext-Qe')
    plt.plot(wl, ext_Qm_ex, color='C4', label='ext-Qm')
    plt.plot(wl, (ext_p_ex + ext_m_ex+ext_Qm_ex+ext_Q_ex), c='C6', dashes=[2,2], label='multipole sum')

    plt.plot(wl, cext_Mie, color='k', lw=1, label='Mie', dashes=[1,1])

    plt.xlabel('wavelength (nm)')
    plt.ylabel('extinction cross section (nm^2)')
    ylim = plt.ylim()
    plt.legend()


    plt.subplot(122, title=name)
    plt.plot(wl, sc, lw=1, label='full scat.')
    plt.plot(wl, sc_p, c='C1', label='p')
    plt.plot(wl, sc_m, c='C2', label='m')
    plt.plot(wl, sc_Qe, c='C3', label='Qe')
    plt.plot(wl, sc_Qm, c='C4', label='Qm')
    plt.plot(wl, (sc_p+sc_m+sc_Qm+sc_Qe), c='C6', dashes=[2,2], label='multipole sum')
    plt.plot(wl, csca_Mie, color='k', lw=1, label='Mie', dashes=[1,1])

    plt.xlabel('wavelength (nm)')
    plt.ylabel('scattering cross section (nm^2)')
    plt.ylim(ylim)
    plt.legend()
    plt.show()

structure initialization - automatic mesh detection: hex
structure initialization - consistency check: 2517/2517 dipoles valid
lossless: n=3, D=220nm - Ndp:2517

/home/hans/.local/lib/python3.8/site-packages/numba/core/dispatcher.py:289: UserWarning: Numba extension module 'numba_scipy' failed to load due to 'VersionConflict((scipy 1.7.1 (/home/hans/.local/lib/python3.8/site-packages), Requirement.parse('scipy<=1.6.2,>=0.16')))'.
  entrypoints.init_all()
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

structure initialization - automatic mesh detection: hex
structure initialization - consistency check: 4247/4247 dipoles valid
lossy: Germanium, D=175nm - Ndp:4247

/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

structure initialization - automatic mesh detection: hex
structure initialization - consistency check: 1555/1555 dipoles valid
plasmonic: Gold, D=50nm - Ndp:1555

/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/home/hans/.local/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)