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

Jump to bottom

revamp test_refl_angular.py to validate results against theory rather than hard-coded values #2017

Merged
merged 1 commit into from
Apr 7, 2022
Merged

revamp test_refl_angular.py to validate results against theory rather than hard-coded values #2017

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
15 changes: 9 additions & 6 deletions doc/docs/Python_Tutorials/Basics.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -433,7 +433,7 @@ A 1d cell must be used since a higher-dimensional cell will introduce [artificia

Creating an oblique planewave source typically requires specifying two parameters: (1) for periodic structures, the Bloch-periodic wavevector $\vec{k}$ via [`k_point`](../FAQ.md#how-does-k_point-define-the-phase-relation-between-adjacent-unit-cells), and (2) the source amplitude function `amp_func` for setting the $e^{i\vec{k} \cdot \vec{r}}$ spatial dependence ($\vec{r}$ is the position vector). Since we have a 1d cell and the source is at a single point, it is not necessary to specify the source amplitude (see this [2d example](https://github.com/NanoComp/meep/blob/master/python/examples/pw-source.py) for how this is done). The magnitude of the Bloch-periodic wavevector is specified according to the dispersion relation formula for a planewave in homogeneous media with index $n$: $\omega=c|\vec{k}|/n$. As the source in this example is incident from air, $|\vec{k}|$ is simply equal to the frequency $\omega$. Note that specifying $\vec{k}$ corresponds to a single frequency. Any broadband source is therefore incident at a specified angle for only a *single* frequency. This is described in more detail in Section 4.5 ("Efficient Frequency-Angle Coverage") in [Chapter 4](https://arxiv.org/abs/1301.5366) ("Electromagnetic Wave Source Conditions") of the book [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707). In this example, $\omega$ is set to the minimum frequency of the pulse to produce propagating fields at all pulse frequencies. In general, any pulse frequencies which are *less* than any non-zero component of $\vec{k}$ will result in *evanescent* fields (which carry zero power to the far field, and computationally will yield exponentially small power).

In this example, the plane of incidence which contains $\vec{k}$ and the surface normal vector is chosen to be $xz$. The source angle $\theta$ is defined in degrees in the counterclockwise (CCW) direction around the $y$ axis with 0 degrees along the $+z$ axis. In Meep, a 1d cell is defined along the $z$ direction. When $\vec{k}$ is not set, only the $E_x$ and $H_y$ field components are permitted. A non-zero $\vec{k}$ results in a 3d simulation where all field components are included and are complex valued (note that the fields are real, by default). A current source with $E_x$ polarization lies within the plane of incidence and corresponds to the convention of $\mathcal{P}$-polarization. In order to model the $\mathcal{S}$-polarization, we must use an $E_y$ source. This example involves just the $\mathcal{P}$-polarization.
In this example, the plane of incidence which contains $\vec{k}$ and the surface normal vector is chosen to be $xz$. The source angle $\theta$ is defined in degrees in the counterclockwise (CCW) direction around the $y$ axis with 0 degrees along the $+z$ axis. In Meep, a 1d cell is defined along the $z$ direction. When $\vec{k}$ is not set, only the $E_x$ and $H_y$ field components are permitted. A non-zero $\vec{k}$ results in a 3d simulation where all field components are included and are complex valued (note that the fields are real, by default). A current source with $E_x$ polarization lies within the plane of incidence and corresponds to the convention of $\mathcal{P}$ polarization. In order to model the $\mathcal{S}$ polarization, we must use an $E_y$ source. This example involves just the $\mathcal{P}$ polarization.

The simulation script is [examples/refl-angular.py](https://github.com/NanoComp/meep/blob/master/python/examples/refl-angular.py). The notebook is [examples/refl-angular.ipynb](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/refl-angular.ipynb).

Expand All @@ -442,13 +442,12 @@ import meep as mp
import argparse
import math

def main(args):

def main(args):
resolution = args.res

dpml = 1.0 # PML thickness
sz = 10 # size of cell (without PMLs)
sz = 10 + 2*dpml
sz = 10 + 2*dpml # size of cell (without PMLs)
cell_size = mp.Vector3(0,0,sz)
pml_layers = [mp.PML(dpml)]

Expand All @@ -472,7 +471,9 @@ def main(args):
else:
dimensions = 3

sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), component=mp.Ex, center=mp.Vector3(0,0,-0.5*sz+dpml))]
sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
component=mp.Ex,
center=mp.Vector3(0,0,-0.5*sz+dpml))]

