Generalized Polarizabilties - farfield patterns

New functionalities in ``multipole`` added in v1.1.3

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

Here we demonstrate how to use the generalized polarizabilities (GP) to calculate scattering farfield radiation patterns. The GP formalism is based on the exact multipole expansion [1]. It is described in [2].

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

[2] Majorel et al. Generalized polarizabilites for an exact multipole analysis of complex nanostructures under inhomogeneous illumination. arXiv (2022)

Setup and run simulation

from __future__ import print_function, division

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


# =============================================================================
# sim setup
# =============================================================================
method = 'lu'
#method = 'cupy'


## ----- set up single structure
step = 25
mat = materials.dummy(4)

geo = structures.rect_wire(step, L=8, W=8, H=8, mesh='cube')
struct = structures.struct(step, geo, mat)
struct = structures.center_struct(struct, which_axis=['x','y','z'])


## ----- illumination
wavelengths = [595, 610, 700, 800, 880, 950, 1020, 1200]
field_generator = fields.plane_wave
field_kwargs = dict(E_s=0, E_p=1, 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))


## ----- environment: vacuum
n1 = 1.0
dyads = propagators.DyadsQuasistatic123(n1)


## ----- simulation1
sim = core.simulation(struct=struct, efield=efield, dyads=dyads)
visu.structure(sim, tit='N dipoles = {}'.format(len(sim.struct.geometry)))



# =============================================================================
# run full simulation for comparison
# =============================================================================
sim.scatter(method=method)

structure initialization - automatic mesh detection: cube
structure initialization - consistency check: 512/512 dipoles valid

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

timing for wl=595.00nm - setup: EE 3098.2ms, inv.: 48.6ms, repropa.: 3051.9ms (1 field configs), tot: 6198.8ms
timing for wl=610.00nm - setup: EE 21.0ms, inv.: 127.0ms, repropa.: 16.1ms (1 field configs), tot: 164.3ms
timing for wl=700.00nm - setup: EE 30.5ms, inv.: 64.9ms, repropa.: 24.2ms (1 field configs), tot: 119.7ms
timing for wl=800.00nm - setup: EE 42.1ms, inv.: 120.9ms, repropa.: 20.5ms (1 field configs), tot: 183.6ms
timing for wl=880.00nm - setup: EE 50.6ms, inv.: 106.6ms, repropa.: 16.1ms (1 field configs), tot: 173.3ms
timing for wl=950.00nm - setup: EE 29.2ms, inv.: 81.7ms, repropa.: 16.0ms (1 field configs), tot: 127.0ms
timing for wl=1020.00nm - setup: EE 29.7ms, inv.: 52.9ms, repropa.: 16.1ms (1 field configs), tot: 98.8ms
timing for wl=1200.00nm - setup: EE 46.4ms, inv.: 42.9ms, repropa.: 16.1ms (1 field configs), tot: 105.5ms

1

calculate farfield patterns and plot

Now we calculate the farf-field patterns via

• the full simulations

• the full sim. multipole expansion

• the generalized polarizabilities

# =============================================================================
# multipole far-field pattern
# =============================================================================
test_wl_list  = [595, 610, 700, 800, 880, 950, 1020,1200]

for wl_eval in test_wl_list:
    field_index = tools.get_closest_field_index(sim, dict(wavelength=wl_eval))
    field_params = tools.get_field_indices(sim)[field_index]
    wl = field_params['wavelength']

    ## 2D radiation pattern config
    r       = 100000.   # nm
    tetamin = 0
    tetamax = 2*np.pi
    phimin  = 0
    phimax  = 0
    Nteta   = 200
    Nphi    = 1

    ## full sim scattering pattern
    tetal, phil, Is_sim, It_sim, I0_sim = linear.farfield(
                                        sim, field_index=field_index,
                                        phimin=phimin, phimax=phimax,
                                        tetamin=tetamin, tetamax=tetamax, r=r,
                                        Nteta=Nteta, Nphi=Nphi)

    ## exact multipoles from internal field
    tetal, phil, Is_mp, It_mp, I0_mp = multipole.farfield(
                                        sim,
                                        field_index=field_index,
                                        phimin=phimin, phimax=phimax,
                                        tetamin=tetamin, tetamax=tetamax, r=r,
                                        Nteta=Nteta, Nphi=Nphi,
                                        which_moments=['p', 'm', 'qe', 'qm'])

    ## exact multipoles from generalized polarizability
    p = multipole.eval_generalized_polarizability_p(sim, field_index, method=method)
    m = multipole.eval_generalized_polarizability_m(sim, field_index, method=method)
    qe = multipole.eval_generalized_polarizability_qe(sim, field_index, method=method)
    qm = multipole.eval_generalized_polarizability_qm(sim, field_index, method=method)

    tetal, phil, Is_mp_gp = multipole._multipole_farfield(
                                        wavelength=wl, eps_env=n1**2,
                                        r0=sim.r0,
                                        p=[p], m=[m], qe=[qe], qm=[qm],
                                        r_p=None, r_m=None, r_qe=None, r_qm=None,
                                        phimin=phimin, phimax=phimax,
                                        tetamin=tetamin, tetamax=tetamax, r=r,
                                        Nteta=Nteta, Nphi=Nphi)

    ## --- plot pattern section
    plt.figure(figsize=(4,4))
    plt.subplot(polar=True)
    plt.title(f"$\lambda$={wl}nm", x=0.75, y=1.07)
    plt.gca().set_theta_zero_location("S")

    plt.plot(np.pi - tetal.T[0], Is_sim.T[0], color='C0',dashes=[3,2], lw=2, label='full sim')

    plt.plot(np.pi - tetal.T[0], Is_mp, color='C2',dashes=[1,1], lw=1, label='exact multipoles')

    plt.plot((np.pi - tetal.T[0])[::-2], Is_mp_gp[::-2],   # plot every second point
             color='C2', marker='x', markersize=3, lw=0, label='via general. pola.')

    plt.legend(loc=[-.05, 1], fontsize=8)
    plt.tight_layout()
    plt.show()

wl=595nm. calc. K: 2313.1ms.  electric... 925.4ms.  magnetic... 165.9ms. Done.

wl=610nm. calc. K: 1688.2ms.  electric... 934.4ms.  magnetic... 172.0ms. Done.

wl=700nm. calc. K: 832.9ms.  electric... 894.1ms.  magnetic... 169.7ms. Done.

wl=800nm. calc. K: 152.0ms.  electric... 913.2ms.  magnetic... 171.4ms. Done.

wl=880nm. calc. K: 1113.6ms.  electric... 921.2ms.  magnetic... 167.3ms. Done.

wl=950nm. calc. K: 883.5ms.  electric... 902.9ms.  magnetic... 166.0ms. Done.

wl=1020nm. calc. K: 864.6ms.  electric... 908.8ms.  magnetic... 164.2ms. Done.

wl=1200nm. calc. K: 247.1ms.  electric... 926.2ms.  magnetic... 163.1ms. Done.