Python implementation for medium evaluations

[smartalecH]

The python portion of #475 that is used to evaluate material dispersion profiles at a given frequency/wavelength. Useful for making sure your profile really is what you expect. Also great for doing several MPB simulations across a band.

Here is a simple test script that pulls the profiles of Ag, Au, Si, and Cu and compares them to their original sources (stored in the attached h5 file):

import numpy as np
import meep as mp
from meep.materials import Ag, Au, Cu, Si
from matplotlib import pyplot as plt
import h5py


dataNames     = ['Ag','Au','Cu','Si']
testMaterials = [Ag,Au,Cu,Si]
hf = h5py.File('materialData.hdf5', 'r')

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

for iter in range(len(testMaterials)):
    # Pull Meep data
    M = testMaterials[iter]
    f0 = M.valid_freq_range[0]
    f1 = M.valid_freq_range[1]
    wavelength = np.linspace(1/f0,1/f1,1000)
    freq = 1/wavelength
    eps  = M.get_eps(freq)
    n = [np.real(np.sqrt(num.x)) for num in eps]
    k = [np.imag(np.sqrt(num.x)) for num in eps]

    # Pull refractiveindex.info data
    ref = hf.get(dataNames[iter])
    wavelength_ref = ref[:,0]
    n_ref = ref[:,1]
    if M !=Si:
        k_ref = ref[:,2]

    linewidth = 2
    # Set up plot
    plt.subplot(421 + 2*iter)
    plt.plot(wavelength_ref,n_ref,color='Red',linewidth=linewidth,label='ref')
    plt.plot(wavelength,n,'--',color='Blue',linewidth=linewidth,label='meep')
    plt.title("{} (real)".format(dataNames[iter]))
    plt.legend()
    plt.xlabel('Wavelength ($\mu$m)')
    plt.grid(True)
    plt.ylabel('n')

    if M !=Si:
        plt.subplot(421 + 2*iter+1)
        plt.plot(wavelength_ref,k_ref,color='Red',linewidth=linewidth,label='ref')
        plt.plot(wavelength,k,'--',color='Blue',linewidth=linewidth,label='meep')
        plt.title("{} (imag)".format(dataNames[iter]))
        plt.legend()
        plt.grid(True)
        plt.xlabel('Wavelength ($\mu$m)')
        plt.ylabel('k')
    

plt.tight_layout()
plt.savefig('materialEval.png')
plt.show()

materialData.hdf5.zip

[smartalecH]

Thanks for the suggestions! The function epsilon() now outputs a list of matrix objects corresponding to the dielectric tensor at that particular frequency.

Here is a test case showing both the ordinary and extraordinary axes for lithium niobate and barium borate:

hf = h5py.File('materialData_ani.hdf5', 'r')

dataNames = ["BaB2","LiNO"]
testMaterials = [BaB2O4,LiNbO3]
names = ["BaB2O4","LiNbO3"]
plt.figure(figsize=(4,4))
for iter in range(len(testMaterials)):
    # Pull Meep data
    M = testMaterials[iter]
    f0 = M.valid_freq_range[0]
    f1 = M.valid_freq_range[1]
    wavelength = np.linspace(1/f0,1/f1,1000)
    freq = 1/wavelength
    eps  = M.epsilon(freq.tolist())
    ne = [np.real(np.sqrt(num.c3.z)) for num in eps]
    no = [np.real(np.sqrt(num.c1.x)) for num in eps]

    # Pull refractiveindex.info data
    ref = hf.get(dataNames[iter] + '_e')
    wavelength_ref = ref[:,0]
    ne_ref = ref[:,1]
    ref = hf.get(dataNames[iter] + '_o')
    wavelength_ref = ref[:,0]
    no_ref = ref[:,1]

    linewidth = 2
    # Set up plot
    plt.subplot(221 + 2*iter)
    plt.plot(wavelength_ref,ne_ref,color='Red',linewidth=linewidth,label='ref')
    plt.plot(wavelength,ne,'--',color='Blue',linewidth=linewidth,label='meep')
    plt.title("{}".format(dataNames[iter]))
    plt.legend()
    plt.xlabel('Wavelength ($\mu$m)')
    plt.grid(True)
    plt.ylabel('n$_e$')


    plt.subplot(221 + 2*iter+1)
    plt.plot(wavelength_ref,no_ref,color='Red',linewidth=linewidth,label='ref')
    plt.plot(wavelength,no,'--',color='Blue',linewidth=linewidth,label='meep')
    plt.title("{}".format(dataNames[iter]))
    plt.legend()
    plt.grid(True)
    plt.xlabel('Wavelength ($\mu$m)')
    plt.ylabel('n$_o$')
    

plt.tight_layout()
plt.savefig('materialEval_ani.png')

plt.show()

The previous simulations also still work.

I still haven't refactored the function into _get_epsmu(tensor, susceptibilities, conductivity) quite yet. Where would you like that function to reside?

Also, the Matrix class seems to have trouble accessing its elements with the bracket notation, even with the __getitem__ overload. That's why I resorted to the odd, num.c3.z notation. The class is also missing copy overloads. We may want to look into that.

materialData_ani.zip

[stevengj]

Here is some anisotropic data on liquid crystal, for example: https://www.tandfonline.com/doi/full/10.1080/15421400600655816?scroll=top&needAccess=true&

[stevengj]

For a test, just comparing the output to two hard-coded values should be sufficient.

[stevengj]

To handle scalars and numpy arrays, you could do something like:

if hasattr(freq, "__iter__"):
    freqs = freq
else:
    freqs = [freq]

then, at the end, you can do something like

if hasattr(freq, "__iter__"):
    return result
else:
    return result[1]

You probably also want to make the result an Nx3x3 numpy array rather than a Python list.

[smartalecH]

I just pushed an update that refactored the method to use a common _get_epsmu() function. I vectorized the computation using numpy arrays (I'm not sure there's actually any gains here).

The function can take either a scalar, list, or numpy array. It computes everything and outputs a numpy array with size (N,3,3) where N is the number of frequency points passed.

I added a test script that tests Si, LiNO3, and Au. I'm still looking for a decent material to test nondiagonal components.

I also added one more Lorentzian to the Au and Ag entries to match the cited publication.

[smartalecH]

My latest commit leverages broadcasting to streamline the computation.

I've also modified the Matrix constructor to allow the user to either specify three column vectors (as in the current implementation) or to specify a diagonal vector and an off-diagonal vector. If we decide to keep this, we should modify the docs too.

[oskooi]

Note that you need to add the test to python/Makefile.am; otherwise it is not part of the test suite and will not run during make check.

[oskooi]

This additional Lorentzian term corresponds to a resonance at a wavelength of 0.061 μm (or 61 nm) that is so far outside the optical range that it should not affect results but it will increase the computational expense due to the additional storage and time stepping of auxiliary fields. (This is why it was originally left out.) Is it really necessary?

[smartalecH]

Right, but it seems to be a pretty broad resonance because even though it's so far away, it does alter the shape of the resonances between 200 nm and 300 nm. After that, however, it's negligible.

Here's a comparison:

The 5 resonance implementation matches refractiveindex.info's raw data points.

[oskooi]

This Lorentzian resonance is at 0.093 μm (or 93 nm). Similar to the additional Ag term; is it really necessary or should we just reduce the tolerance of the test slightly (i.e., matching the results to 4 decimal places rather than 6)?

[smartalecH]

Same as before. Let me know if you want me to remove it.