sim = mp.Simulation(cell_size=cell_size,
boundary_layers=pml_layers,
Expand All @@ -491,7 +492,9 @@ def main(args):
sim.reset_meep()

# add a block with n=3.5 for the air-dielectric interface
geometry = [mp.Block(mp.Vector3(mp.inf,mp.inf,0.5*sz), center=mp.Vector3(0,0,0.25*sz), material=mp.Medium(index=3.5))]
geometry = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,0.5*sz),
center=mp.Vector3(0,0,0.25*sz),
material=mp.Medium(index=3.5))]

sim = mp.Simulation(cell_size=cell_size,
geometry=geometry,
Expand Down
22 changes: 13 additions & 9 deletions python/examples/refl-angular.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,12 @@
import argparse
import math


def main(args):
resolution = args.res

dpml = 1.0 # PML thickness
sz = 10 # size of computational cell (without PMLs)
sz = 10 + 2*dpml
sz = 10 + 2*dpml # size of computational cell (without PMLs)
cell_size = mp.Vector3(0,0,sz)
pml_layers = [mp.PML(dpml)]

Expand All @@ -18,20 +18,22 @@ def main(args):
fcen = 0.5*(fmin+fmax) # center frequency
df = fmax-fmin # frequency width
nfreq = 50 # number of frequency bins

# rotation angle (in degrees) of source: CCW around Y axis, 0 degrees along +Z axis
theta_r = math.radians(args.theta)

# plane of incidence is XZ
k = mp.Vector3(math.sin(theta_r),0,math.cos(theta_r)).scale(fmin)

# if normal incidence, force number of dimensions to be 1
if theta_r == 0:
dimensions = 1
else:
dimensions = 3

sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), component=mp.Ex, center=mp.Vector3(0,0,-0.5*sz+dpml))]

sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
component=mp.Ex,
center=mp.Vector3(0,0,-0.5*sz+dpml))]

sim = mp.Simulation(cell_size=cell_size,
boundary_layers=pml_layers,
Expand All @@ -42,15 +44,17 @@ def main(args):

refl_fr = mp.FluxRegion(center=mp.Vector3(0,0,-0.25*sz))
refl = sim.add_flux(fcen, df, nfreq, refl_fr)

sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(0,0,-0.5*sz+dpml), 1e-9))

empty_flux = mp.get_fluxes(refl)
empty_data = sim.get_flux_data(refl)
sim.reset_meep()

# add a block with n=3.5 for the air-dielectric interface
geometry = [mp.Block(mp.Vector3(mp.inf,mp.inf,0.5*sz), center=mp.Vector3(0,0,0.25*sz), material=mp.Medium(index=3.5))]
geometry = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,0.5*sz),
center=mp.Vector3(0,0,0.25*sz),
material=mp.Medium(index=3.5))]

