# -*- coding: utf-8 -*- # polarization grating from C. Oh and M.J. Escuti, Optics Letters, Vol. 33, No. 20, pp. 2287-9, 2008 # note: reference uses z as the propagation direction and y as the out-of-plane direction; this script uses x and z, respectively import meep as mp import numpy as np import math import matplotlib.pyplot as plt resolution = 50 # pixels/μm dpml = 1.0 # PML thickness dsub = 1.0 # substrate thickness dpad = 1.0 # padding thickness k_point = mp.Vector3(0,0,0) pml_layers = [mp.PML(thickness=dpml,direction=mp.X)] n_0 = 1.55 delta_n = 0.159 epsilon_diag = mp.Matrix(mp.Vector3(n_0**2,0,0),mp.Vector3(0,n_0**2,0),mp.Vector3(0,0,(n_0+delta_n)**2)) wvl = 0.54 # center wavelength fcen = 1/wvl # center frequency def pol_grating(d,ph,gp,nmode): sx = dpml+dsub+d+d+dpad+dpml sy = gp cell_size = mp.Vector3(sx,sy,0) # twist angle of nematic director; from equation 1b def phi(p): xx = p.x-(-0.5*sx+dpml+dsub) if (xx >= 0) and (xx <= d): return math.pi*p.y/gp + ph*xx/d else: return math.pi*p.y/gp - ph*xx/d + 2*ph # return the anisotropic permittivity tensor for a uniaxial, twisted nematic liquid crystal def lc_mat(p): # rotation matrix for rotation around x axis Rx = mp.Matrix(mp.Vector3(1,0,0),mp.Vector3(0,math.cos(phi(p)),math.sin(phi(p))),mp.Vector3(0,-math.sin(phi(p)),math.cos(phi(p)))) lc_epsilon = Rx * epsilon_diag * Rx.transpose() lc_epsilon_diag = mp.Vector3(lc_epsilon[0].x,lc_epsilon[1].y,lc_epsilon[2].z) lc_epsilon_offdiag = mp.Vector3(lc_epsilon[1].x,lc_epsilon[2].x,lc_epsilon[2].y) return mp.Medium(epsilon_diag=lc_epsilon_diag,epsilon_offdiag=lc_epsilon_offdiag) geometry = [mp.Block(center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub)),size=mp.Vector3(dpml+dsub,mp.inf,mp.inf),material=mp.Medium(index=n_0)), mp.Block(center=mp.Vector3(-0.5*sx+dpml+dsub+d),size=mp.Vector3(2*d,mp.inf,mp.inf),material=lc_mat)] # linear-polarized planewave pulse source src_pt = mp.Vector3(-0.5*sx+dpml+0.3*dsub,0,0) sources = [mp.Source(mp.GaussianSource(fcen,fwidth=0.05*fcen), component=mp.Ez, center=src_pt, size=mp.Vector3(0,sy,0)), mp.Source(mp.GaussianSource(fcen,fwidth=0.05*fcen), component=mp.Ey, center=src_pt, size=mp.Vector3(0,sy,0))] sim = mp.Simulation(resolution=resolution, cell_size=cell_size, boundary_layers=pml_layers, k_point=k_point, sources=sources, default_material=mp.Medium(index=n_0)) tran_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0) tran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0))) sim.run(until_after_sources=100) input_flux = mp.get_fluxes(tran_flux) sim.reset_meep() sim = mp.Simulation(resolution=resolution, cell_size=cell_size, boundary_layers=pml_layers, k_point=k_point, sources=sources, geometry=geometry) tran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0))) sim.run(until_after_sources=300) res1 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.ODD_Z+mp.EVEN_Y) res2 = sim.get_eigenmode_coefficients(tran_flux, range(1,nmode+1), eig_parity=mp.EVEN_Z+mp.ODD_Y) angles = [math.degrees(math.acos(kdom.x/fcen)) for kdom in res1.kdom] return input_flux[0], angles, res1.alpha[:,0,0], res2.alpha[:,0,0]; ph_uniaxial = 0 # chiral layer twist angle for uniaxial grating ph_twisted = 70 # chiral layer twist angle for bilayer grating gp = 6.5 # grating period nmode = 5 # number of mode coefficients to compute dd = np.arange(0.1,3.5,0.1) # chiral layer thickness m0_uniaxial = np.zeros(dd.size) m1_uniaxial = np.zeros(dd.size) ang_uniaxial = np.zeros(dd.size) m0_twisted = np.zeros(dd.size) m1_twisted = np.zeros(dd.size) ang_twisted = np.zeros(dd.size) for k in range(len(dd)): input_flux, angles, coeffs1, coeffs2 = pol_grating(0.5*dd[k],math.radians(ph_uniaxial),gp,nmode) tran = (abs(coeffs1)**2+abs(coeffs2)**2)/input_flux for m in range(nmode): print("tran (uniaxial):, {}, {:.2f}, {:.5f}".format(m,angles[m],tran[m])) m0_uniaxial[k] = tran[0] m1_uniaxial[k] = tran[1] ang_uniaxial[k] = angles[1] input_flux, angles, coeffs1, coeffs2 = pol_grating(dd[k],math.radians(ph_twisted),gp,nmode) tran = (abs(coeffs1)**2+abs(coeffs2)**2)/input_flux for m in range(nmode): print("tran (twisted):, {}, {:.2f}, {:.5f}".format(m,angles[m],tran[m])) m0_twisted[k] = tran[0] m1_twisted[k] = tran[1] ang_twisted[k] = angles[1] cos_angles = [math.cos(math.radians(t)) for t in ang_uniaxial] tran = m0_uniaxial+2*m1_uniaxial eff_m0 = m0_uniaxial/tran eff_m1 = (2*m1_uniaxial/tran)/cos_angles phase = delta_n*dd/wvl eff_m0_analytic = [math.cos(math.pi*p)**2 for p in phase] eff_m1_analytic = [math.sin(math.pi*p)**2 for p in phase] plt.figure(dpi=150) plt.subplot(1,2,1) plt.plot(phase,eff_m0,'bo-',clip_on=False,label='0th order (meep)') plt.plot(phase,eff_m0_analytic,'b--',clip_on=False,label='0th order (analytic)') plt.plot(phase,eff_m1,'ro-',clip_on=False,label='±1 orders (meep)') plt.plot(phase,eff_m1_analytic,'r--',clip_on=False,label='±1 orders (analytic)') plt.axis([0, 1.0, 0, 1]) plt.xticks([t for t in np.arange(0,1.2,0.2)]) plt.xlabel("phase delay Δnd/λ") plt.ylabel("diffraction efficiency @ λ = 0.54 μm") plt.legend(loc='center') plt.title("homogeneous uniaxial grating") cos_angles = [math.cos(math.radians(t)) for t in ang_twisted] tran = m0_twisted+2*m1_twisted eff_m0 = m0_twisted/tran eff_m1 = (2*m1_twisted/tran)/cos_angles plt.subplot(1,2,2) plt.plot(phase,eff_m0,'bo-',clip_on=False,label='0th order (meep)') plt.plot(phase,eff_m1,'ro-',clip_on=False,label='±1 orders (meep)') plt.axis([0, 1.0, 0, 1]) plt.xticks([t for t in np.arange(0,1.2,0.2)]) plt.xlabel("phase delay Δnd/λ") plt.ylabel("diffraction efficiency @ λ = 0.54 μm") plt.legend(loc='center') plt.title("bilayer twisted-nematic grating") plt.show()