Tutorial: Own structures and materials

01/2021: updated to pyGDM v1.1+

This tutorial demonstrates how user-defined structure models and material classes can be implemented in pyGDM.

Note: Structures can also be loaded from text-files as a list of (x,y,z) coordinates. Material dispersions can also be loaded from files containing tabulated data via the “materials.fromFile” class.

First load pyGDM

[1]:
import numpy as np
import matplotlib.pyplot as plt

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

Implement a structure generator

For this example we implement a simple structure generator, which will return the coordinates of a volume discretized cube. The important part is the definition of “condition_cube”, which will be called to choose whether a meshpoint coordinate lies inside or outside the volume of the structure.

[2]:
from pyGDM2.structures import _meshCubic, _meshHexagonalCompact

def structure_cube(step, L, mesh='cube', orientation_hex=1):
    """Generate volume discretized cube. Side length 'L', stepsize 'step'"""

    ## condition for a point being inside the cube volume
    def condition_cube(xi,yi,zi):
        return (0<xi<=L and 0<yi<=L and 0<zi<=L)

    ## call mesher
    if mesh == 'cube':
        sp = _meshCubic( (-2, int(L)+2), (-2, int(L)+2), (0, int(L)+1),
                        condition_cube)
    elif mesh == 'hex':
        sp = _meshHexagonalCompact( (-2, int(L)+2), (-2, int(L)+2), (-2, int(L)+2),
                                    condition_cube, orientation_hex)

    ## convert to numpy array and scale with discretization stepsize
    sp = np.array(sp, dtype=np.float)
    sp *= step

    return sp

Let’s test the new structure model.

[3]:
step = 10.0
geometry = structure_cube(step, L=5, mesh='cube')
geometry_hex = structure_cube(step, L=5, mesh='hex')


plt.subplot(121, aspect='equal')
pt = visu.structure(geometry, tit='cube, cubic mesh', show=0)
plt.subplot(122, aspect='equal')
pt = visu.structure(geometry_hex, tit='hex. mesh also works', show=1)
/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...")
../_images/tutorials_extend_01_own_structure_and_material_5_1.png

Implement a material dispersion class

To demonstrate the material class, we implement a simple material with constant refractive index of \(n=2\).

Note that the complex permittivity (\(\epsilon\)) is expected by pyGDM.

[4]:
class refindex_n2(object):
    def __init__(self):
        self.__name__ = 'constant index, n=2'

    def epsilon(self, wavelength):
        """return the index n=2 (epsilon=4)"""
        eps = complex(4.0)
        return eps


## test the dispersion
material_n2 = refindex_n2()
la = np.linspace(400, 1000, 50)


plt.plot(la, [material_n2.epsilon(l).real for l in la], label='real part')
plt.plot(la, [material_n2.epsilon(l).imag for l in la], label='imag. part')

plt.xlabel("wavelength (nm)")
plt.ylabel("epsilon")
plt.legend()
plt.show()
../_images/tutorials_extend_01_own_structure_and_material_7_0.png

Test the own model in a simulation

Now let’s test the own model and dispersion in a simple simulation

[5]:
## setup structure and disperision
struct = structures.struct(step, geometry, material_n2)

field_generator = fields.planewave
efield = fields.efield(field_generator, wavelengths=[400],
                       kwargs=dict(theta = [0.0]))

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

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

## run the scattering simulation
sim.scatter()

## plot the near-field (index 1: total electric field) in a plane close to the cube
r_probe = tools.generate_NF_map(-500,500,50, -500,500,50, Z0=100)
nearfield = linear.nearfield(sim, field_index=0, r_probe=r_probe)[1]

im = visu.vectorfield_color(nearfield, show=0)
plt.colorbar(im)
plt.show()
structure initialization - automatic mesh detection: cube
structure initialization - consistency check: 125/125 dipoles valid
/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=400.00nm - setup: EE 3115.7ms, inv.: 149.0ms, repropa.: 862.8ms (1 field configs), tot: 4127.9ms
../_images/tutorials_extend_01_own_structure_and_material_9_3.png