1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
|
# -*- coding: utf-8 -*-
import meep as mp
from meep import mpb
import numpy as np
import matplotlib.pyplot as plt
resolution = 128 # pixels/μm
Si = mp.Medium(index=3.45)
syz = 10
geometry_lattice = mp.Lattice(size=mp.Vector3(0,syz,syz))
k_points = [mp.Vector3(0.5)]
num_bands = 1
tolerance = 1e-9
a = 1.0 # waveguide width
def parallel_waveguide(s,yodd):
geometry = [mp.Block(center=mp.Vector3(0,-0.5*(s+a),0),
size=mp.Vector3(mp.inf,a,a),
material=Si),
mp.Block(center=mp.Vector3(0,0.5*(s+a),0),
size=mp.Vector3(mp.inf,a,a),
material=Si)]
ms = mpb.ModeSolver(resolution=resolution,
k_points=k_points,
geometry_lattice=geometry_lattice,
geometry=geometry,
num_bands=num_bands,
tolerance=tolerance)
if yodd:
ms.run_yodd_zodd()
else:
ms.run_yeven_zodd()
f = ms.get_freqs()[0]
vg = ms.compute_group_velocity_component(mp.Vector3(1,0,0))[0]
return f,vg
ss = np.arange(0.1,1.2,0.2)
f_odd = np.zeros(len(ss))
vg_odd = np.zeros(len(ss))
f_even = np.zeros(len(ss))
vg_even = np.zeros(len(ss))
for j in range(len(ss)):
f_odd[j], vg_odd[j] = parallel_waveguide(ss[j],True)
f_even[j], vg_even[j] = parallel_waveguide(ss[j],False)
ds = ss[1]-ss[0]
def compute_force(f,vg):
f_avg = 0.5*(f[:-1]+f[1:])
df = f[1:]-f[:-1]
vg_avg = 0.5*(vg[:-1]+vg[1:])
return np.multiply(np.multiply(-1/f_avg,df/ds), 1/vg_avg)
force_odd = compute_force(f_odd,vg_odd)
force_even = compute_force(f_even,vg_even)
plt.plot(ss[:-1],force_odd,'b-',label='antisymmetric')
plt.plot(ss[:-1],force_even,'r-',label='symmetric')
plt.xlabel("waveguide separation s/a")
plt.ylabel("optical force (F/L)(ac/P)")
plt.legend(loc='upper right')
plt.xticks(np.arange(0,1.2,0.2))
plt.yticks(np.arange(-1.5,1.0,0.5))
plt.show()
|