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

Add unit test for broadband sources at fixed angle in a periodic cell (BFAST) #2722

Merged
merged 3 commits into from
Nov 25, 2023
Merged

Add unit test for broadband sources at fixed angle in a periodic cell (BFAST) #2722

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5122048
add unit test for broadband sources at fixed angle (BFAST)
oskooi Nov 20, 2023
ee61528
include index of source medium in bfast_k_bar
oskooi Nov 20, 2023
685d9bd
set Courant factor using analytic expression involving bfast_k_bar
oskooi Nov 25, 2023
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
213 changes: 140 additions & 73 deletions python/tests/test_refl_angular.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
import math
from typing import List, Tuple
import unittest

import numpy as np
Expand All @@ -8,45 +9,75 @@
import meep as mp


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

cls.n1 = 1.4 # refractive index of medium 1
cls.n2 = 3.5 # refractive index of medium 2

cls.dpml = 1.0
cls.dz = 7.0
cls.sz = cls.dz + 2 * cls.dpml

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

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

# 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)
)
cls.t_pml = 1.0
cls.length_z = 7.0
cls.size_z = cls.length_z + 2 * cls.t_pml

cls.wavelength_min = 0.4
cls.wavelength_max = 0.8
cls.frequency_min = 1 / cls.wavelength_max
cls.frequency_max = 1 / cls.wavelength_min
cls.frequency_center = 0.5 * (cls.frequency_min + cls.frequency_max)
cls.frequency_width = cls.frequency_max - cls.frequency_min
cls.num_freq = 11

def reflectance_angular(
self, theta_deg: float, need_bfast: bool
) -> Tuple[List, List, np.ndarray]:
"""Computes properties of the incident and reflected planewave.

Args:
theta_deg: angle of incident planewave.
need_bfast: whether to use the same angle for the incident planewave
for all frequencies. If False, the incident angle is frequency
dependent.

Returns:
A 3-tuple comprising the frequencies of the incident planewave,
angles of the incident planewave, and the reflectance.
"""
theta_rad = math.radians(theta_deg)

if need_bfast:
bfast_k_bar = (self.n1 * np.sin(theta_rad), 0, 0)

Courant = (1 - bfast_k_bar[0]) / 3**0.5

k = mp.Vector3()
else:
bfast_k_bar = (0, 0, 0)

Courant = 0.5

# Wavevector in source medium with refractive index n1.
# Plane of incidence is XZ. Rotation is counter clockwise about
# Y axis. A rotation angle of zero is the +Z axis.
k = (
mp.Vector3(0, 0, 1)
.rotate(mp.Vector3(0, 1, 0), theta_rad)
.scale(self.n1 * self.frequency_min)
)

dimensions = 1 if theta == 0 else 3
cell_size = mp.Vector3(z=self.sz)
pml_layers = [mp.PML(self.dpml)]
dimensions = 1 if theta_deg == 0 else 3
cell_size = mp.Vector3(z=self.size_z)
pml_layers = [mp.PML(self.t_pml)]

# P polarization.
source_component = mp.Ex

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),
mp.GaussianSource(self.frequency_center, fwidth=self.frequency_width),
component=source_component,
center=mp.Vector3(z=-0.5 * self.size_z + self.t_pml),
)
]

Expand All @@ -58,25 +89,31 @@ def refl_angular(self, theta):
sources=sources,
boundary_layers=pml_layers,
k_point=k,
need_bfast=need_bfast,
bfast_k_bar=bfast_k_bar,
Courant=Courant,
)

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)
monitor_point = -0.5 * self.size_z + self.t_pml + 0.25 * self.length_z
monitor_region = mp.FluxRegion(center=mp.Vector3(z=monitor_point))
flux_monitor = sim.add_flux(
self.frequency_center, self.frequency_width, self.num_freq, monitor_region
)