sim = mp.Simulation(cell_size=cell_size,
geometry=geometry,
Expand All @@ -70,7 +74,7 @@ def main(args):

for i in range(nfreq):
print("refl:, {}, {}, {}, {}".format(k.x,1/freqs[i],math.degrees(math.asin(k.x/freqs[i])),-refl_flux[i]/empty_flux[i]))

if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('-res', type=int, default=200, help='resolution (default: 200 pixels/um)')
Expand Down
195 changes: 101 additions & 94 deletions python/tests/test_refl_angular.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,127 +1,134 @@
import meep as mp
from utils import ApproxComparisonTestCase
import math
# -*- coding: utf-8 -*-

import unittest
import parameterized
import numpy as np
import math
import meep as mp
from utils import ApproxComparisonTestCase


class TestReflAngular(ApproxComparisonTestCase):
@classmethod
def setUpClass(cls):
cls.resolution = 400 # pixels/μm

def test_refl_angular(self):
resolution = 100
cls.n1 = 1.4 # refractive index of medium 1
cls.n2 = 3.5 # refractive index of medium 2

dpml = 1.0
sz = 10
sz = sz + 2 * dpml
cell_size = mp.Vector3(z=sz)
pml_layers = [mp.PML(dpml)]
cls.dpml = 1.0
cls.dz = 7.0
cls.sz = cls.dz+2*cls.dpml

wvl_min = 0.4
wvl_max = 0.8
fmin = 1 / wvl_max
fmax = 1 / wvl_min
fcen = 0.5 * (fmin + fmax)
df = fmax - fmin
nfreq = 50
cls.wvl_min = 0.4
cls.wvl_max = 0.8
cls.fmin = 1/cls.wvl_max
cls.fmax = 1/cls.wvl_min
cls.fcen = 0.5*(cls.fmin+cls.fmax)
cls.df = cls.fmax-cls.fmin
cls.nfreq = 11

theta_r = math.radians(0)
def refl_angular(self, theta):
theta_r = math.radians(theta)

k = mp.Vector3(math.sin(theta_r), 0, math.cos(theta_r)).scale(fcen)
# wavevector (in source medium); plane of incidence is XZ
k = mp.Vector3(0,0,1).rotate(mp.Vector3(0,1,0),theta_r).scale(self.n1*self.fmin)

dimensions = 1
if theta == 0:
dimensions = 1
else:
dimensions = 3

sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ex,
center=mp.Vector3(z=-0.5 * sz + dpml))]
cell_size = mp.Vector3(z=self.sz)
pml_layers = [mp.PML(self.dpml)]

sim = mp.Simulation(cell_size=cell_size,
boundary_layers=pml_layers,
sources=sources,
k_point=k,
sources = [mp.Source(mp.GaussianSource(self.fcen, fwidth=self.df),
component=mp.Ex, # P polarization
center=mp.Vector3(z=-0.5*self.sz+self.dpml))]

sim = mp.Simulation(resolution=self.resolution,
cell_size=cell_size,
dimensions=dimensions,
resolution=resolution)
default_material=mp.Medium(index=self.n1),
sources=sources,
boundary_layers=pml_layers,
k_point=k)

refl_fr = mp.FluxRegion(center=mp.Vector3(z=-0.25 * sz))
refl = sim.add_flux(fcen, df, nfreq, refl_fr, decimation_factor=1)
mon_pt = -0.5*self.sz+self.dpml+0.25*self.dz
refl_fr = mp.FluxRegion(center=mp.Vector3(z=mon_pt))
refl = sim.add_flux(self.fcen,
self.df,
self.nfreq,
refl_fr)

sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(z=-0.5 * sz + dpml), 1e-9))
termination_cond = mp.stop_when_fields_decayed(50,
mp.Ex,
mp.Vector3(z=mon_pt),
1e-9)
sim.run(until_after_sources=termination_cond)

empty_data = sim.get_flux_data(refl)
empty_flux = mp.get_fluxes(refl)
sim.reset_meep()

geometry = [mp.Block(mp.Vector3(mp.inf, mp.inf, 0.5 * sz), center=mp.Vector3(z=0.25 * sz),
material=mp.Medium(index=3.5))]
geometry = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,0.5*self.sz),
center=mp.Vector3(z=0.25*self.sz),
material=mp.Medium(index=self.n2))]

sim = mp.Simulation(cell_size=cell_size,
geometry=geometry,
boundary_layers=pml_layers,
sim = mp.Simulation(resolution=self.resolution,
cell_size=cell_size,
dimensions=dimensions,
default_material=mp.Medium(index=self.n1),
sources=sources,
boundary_layers=pml_layers,
k_point=k,
dimensions=dimensions,
resolution=resolution)
geometry=geometry)

refl = sim.add_flux(fcen, df, nfreq, refl_fr, decimation_factor=1)
refl = sim.add_flux(self.fcen,
self.df,
self.nfreq,
refl_fr)
sim.load_minus_flux_data(refl, empty_data)

sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(z=-0.5 * sz + dpml), 1e-9))
sim.run(until_after_sources=termination_cond)

refl_flux = mp.get_fluxes(refl)
freqs = mp.get_flux_freqs(refl)

