Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Python implementation for medium evaluations #862

Merged
merged 8 commits into from
May 14, 2019
Merged

Python implementation for medium evaluations #862

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
6051db7
Python implementation for medium evaluations
smartalecH May 7, 2019
61f95d6
Enabled full Chi1 tensor, small naming conventions
smartalecH May 7, 2019
4345afe
refactored method with test and updated materials
smartalecH May 9, 2019
b633fa2
refactored matrix init and streamlined computation with broadcasting
smartalecH May 10, 2019
1cd1e11
squeeze output for scalar freq inputs
smartalecH May 10, 2019
371f9cc
Add test to makefile
smartalecH May 10, 2019
a7dd4cd
complex hermitian Matrix init
smartalecH May 10, 2019
373795c
docs
smartalecH May 14, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
53 changes: 53 additions & 0 deletions python/geom.py
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,9 @@ def __init__(self, epsilon_diag=Vector3(1, 1, 1),
H_chi2_diag=Vector3(),
H_chi3_diag=Vector3(),
D_conductivity_diag=Vector3(),
D_conductivity_offdiag=Vector3(),
B_conductivity_diag=Vector3(),
B_conductivity_offdiag=Vector3(),
epsilon=None,
index=None,
mu=None,
Expand Down Expand Up @@ -211,7 +213,9 @@ def __init__(self, epsilon_diag=Vector3(1, 1, 1),
self.H_chi2_diag = H_chi2_diag
self.H_chi3_diag = H_chi3_diag
self.D_conductivity_diag = D_conductivity_diag
self.D_conductivity_offdiag = D_conductivity_offdiag
self.B_conductivity_diag = B_conductivity_diag
self.B_conductivity_offdiag = D_conductivity_offdiag
self.valid_freq_range = valid_freq_range

def transform(self, m):
Expand All @@ -235,6 +239,39 @@ def transform(self, m):
for s in self.H_susceptibilities:
s.transform(m)

def epsilon(self,freq):
return self._get_epsmu(self.epsilon_diag, self.epsilon_offdiag, self.E_susceptibilities, self.D_conductivity_diag, self.D_conductivity_offdiag, freq)

def mu(self,freq):
return self._get_epsmu(self.mu_diag, self.mu_offdiag, self.H_susceptibilities, self.B_conductivity_diag, self.B_conductivity_offdiag, freq)

def _get_epsmu(self, diag, offdiag, susceptibilities, conductivity_diag, conductivity_offdiag, freq):
# Clean the input
if type(freq) is not np.ndarray:
freqs = np.array(freq)
else:
freqs = np.squeeze(freq)

numFreqs = freqs.size

# Initialize with instantaneous dielectric tensor
epsmu = np.tile(np.array(Matrix(mp.Vector3(diag.x, offdiag.x, offdiag.y),
mp.Vector3(offdiag.x, diag.y, offdiag.z),
mp.Vector3(offdiag.y, offdiag.z, diag.z))),(numFreqs,1,1))
stevengj marked this conversation as resolved.
Show resolved Hide resolved

# Iterate through susceptibilities
for i_sus in range(len(susceptibilities)):
epsmu = epsmu + susceptibilities[i_sus].eval_susceptibility(freqs)

# Account for conductivity term
conductivity = np.tile(np.array(Matrix(mp.Vector3(conductivity_diag.x, conductivity_offdiag.x, conductivity_offdiag.y),
mp.Vector3(conductivity_offdiag.x, conductivity_diag.y, conductivity_offdiag.z),
mp.Vector3(conductivity_offdiag.y, conductivity_offdiag.z, conductivity_diag.z))),(numFreqs,1,1))
epsmu = (1 + 1j/np.tile(freqs,(3,3,1)).T * conductivity) * epsmu
stevengj marked this conversation as resolved.
Show resolved Hide resolved

# Convert list matrix to 3D numpy array size [freqs,3,3]
return np.array(epsmu)
stevengj marked this conversation as resolved.
Show resolved Hide resolved


class Susceptibility(object):

Expand All @@ -257,6 +294,14 @@ def __init__(self, frequency=0.0, gamma=0.0, **kwargs):
super(LorentzianSusceptibility, self).__init__(**kwargs)
self.frequency = frequency
self.gamma = gamma

def eval_susceptibility(self,freq):

sigma = np.tile(np.array(Matrix(mp.Vector3(self.sigma_diag.x, self.sigma_offdiag.x, self.sigma_offdiag.y),
mp.Vector3(self.sigma_offdiag.x, self.sigma_diag.y, self.sigma_offdiag.z),
mp.Vector3(self.sigma_offdiag.y, self.sigma_offdiag.z, self.sigma_diag.z))),(freq.size,1,1))

return np.tile(self.frequency*self.frequency / (self.frequency*self.frequency - freq*freq - 1j*self.gamma*freq),(3,3,1)).T * sigma


class DrudeSusceptibility(Susceptibility):
Expand All @@ -265,6 +310,14 @@ def __init__(self, frequency=0.0, gamma=0.0, **kwargs):
super(DrudeSusceptibility, self).__init__(**kwargs)
self.frequency = frequency
self.gamma = gamma

def eval_susceptibility(self,freq):

sigma = np.tile(np.array(Matrix(mp.Vector3(self.sigma_diag.x, self.sigma_offdiag.x, self.sigma_offdiag.y),
mp.Vector3(self.sigma_offdiag.x, self.sigma_diag.y, self.sigma_offdiag.z),
mp.Vector3(self.sigma_offdiag.y, self.sigma_offdiag.z, self.sigma_diag.z))),(freq.size,1,1))

return np.tile(-self.frequency*self.frequency / (freq*(freq + 1j*self.gamma)),(3,3,1)).T * sigma


class NoisyLorentzianSusceptibility(LorentzianSusceptibility):
Expand Down
14 changes: 12 additions & 2 deletions python/materials.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -258,12 +258,17 @@
Ag_frq4 = 9.083*eV_um_scale # 0.137 um
Ag_gam4 = 0.916*eV_um_scale
Ag_sig4 = Ag_f4*Ag_plasma_frq**2/Ag_frq4**2
Ag_f5 = 5.646
Ag_frq5 = 20.29*eV_um_scale
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Collaborator Author

@smartalecH smartalecH May 10, 2019 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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:

copareResonance

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

Ag_gam5 = 2.419*eV_um_scale
Ag_sig5 = Ag_f5*Ag_plasma_frq**2/Ag_frq5**2

Ag_susc = [mp.DrudeSusceptibility(frequency=Ag_frq0, gamma=Ag_gam0, sigma=Ag_sig0),
mp.LorentzianSusceptibility(frequency=Ag_frq1, gamma=Ag_gam1, sigma=Ag_sig1),
mp.LorentzianSusceptibility(frequency=Ag_frq2, gamma=Ag_gam2, sigma=Ag_sig2),
mp.LorentzianSusceptibility(frequency=Ag_frq3, gamma=Ag_gam3, sigma=Ag_sig3),
mp.LorentzianSusceptibility(frequency=Ag_frq4, gamma=Ag_gam4, sigma=Ag_sig4)]
mp.LorentzianSusceptibility(frequency=Ag_frq4, gamma=Ag_gam4, sigma=Ag_sig4),
mp.LorentzianSusceptibility(frequency=Ag_frq5, gamma=Ag_gam5, sigma=Ag_sig5)]

Ag = mp.Medium(epsilon=1.0, E_susceptibilities=Ag_susc, valid_freq_range=metal_range)

Expand Down Expand Up @@ -291,12 +296,17 @@
Au_frq4 = 4.304*eV_um_scale # 0.288 um
Au_gam4 = 2.494*eV_um_scale
Au_sig4 = Au_f4*Au_plasma_frq**2/Au_frq4**2
Au_f5 = 4.384
Au_frq5 = 13.32*eV_um_scale
Copy link
Collaborator

@oskooi oskooi May 10, 2019 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Au_gam5 = 2.214*eV_um_scale
Au_sig5 = Au_f5*Au_plasma_frq**2/Au_frq5**2

Au_susc = [mp.DrudeSusceptibility(frequency=Au_frq0, gamma=Au_gam0, sigma=Au_sig0),
mp.LorentzianSusceptibility(frequency=Au_frq1, gamma=Au_gam1, sigma=Au_sig1),
mp.LorentzianSusceptibility(frequency=Au_frq2, gamma=Au_gam2, sigma=Au_sig2),
mp.LorentzianSusceptibility(frequency=Au_frq3, gamma=Au_gam3, sigma=Au_sig3),
mp.LorentzianSusceptibility(frequency=Au_frq4, gamma=Au_gam4, sigma=Au_sig4)]
mp.LorentzianSusceptibility(frequency=Au_frq4, gamma=Au_gam4, sigma=Au_sig4),
mp.LorentzianSusceptibility(frequency=Au_frq5, gamma=Au_gam5, sigma=Au_sig5)]