termination_cond = mp.stop_when_fields_decayed(
50, mp.Ex, mp.Vector3(z=mon_pt), 1e-9
termination_criteria = mp.stop_when_fields_decayed(
50, source_component, mp.Vector3(z=monitor_point), 1e-6
)
sim.run(until_after_sources=termination_cond)
sim.run(until_after_sources=termination_criteria)

empty_data = sim.get_flux_data(flux_monitor)
empty_flux = mp.get_fluxes(flux_monitor)

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

geometry = [
mp.Block(
size=mp.Vector3(mp.inf, mp.inf, 0.5 * self.sz),
center=mp.Vector3(z=0.25 * self.sz),
size=mp.Vector3(mp.inf, mp.inf, 0.5 * self.size_z),
center=mp.Vector3(z=0.25 * self.size_z),
material=mp.Medium(index=self.n2),
)
]
Expand All @@ -89,58 +126,88 @@ def refl_angular(self, theta):
sources=sources,
boundary_layers=pml_layers,
k_point=k,
need_bfast=need_bfast,
bfast_k_bar=bfast_k_bar,
Courant=Courant,
geometry=geometry,
)

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=termination_cond)

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

Rs = -np.array(refl_flux) / np.array(empty_flux)

thetas = [math.asin(k.x / (self.n1 * freqs[i])) for i in range(self.nfreq)]
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)
flux_monitor = sim.add_flux(
self.frequency_center, self.frequency_width, self.num_freq, monitor_region
)
sim.load_minus_flux_data(flux_monitor, empty_data)

sim.run(until_after_sources=termination_criteria)

flux_monitor_flux = mp.get_fluxes(flux_monitor)
freqs = mp.get_flux_freqs(flux_monitor)

reflectance = -np.array(flux_monitor_flux) / np.array(empty_flux)

if need_bfast:
theta_in_rad = [theta_rad] * self.num_freq
else:
# Returns the angle of the incident planewave in medium n1 based
# on its frequency given a fixed wavevector component in X.
theta_in_rad = [
math.asin(k.x / (self.n1 * freqs[i])) for i in range(self.num_freq)
]

return freqs, theta_in_rad, reflectance

@parameterized.parameterized.expand([(0, False), (20.6, False), (35.7, True)])
def test_reflectance_angular(self, theta_deg: float, need_bfast: bool):
(
frequency_meep,
theta_in_rad_meep,
reflectance_meep,
) = self.reflectance_angular(theta_deg, need_bfast)

# Returns angle of refracted planewave in medium n2 given
# an incident planewave in medium n1 at angle theta_in_rad.
theta_out = lambda theta_in_rad: math.asin(
self.n1 * math.sin(theta_in_rad) / 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: (
# Returns Fresnel reflectance for P polarization in medium n2
# for an incident planewave in medium n1 at angle theta_in_rad.
reflectance_fresnel = lambda theta_in_rad: (
math.fabs(
(self.n1 * math.cos(theta_out(theta_in)) - self.n2 * math.cos(theta_in))
(
self.n1 * math.cos(theta_out(theta_in_rad))
- self.n2 * math.cos(theta_in_rad)
)
/ (
self.n1 * math.cos(theta_out(theta_in))
+ self.n2 * math.cos(theta_in)
self.n1 * math.cos(theta_out(theta_in_rad))
+ self.n2 * math.cos(theta_in_rad)
)
)
** 2
)

Ranalytic = np.empty((self.nfreq,))
reflectance_analytic = np.empty((self.num_freq,))
print(
"refl:, wavelength (μm), incident angle (°), reflectance (Meep), reflectance (analytic), error"
"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]
for i in range(self.num_freq):
reflectance_analytic[i] = reflectance_fresnel(theta_in_rad_meep[i])
err = (
abs(reflectance_meep[i] - reflectance_analytic[i])
/ reflectance_analytic[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
1 / frequency_meep[i],
math.degrees(theta_in_rad_meep[i]),
reflectance_meep[i],
reflectance_analytic[i],
err,
)
)

tol = 0.005 if mp.is_single_precision() else 0.004
self.assertClose(Rmeep, Ranalytic, epsilon=tol)
tol = 0.03
self.assertClose(reflectance_meep, reflectance_analytic, epsilon=tol)


if __name__ == "__main__":
Expand Down