expected = [
(1.25, -1.123696883299492e-6),
(1.2755102040816326, -2.5749667658387866e-6),
(1.3010204081632653, -5.70480204599006e-6),
(1.3265306122448979, -1.2220464827582253e-5),
(1.3520408163265305, -2.531247480206961e-5),
(1.3775510204081631, -5.069850309492639e-5),
(1.4030612244897958, -9.819256552437341e-5),
(1.4285714285714284, -1.8390448659017395e-4),
(1.454081632653061, -3.330762066794769e-4),
(1.4795918367346936, -5.833650417163753e-4),
(1.5051020408163263, -9.8807834237052e-4),
(1.5306122448979589, -0.001618472171445976),
(1.5561224489795915, -0.0025638388059825985),
(1.5816326530612241, -0.003927863989816029),
(1.6071428571428568, -0.005819831283556752),
(1.6326530612244894, -0.008339881000982728),
(1.658163265306122, -0.011558769654206626),
(1.6836734693877546, -0.015494308354153143),
(1.7091836734693873, -0.02008850084337135),
(1.73469387755102, -0.025190871516857616),
(1.7602040816326525, -0.030553756123198477),
(1.7857142857142851, -0.03584404966066722),
(1.8112244897959178, -0.040672967700428275),
(1.8367346938775504, -0.04464118393086191),
(1.862244897959183, -0.047392712128477496),
(1.8877551020408156, -0.048667403362887635),
(1.9132653061224483, -0.048341494285878264),
(1.938775510204081, -0.04644739000778679),
(1.9642857142857135, -0.043168390293742316),
(1.9897959183673462, -0.0388094755730579),
(2.0153061224489788, -0.03375052221907117),
(2.0408163265306114, -0.02839209067703472),
(2.066326530612244, -0.023104245646230648),
(2.0918367346938767, -0.01818725699718267),
(2.1173469387755093, -0.013849270759480073),
(2.142857142857142, -0.010201733597436358),
(2.1683673469387745, -0.007269616609175294),
(2.193877551020407, -0.005011210495189995),
(2.21938775510204, -0.0033417192031464896),
(2.2448979591836724, -0.0021557351734376256),
(2.270408163265305, -0.0013453062176115673),
(2.2959183673469377, -8.121742663131631e-4),
(2.3214285714285703, -4.7433135191915683e-4),
(2.346938775510203, -2.6799188013374266e-4),
(2.3724489795918355, -1.464781343401766e-4),
(2.397959183673468, -7.745339273024636e-5),
(2.423469387755101, -3.9621374769542025e-5),
(2.4489795918367334, -1.9608458558430508e-5),
(2.474489795918366, -9.38818477949983e-6),
(2.4999999999999987, -4.3484671364929225e-6),
]

tol = 1e-7 if mp.is_single_precision() else 1e-8
self.assertClose(expected, list(zip(freqs, refl_flux)), epsilon=tol)
Rs = -np.array(refl_flux)/np.array(empty_flux)

thetas = []
for i in range(self.nfreq):
thetas.append(math.asin(k.x/(self.n1*freqs[i])))

return freqs, thetas, Rs

@parameterized.parameterized.expand([(0,), (20.6,)])
def test_refl_angular(self,theta):
fmeep, tmeep, Rmeep = self.refl_angular(theta)

# angle of refracted planewave in medium n2 for an
# incident planewave in medium n1 at angle theta_in
theta_out = lambda theta_in: math.asin(self.n1*math.sin(theta_in)/self.n2)

# Fresnel reflectance for P polarization in medium n2 for
# an incident planewave in medium n1 at angle theta_in
Rfresnel = lambda theta_in: (math.fabs((self.n1*math.cos(theta_out(theta_in))-self.n2*math.cos(theta_in)) /
(self.n1*math.cos(theta_out(theta_in))+self.n2*math.cos(theta_in)))**2)

Ranalytic = np.empty((self.nfreq,))
print("refl:, wavelength (μm), incident angle (°), reflectance (Meep), reflectance (analytic), error")
for i in range(self.nfreq):
Ranalytic[i] = Rfresnel(tmeep[i])
err = abs(Rmeep[i]-Ranalytic[i])/Ranalytic[i]
print("refl:, {:4.2f}, {:4.2f}, {:8.6f}, {:8.6f}, {:6.4f}".format(1/fmeep[i],
math.degrees(tmeep[i]),
Rmeep[i],
Ranalytic[i],
err))

tol = 0.004
self.assertClose(Rmeep, Ranalytic, epsilon=tol)


if __name__ == '__main__':
Expand Down