Au = mp.Medium(epsilon=1.0, E_susceptibilities=Au_susc, valid_freq_range=metal_range)

Expand Down
46 changes: 46 additions & 0 deletions python/tests/medium_evaluations.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,46 @@

# medium_evaluations.py - Tests the evaluation of material permitivity profiles.
# Checks materials with lorentizian, drude, and non uniform diagonals.
# The extracted values are compared against actual datapoints pulled from
# refractiveindex.info.
# TODO:
# * check materials with off diagonal components
# * check magnetic profiles

from __future__ import division

import unittest
import meep as mp
import numpy as np


class TestMediumEvaluations(unittest.TestCase):

def test_medium_evaluations(self):
from meep.materials import Si, Au, LiNbO3

# Check Silicon
w0 = 1/1.357
w1 = 1/11.04
eps = Si.epsilon([w0,w1])
self.assertAlmostEqual(np.real(np.sqrt(eps[0,0,0])), 3.4975134228247, places=6)
self.assertAlmostEqual(np.real(np.sqrt(eps[1,0,0])), 3.4175012372164, places=6)

# Check Gold
w0 = 1/0.24797
w1 = 1/6.1992
eps = Au.epsilon([w0,w1])
self.assertAlmostEqual(np.real(np.sqrt(eps[0,0,0])), 1.0941, places=4)
self.assertAlmostEqual(np.real(np.sqrt(eps[1,0,0])), 5.0974, places=4)

# Check Lithium Niobate
w0 = 1/0.4
w1 = 1/5.0
eps = LiNbO3.epsilon([w0,w1])
self.assertAlmostEqual(np.real(np.sqrt(eps[0,0,0])), 2.4392710888511, places=6)
self.assertAlmostEqual(np.real(np.sqrt(eps[1,0,0])), 2.0508048794821, places=6)
self.assertAlmostEqual(np.real(np.sqrt(eps[0,2,2])), 2.332119801394, places=6)
self.assertAlmostEqual(np.real(np.sqrt(eps[1,2,2])), 2.0025312699385, places=6)

if __name__ == '__main__':
unittest.main()