loss rate calculation is off
alexysong opened this issue · 3 comments
I suspect it is either because of multiple layers, or because cladding is not 1. An example:
import numpy as np
import legume
from legume.phc import PhotCryst, Lattice
from legume.shapes import Square, Circle, Poly
from legume.gme import GuidedModeExp
from legume.utils import plot_reciprocal
import matplotlib.pyplot as plt
# Initialize a lattice (can be 'square', 'hexagonal', or defined by primitive vectors)
lattice = Lattice('square')
# Initialize a PhC (optional permittivity in the lower and upper cladding)
phc = PhotCryst(lattice, eps_l=3.2826**2, eps_u=3.2826**2)
phc.add_layer(d=1.007, eps_b=3.5032**2) # PhC layer
# A triangle
triangle = Poly(x_edges=[0., 0.5256, 0.], y_edges=[0., 0., 0.5256], eps=1.)
phc.layers[-1].add_shape(triangle)
phc.add_layer(d=0.375, eps_b=3.5032**2) # sch layer
phc.add_layer(d=0.495, eps_b=3.3955**2) # active layer
# We can plot an overview of what we've built so far
phc.plot_overview()
# Initialize GME and plot what we're simulating
gmax = 5 # plane wave truncation
gme = legume.GuidedModeExp(phc, gmax=gmax)
gme.plot_overview_ft()
n_k = 21
path = phc.lattice.bz_path(['G', np.array([np.pi/2., 0])], [n_k-1])
# See the legume README to understand the options below
neig = 12
options = {'gmode_inds': np.arange(8), 'gmode_npts': 4000, 'numeig': neig, 'verbose': True}
gme.run(kpoints=path.kpoints, **options)
ax = legume.viz.bands(gme)
ks = np.arange(n_k)
freqs_i = []
for kind in ks:
# kind = 0 # index of the k point
minds = np.arange(0, neig) # mode indexes whose loss rates will be computed
(freqs_im, coup_l, coup_u) = gme.compute_rad(kind=kind, minds=minds)
freqs_i.append(freqs_im)
# Note that states below light line still have zero losses
freqs_i = np.array(freqs_i)
qf = gme.freqs/2./freqs_i
# print("Frequencies of first 4 modes at kx = %f are: " % gme.kpoints[0, kind], gme.freqs[kind, minds])
# print("Loss rates of first 4 modes at kx = %f are : " % gme.kpoints[0, kind], freqs_im)
f1 = plt.figure()
ax = f1.add_subplot(111)
ax.plot(np.log10(qf))
Expect: Q factor at gamma to be in the 10^5 range.
Actual outcome: Q factor in 10^2 range.
Update: I think it's probably related to the cladding being not air. If I delete these two lines:
phc.add_layer(d=0.375, eps_b=3.5032**2) # sch layer
phc.add_layer(d=0.495, eps_b=3.3955**2) # active layer
I am left with a one layer structure, with claddings being not air. The calculated loss rate is still off.
So the code above doesn't run at all, because currently loss rates are only implemented for even gmode_inds. So I can't really reproduce what you're seeing.
Maybe when I'm done with the loss rates fully, we can have a look into that again.
P. S. Note that I've never tried GME with such a small contrast between slab and claddings, but I hope that it will still work well - we'll see.