From 6051db7e6b81588d96f513ba0eefded79455f6ff Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Tue, 7 May 2019 12:42:40 -0600 Subject: [PATCH 1/8] Python implementation for medium evaluations --- python/geom.py | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) mode change 100644 => 100755 python/geom.py diff --git a/python/geom.py b/python/geom.py old mode 100644 new mode 100755 index a3dc12715..1c2627b12 --- a/python/geom.py +++ b/python/geom.py @@ -234,6 +234,27 @@ def transform(self, m): for s in self.H_susceptibilities: s.transform(m) + def get_eps(self,freq): + + # Clean the input + freq = np.squeeze([freq]) + + # Compensate for scalar inputs + if len(freq.shape) == 0: + freq = np.array([freq]) + + # Initialize with instantaneous dielectric tensor, use numpy arrays for + # convenience and speed. + eps = np.array([self.epsilon_diag.x,self.epsilon_diag.y,self.epsilon_diag.z],dtype='complex128') + + # Iterate through and sum up susceptibilities + for k in range(len(self.E_susceptibilities)): + eps = eps + self.E_susceptibilities[k].eval_susceptibility(freq) + + # Account for conductivity term + eps = (1 + 1j*np.array( + [self.D_conductivity_diag.x,self.D_conductivity_diag.y,self.D_conductivity_diag.z])/freq.reshape(-1,1)) * eps + return [Vector3(eps[row,0],eps[row,1],eps[row,2]) for row in range(freq.shape[0])] class Susceptibility(object): @@ -257,6 +278,17 @@ 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): + + # Use numpy arrays for element-wise multiplication later + sigma = np.array([self.sigma_diag.x,self.sigma_diag.y,self.sigma_diag.z]) + eps = np.zeros((freq.shape[0],3),dtype='complex128') + + # Iterate through dimensions + for k in range(3): + eps[:,k] = (sigma[k]*self.frequency*self.frequency) / (self.frequency*self.frequency - freq*freq - 1j*self.gamma*freq) + return eps class DrudeSusceptibility(Susceptibility): @@ -265,6 +297,17 @@ 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): + + # Use numpy arrays for element-wise multiplication later + sigma = np.array([self.sigma_diag.x,self.sigma_diag.y,self.sigma_diag.z]) + eps = np.zeros((freq.shape[0],3),dtype='complex128') + + # Iterate through dimensions + for k in range(3): + eps[:,k] = (-sigma[k]*self.frequency*self.frequency) / (freq*(freq + 1j*self.gamma)) + return eps class NoisyLorentzianSusceptibility(LorentzianSusceptibility): From 61f95d61bd14cc64d7026d3398210a933fcfb2ff Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Tue, 7 May 2019 15:11:11 -0600 Subject: [PATCH 2/8] Enabled full Chi1 tensor, small naming conventions --- python/geom.py | 65 +++++++++++++++++++++++++------------------------- 1 file changed, 33 insertions(+), 32 deletions(-) diff --git a/python/geom.py b/python/geom.py index 1c2627b12..e5cb63071 100755 --- a/python/geom.py +++ b/python/geom.py @@ -211,6 +211,7 @@ 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 = Vector3() self.B_conductivity_diag = B_conductivity_diag self.valid_freq_range = valid_freq_range @@ -234,27 +235,33 @@ def transform(self, m): for s in self.H_susceptibilities: s.transform(m) - def get_eps(self,freq): - # Clean the input - freq = np.squeeze([freq]) - - # Compensate for scalar inputs - if len(freq.shape) == 0: - freq = np.array([freq]) + def epsilon(self,freq): - # Initialize with instantaneous dielectric tensor, use numpy arrays for - # convenience and speed. - eps = np.array([self.epsilon_diag.x,self.epsilon_diag.y,self.epsilon_diag.z],dtype='complex128') - - # Iterate through and sum up susceptibilities - for k in range(len(self.E_susceptibilities)): - eps = eps + self.E_susceptibilities[k].eval_susceptibility(freq) + # Clean the input + if not isinstance(freq, list): + raise ValueError('Frequencies must be supplied as a list') - # Account for conductivity term - eps = (1 + 1j*np.array( - [self.D_conductivity_diag.x,self.D_conductivity_diag.y,self.D_conductivity_diag.z])/freq.reshape(-1,1)) * eps - return [Vector3(eps[row,0],eps[row,1],eps[row,2]) for row in range(freq.shape[0])] + # Initialize with instantaneous dielectric tensor + eps = [Matrix(mp.Vector3(self.epsilon_diag.x, self.epsilon_offdiag.x, self.epsilon_offdiag.y), + mp.Vector3(self.epsilon_offdiag.x, self.epsilon_diag.y, self.epsilon_offdiag.z), + mp.Vector3(self.epsilon_offdiag.y, self.epsilon_offdiag.z, self.epsilon_diag.z)) for i in range(len(freq))] + + # Iterate through frequencies + for i_freq in range(len(freq)): + # Iterate through susceptibilities + for i_sus in range(len(self.E_susceptibilities)): + eps[i_freq] = eps[i_freq] + self.E_susceptibilities[i_sus].eval_susceptibility(freq[i_freq]) + # Account for conductivity term + D_conductivity = Matrix(mp.Vector3(self.D_conductivity_diag.x, self.D_conductivity_offdiag.x, self.D_conductivity_offdiag.y), + mp.Vector3(self.D_conductivity_offdiag.x, self.D_conductivity_diag.y, self.D_conductivity_offdiag.z), + mp.Vector3(self.D_conductivity_offdiag.y, self.D_conductivity_offdiag.z, self.D_conductivity_diag.z)) + ones_matrix = Matrix(mp.Vector3(1, 1, 1), + mp.Vector3(1, 1, 1), + mp.Vector3(1, 1, 1)) + eps[i_freq] = (ones_matrix + 1j/freq[i_freq] * D_conductivity) * eps[i_freq] + + return eps class Susceptibility(object): @@ -281,14 +288,11 @@ def __init__(self, frequency=0.0, gamma=0.0, **kwargs): def eval_susceptibility(self,freq): - # Use numpy arrays for element-wise multiplication later - sigma = np.array([self.sigma_diag.x,self.sigma_diag.y,self.sigma_diag.z]) - eps = np.zeros((freq.shape[0],3),dtype='complex128') + sigma = 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)) - # Iterate through dimensions - for k in range(3): - eps[:,k] = (sigma[k]*self.frequency*self.frequency) / (self.frequency*self.frequency - freq*freq - 1j*self.gamma*freq) - return eps + return self.frequency*self.frequency / (self.frequency*self.frequency - freq*freq - 1j*self.gamma*freq) * sigma class DrudeSusceptibility(Susceptibility): @@ -300,14 +304,11 @@ def __init__(self, frequency=0.0, gamma=0.0, **kwargs): def eval_susceptibility(self,freq): - # Use numpy arrays for element-wise multiplication later - sigma = np.array([self.sigma_diag.x,self.sigma_diag.y,self.sigma_diag.z]) - eps = np.zeros((freq.shape[0],3),dtype='complex128') + sigma = 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)) - # Iterate through dimensions - for k in range(3): - eps[:,k] = (-sigma[k]*self.frequency*self.frequency) / (freq*(freq + 1j*self.gamma)) - return eps + return -self.frequency*self.frequency / (freq*(freq + 1j*self.gamma)) * sigma class NoisyLorentzianSusceptibility(LorentzianSusceptibility): From 4345afeea73543031ba293e936070dfc71e77c6a Mon Sep 17 00:00:00 2001 From: smartalecH Date: Thu, 9 May 2019 15:44:39 -0600 Subject: [PATCH 3/8] refactored method with test and updated materials --- python/geom.py | 65 +++++++++++++++++------------- python/materials.py | 14 ++++++- python/tests/medium_evaluations.py | 46 +++++++++++++++++++++ 3 files changed, 95 insertions(+), 30 deletions(-) create mode 100644 python/tests/medium_evaluations.py diff --git a/python/geom.py b/python/geom.py index e5cb63071..b94001f3e 100755 --- a/python/geom.py +++ b/python/geom.py @@ -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, @@ -211,8 +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 = Vector3() + 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): @@ -237,31 +240,37 @@ def transform(self, m): 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 not isinstance(freq, list): - raise ValueError('Frequencies must be supplied as a list') + if type(freq) is not np.ndarray: + freqs = np.array(freq) + else: + freqs = np.squeeze(freq) + + numFreqs = freqs.size # Initialize with instantaneous dielectric tensor - eps = [Matrix(mp.Vector3(self.epsilon_diag.x, self.epsilon_offdiag.x, self.epsilon_offdiag.y), - mp.Vector3(self.epsilon_offdiag.x, self.epsilon_diag.y, self.epsilon_offdiag.z), - mp.Vector3(self.epsilon_offdiag.y, self.epsilon_offdiag.z, self.epsilon_diag.z)) for i in range(len(freq))] - - # Iterate through frequencies - for i_freq in range(len(freq)): - # Iterate through susceptibilities - for i_sus in range(len(self.E_susceptibilities)): - eps[i_freq] = eps[i_freq] + self.E_susceptibilities[i_sus].eval_susceptibility(freq[i_freq]) - # Account for conductivity term - D_conductivity = Matrix(mp.Vector3(self.D_conductivity_diag.x, self.D_conductivity_offdiag.x, self.D_conductivity_offdiag.y), - mp.Vector3(self.D_conductivity_offdiag.x, self.D_conductivity_diag.y, self.D_conductivity_offdiag.z), - mp.Vector3(self.D_conductivity_offdiag.y, self.D_conductivity_offdiag.z, self.D_conductivity_diag.z)) - ones_matrix = Matrix(mp.Vector3(1, 1, 1), - mp.Vector3(1, 1, 1), - mp.Vector3(1, 1, 1)) - eps[i_freq] = (ones_matrix + 1j/freq[i_freq] * D_conductivity) * eps[i_freq] + 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)) + + # Iterate through susceptibilities + for i_sus in range(len(susceptibilities)): + epsmu = epsmu + susceptibilities[i_sus].eval_susceptibility(freqs) - return eps + # 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 + + # Convert list matrix to 3D numpy array size [freqs,3,3] + return np.array(epsmu) class Susceptibility(object): @@ -288,11 +297,11 @@ def __init__(self, frequency=0.0, gamma=0.0, **kwargs): def eval_susceptibility(self,freq): - sigma = Matrix(mp.Vector3(self.sigma_diag.x, self.sigma_offdiag.x, self.sigma_offdiag.y), + 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)) + mp.Vector3(self.sigma_offdiag.y, self.sigma_offdiag.z, self.sigma_diag.z))),(freq.size,1,1)) - return self.frequency*self.frequency / (self.frequency*self.frequency - freq*freq - 1j*self.gamma*freq) * sigma + 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): @@ -304,11 +313,11 @@ def __init__(self, frequency=0.0, gamma=0.0, **kwargs): def eval_susceptibility(self,freq): - sigma = Matrix(mp.Vector3(self.sigma_diag.x, self.sigma_offdiag.x, self.sigma_offdiag.y), + 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)) + mp.Vector3(self.sigma_offdiag.y, self.sigma_offdiag.z, self.sigma_diag.z))),(freq.size,1,1)) - return -self.frequency*self.frequency / (freq*(freq + 1j*self.gamma)) * sigma + return np.tile(-self.frequency*self.frequency / (freq*(freq + 1j*self.gamma)),(3,3,1)).T * sigma class NoisyLorentzianSusceptibility(LorentzianSusceptibility): diff --git a/python/materials.py b/python/materials.py index c02af9233..d0b595a21 100644 --- a/python/materials.py +++ b/python/materials.py @@ -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 +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) @@ -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 +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) diff --git a/python/tests/medium_evaluations.py b/python/tests/medium_evaluations.py new file mode 100644 index 000000000..a40e5b3e8 --- /dev/null +++ b/python/tests/medium_evaluations.py @@ -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() From b633fa22a5566a0f7c46a5f1ee1501a8e48854c4 Mon Sep 17 00:00:00 2001 From: smartalecH Date: Thu, 9 May 2019 21:44:09 -0600 Subject: [PATCH 4/8] refactored matrix init and streamlined computation with broadcasting --- python/geom.py | 38 ++++++++++++++------------------------ 1 file changed, 14 insertions(+), 24 deletions(-) diff --git a/python/geom.py b/python/geom.py index b94001f3e..82adb5df1 100755 --- a/python/geom.py +++ b/python/geom.py @@ -252,22 +252,18 @@ def _get_epsmu(self, diag, offdiag, susceptibilities, conductivity_diag, conduct else: freqs = np.squeeze(freq) - numFreqs = freqs.size + freqs = freqs[:,np.newaxis,np.newaxis] # 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)) + epsmu = np.expand_dims(Matrix(diag=diag,offdiag=offdiag),axis=0) # 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 + conductivity = np.expand_dims(Matrix(diag=conductivity_diag,offdiag=conductivity_offdiag),axis=0) + epsmu = (1 + 1j/freqs * conductivity) * epsmu # Convert list matrix to 3D numpy array size [freqs,3,3] return np.array(epsmu) @@ -280,9 +276,7 @@ def __init__(self, sigma_diag=Vector3(), sigma_offdiag=Vector3(), sigma=None): self.sigma_offdiag = sigma_offdiag def transform(self, m): - sigma = 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)) + sigma = Matrix(diag=self.sigma_diag,offdiag=self.sigma_offdiag) new_sigma = (m * sigma * m.transpose()) / abs(m.determinant()) self.sigma_diag = mp.Vector3(new_sigma.c1.x, new_sigma.c2.y, new_sigma.c3.z) self.sigma_offdiag = mp.Vector3(new_sigma.c2.x, new_sigma.c3.x, new_sigma.c3.y) @@ -296,12 +290,8 @@ def __init__(self, frequency=0.0, gamma=0.0, **kwargs): 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 + sigma = np.expand_dims(Matrix(diag=self.sigma_diag,offdiag=self.sigma_offdiag),axis=0) + return self.frequency*self.frequency / (self.frequency*self.frequency - freq*freq - 1j*self.gamma*freq) * sigma class DrudeSusceptibility(Susceptibility): @@ -312,12 +302,8 @@ def __init__(self, frequency=0.0, gamma=0.0, **kwargs): 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 + sigma = np.expand_dims(Matrix(diag=self.sigma_diag,offdiag=self.sigma_offdiag),axis=0) + return -self.frequency*self.frequency / (freq*(freq + 1j*self.gamma)) * sigma class NoisyLorentzianSusceptibility(LorentzianSusceptibility): @@ -489,10 +475,14 @@ def __init__(self, vertices, height, axis=Vector3(z=1), center=None, **kwargs): class Matrix(object): - def __init__(self, c1=Vector3(), c2=Vector3(), c3=Vector3()): + def __init__(self, c1=Vector3(), c2=Vector3(), c3=Vector3(), diag=Vector3(), offdiag=Vector3()): self.c1 = c1 self.c2 = c2 self.c3 = c3 + if c1 == c2 == c3 == Vector3(): + self.c1 = Vector3(diag.x,offdiag.x,offdiag.y) + self.c2 = Vector3(offdiag.x,diag.y,offdiag.z) + self.c3 = Vector3(offdiag.y,offdiag.z,diag.z) def __getitem__(self, i): return self.row(i) From 1cd1e115f3c8e71a094ec582a5f206415a52578b Mon Sep 17 00:00:00 2001 From: smartalecH Date: Thu, 9 May 2019 21:52:05 -0600 Subject: [PATCH 5/8] squeeze output for scalar freq inputs --- python/geom.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/geom.py b/python/geom.py index 82adb5df1..7d08e5a5b 100755 --- a/python/geom.py +++ b/python/geom.py @@ -266,7 +266,7 @@ def _get_epsmu(self, diag, offdiag, susceptibilities, conductivity_diag, conduct epsmu = (1 + 1j/freqs * conductivity) * epsmu # Convert list matrix to 3D numpy array size [freqs,3,3] - return np.array(epsmu) + return np.squeeze(epsmu) class Susceptibility(object): From 371f9ccb9047ef37ed25c22b4ed798f7bd0f0bb8 Mon Sep 17 00:00:00 2001 From: smartalecH Date: Fri, 10 May 2019 07:10:31 -0600 Subject: [PATCH 6/8] Add test to makefile --- python/Makefile.am | 1 + 1 file changed, 1 insertion(+) diff --git a/python/Makefile.am b/python/Makefile.am index fd1197209..521f6b7da 100644 --- a/python/Makefile.am +++ b/python/Makefile.am @@ -53,6 +53,7 @@ TESTS = \ $(KDOM_TEST) \ $(TEST_DIR)/ldos.py \ $(MDPYTEST) \ + $(TEST_DIR)/medium_evaluations.py \ $(MPBPYTEST) \ $(MODE_COEFFS_TEST) \ $(MODE_DECOMPOSITION_TEST) \ From a7dd4cd71cd380ec5a94d7e97e764c986b1963e7 Mon Sep 17 00:00:00 2001 From: smartalecH Date: Fri, 10 May 2019 14:57:47 -0600 Subject: [PATCH 7/8] complex hermitian Matrix init --- python/geom.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/geom.py b/python/geom.py index 7d08e5a5b..2829195a6 100755 --- a/python/geom.py +++ b/python/geom.py @@ -481,8 +481,8 @@ def __init__(self, c1=Vector3(), c2=Vector3(), c3=Vector3(), diag=Vector3(), off self.c3 = c3 if c1 == c2 == c3 == Vector3(): self.c1 = Vector3(diag.x,offdiag.x,offdiag.y) - self.c2 = Vector3(offdiag.x,diag.y,offdiag.z) - self.c3 = Vector3(offdiag.y,offdiag.z,diag.z) + self.c2 = Vector3(np.conj(offdiag.x),diag.y,offdiag.z) + self.c3 = Vector3(np.conj(offdiag.y),np.conj(offdiag.z),diag.z) def __getitem__(self, i): return self.row(i) From 373795c314996b18d9bfbaa6695220918442230a Mon Sep 17 00:00:00 2001 From: smartalecH Date: Tue, 14 May 2019 13:19:55 -0600 Subject: [PATCH 8/8] docs --- doc/docs/Python_User_Interface.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index 88b0050e6..2d5191fd0 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -300,6 +300,14 @@ List of dispersive susceptibilities (see below) added to the permeability μ in — Transforms `epsilon`, `mu`, and `sigma` of any [susceptibilities](#susceptibility) by the 3×3 matrix `M`. If `M` is a [rotation matrix](https://en.wikipedia.org/wiki/Rotation_matrix), then the principal axes of the susceptibilities are rotated by `M`. More generally, the susceptibilities χ are transformed to MχMᵀ/|det M|, which corresponds to [transformation optics](http://math.mit.edu/~stevenj/18.369/coordinate-transform.pdf) for an arbitrary curvilinear coordinate transformation with Jacobian matrix M. The absolute value of the determinant is to prevent inadvertent construction of left-handed materials, which are [problematic in nondispersive media](FAQ.md#why-does-my-simulation-diverge-if-0). +**`epsilon(freq` [ scalar, list, or `numpy` array, ]`)`** +- +Calculates the medium's permittivity tensor as a 3x3 `numpy` array at the specified frequency, `freq`. Either a scalar, list, or `numpy` array of frequencies may be supplied. In the case `N` frequency points, a `numpy` array size Nx3x3 is returned. + +**`mu(freq` [ scalar, list, or `numpy` array, ]`)`** +- +Calculates the medium's permeability tensor as a 3x3 `numpy` array at the specified frequency, `freq`. Either a scalar, list, or `numpy` array of frequencies may be supplied. In the case `N` frequency points, a `numpy` array size Nx3x3 is returned. + **material functions** Any function that accepts a `Medium` instance can also accept a user defined Python function. This allows you to specify the material as an arbitrary function of position. The function must have one argument, the position `Vector3`, and return the material at that point, which should be a Python `Medium` instance. This is accomplished by passing a function to the `material_function` keyword argument in the `Simulation` constructor, or the `material` keyword argument in any `GeometricObject` constructor.