import numpy as np import matplotlib.pyplot as plt from scipy.linalg import toeplitz, eig def fmm1d_te_layer_modes(eps, period, k0, kx, N): epsf = np.fft.fft(eps)/len(eps) # With N Fourier modes per +- side, we have 2N+1 unknowns. The matrix # that contains the k0^2 eps_j terms has eps_0 on the diagonal and then # we need up to eps_2N. epsf_picked = np.concatenate([ epsf[-2*N:], epsf[0:2*N+1] ]) m = np.arange(2*N+1) - N K = k0**2*toeplitz(epsf_picked[2*N:], r = np.flip(epsf_picked[0:2*N+1])) D = np.diag(-(kx + m*2*np.pi/period)**2) beta2, phie = eig(K+D) # If there are negative imaginary parts, multiply by - beta = np.sqrt(beta2) beta = np.where(beta.real+beta.imag >= 0, beta, -beta) return beta, phie def homogeneous_layer_modes(eps_scalar, period, k0, kx, N): m = np.arange(2*N+1) - N kgr = 2*np.pi/period beta = np.sqrt(k0**2*eps_scalar*(1+0j) - (kx + m*kgr)**2) phie = np.eye(2*N+1) return beta, phie def fourier_to_real(coeffs, x, period): # The fmm1d_te_layer_modes function calculates the modes in a coordinate # system that has a zero on the left-hand side of the unit cell. N, Nmodes = coeffs.shape m = np.arange(N) - N//2 field = np.zeros((len(x), Nmodes), dtype='complex128') for iw in range(N): wave = np.exp(1j*2*np.pi/period*m[iw]*(x-period/2)) for im in range(Nmodes): field[:,im] += coeffs[iw,im]*wave return field def transfermatrix_interface(beta1, modes1, beta2, modes2): B1 = np.diag(beta1) B2 = np.diag(beta2) M1 = np.block([[modes1, modes1], [modes1 @ B1, -modes1 @ B1]]) M2 = np.block([[modes2, modes2], [modes2 @ B2, -modes2 @ B2]]) #return np.linalg.inv(M2) @ M1 return np.linalg.solve(M2,M1) def transfermatrix_propagation(beta, d): prop1 = np.diag(np.exp(1j*beta*d)) prop2 = np.diag(np.exp(-1j*beta*d)) z = np.zeros(prop1.shape) return np.block([[prop1, z], [z, prop2]]) def fmm1d_te(wl0, theta, period, eps_in, eps_out, eps, thick, N): # Should return eta_r, eta_t, r, t # The whole structure is # interface eps_in - eps[0] # propagation eps[0] # interface eps[0] - eps[1] # ... # interface eps[-1] - eps_out # We need the modes in each layer as well as on the input and output sides. # Rows of eps: the permittivity in one layer Ntot = 2*N + 1 Nlayers, Nx = eps.shape k0 = 2*np.pi/wl0 kx = k0*np.sqrt(eps_in)*np.sin(theta) # Input and output side modes beta_in, phi_in = homogeneous_layer_modes(eps_in, period, k0, kx, N) beta_out, phi_out = homogeneous_layer_modes(eps_out, period, k0, kx, N) beta = np.zeros((Nlayers, Ntot), dtype='complex128') phi = np.zeros((Nlayers, Ntot, Ntot), dtype='complex128') for il in range(Nlayers): beta[il,:], phi[il,:,:] = fmm1d_te_layer_modes(eps[il,:], period, k0, kx, N) # The transfer matrix of an interface is # inv([VB VB ; VB bB VB bB]) [VA VA ; VA bA VA bA] # where VA is the mode matrix ("phi") and bA is a diagonal matrix containing # the propagation constants. # The transfer matrix of propagation is # [exp(i beta d) 0 ] # [ 0 exp(-i beta d)] T = transfermatrix_interface(beta_in, phi_in, beta[0,:], phi[0,:,:]) for il in range(Nlayers-1): T = transfermatrix_propagation(beta[il,:], thick[il]) @ T T = transfermatrix_interface(beta[il,:], phi[il,:,:], beta[il+1,:], phi[il+1,:,:]) @ T T = transfermatrix_propagation(beta[-1,:], thick[-1]) @ T T = transfermatrix_interface(beta[-1,:], phi[-1,:,:], beta_out, phi_out) @ T # That's the transfer matrix of the whole structure. With some incident # field, we get # [u_tr] [u_inc] # [0] = T [u_ref] # and we can solve the transmitted and reflected fields # u_tr = T11 u_inc + T12 u_ref # 0 = T21 u_inc + T22 u_ref # yielding # u_ref = -inv(T22) T21 u_inc # u_tr = (T11 - T12 inv(T22) T21) u_inc T11 = T[0:Ntot,0:Ntot] T12 = T[0:Ntot,Ntot:] T21 = T[Ntot:,0:Ntot] T22 = T[Ntot:,Ntot:] i_inc = N u_inc = np.zeros(Ntot, dtype='complex128') u_inc[i_inc] = 1.0 u_ref = -np.linalg.solve(T22,T21) @ u_inc u_tr = (T11 - T12 @ np.linalg.solve(T22,T21)) @ u_inc # Those are the amplitudes of the electric field. # Now to calculate the power efficiency through # Sz = 1/(2 eta) |E|^2 * cos(theta) # Sz / Sz_inc = eta_inc/eta |E|^2 cos(theta) / cos(theta_inc) # = eta_inc/eta |E|^2 kz k_inc / (k kz_inc) # = |E|^2 kz/kz_inc kz_inc = k0*np.sqrt(eps_in)*np.cos(theta) eff_tr = abs(u_tr)**2 * beta_out.real / kz_inc eff_ref = abs(u_ref)**2 * beta_in.real / kz_inc return eff_ref, eff_tr, u_ref, u_tr def smatrix_interface(beta1, modes1, beta2, modes2): B1 = np.diag(beta1) B2 = np.diag(beta2) M1 = np.block([[modes1, -modes2], [modes1 @ B1, modes2 @ B2]]) M2 = np.block([[modes2, -modes1], [modes2 @ B2, modes1 @ B1]]) return np.linalg.solve(M2,M1) def smatrix_propagation(beta, d): prop = np.diag(np.exp(1j*beta*d)) z = np.zeros(prop.shape) return np.block([[prop, z],[z,prop]]) def stack_smatrices(S, Q): N = S.shape[0]//2 S11 = S[0:N,0:N] S12 = S[0:N,N:] S21 = S[N:,0:N] S22 = S[N:,N:] Q11 = Q[0:N,0:N] Q12 = Q[0:N,N:] Q21 = Q[N:,0:N] Q22 = Q[N:,N:] I = np.eye(S11.shape[0]) M = np.linalg.inv(I - S12 @ Q21) revM = np.linalg.inv(I - Q21 @ S12) return np.block([[Q11 @ M @ S11 , Q12 + Q11 @ S12 @ revM @ Q22], [S21 + S22 @ Q21 @ M @ S11 , S22 @ revM @ Q22]]) def fmm1d_te_smatrix(wl0, theta, period, eps_in, eps_out, eps, thick, N): Ntot = 2*N + 1 Nlayers, Nx = eps.shape k0 = 2*np.pi/wl0 kx = k0*np.sqrt(eps_in)*np.sin(theta) beta_in, phi_in = homogeneous_layer_modes(eps_in, period, k0, kx, N) beta_out, phi_out = homogeneous_layer_modes(eps_out, period, k0, kx, N) beta = np.zeros((Nlayers, Ntot), dtype='complex128') phi = np.zeros((Nlayers, Ntot, Ntot), dtype='complex128') for il in range(Nlayers): beta[il,:], phi[il,:,:] = fmm1d_te_layer_modes(eps[il,:], period, k0, kx, N) S = smatrix_interface(beta_in, phi_in, beta[0,:], phi[0,:,:]) for il in range(Nlayers-1): S = stack_smatrices(S,smatrix_propagation(beta[il,:], thick[il])) S = stack_smatrices(S,smatrix_interface(beta[il,:], phi[il,:,:], beta[il+1,:], phi[il+1,:,:])) S = stack_smatrices(S, smatrix_propagation(beta[-1,:], thick[-1])) S = stack_smatrices(S, smatrix_interface(beta[-1,:], phi[-1,:,:], beta_out, phi_out)) # Now, # [u_tr] [u_in] # [u_ref] = S [0] S11 = S[0:Ntot,0:Ntot] S12 = S[0:Ntot,Ntot:] S21 = S[Ntot:,0:Ntot] S22 = S[Ntot:,Ntot:] i_inc = N u_inc = np.zeros(Ntot, dtype='complex128') u_inc[i_inc] = 1.0 u_tr = S11 @ u_inc u_ref = S21 @ u_inc kz_inc = k0*np.sqrt(eps_in)*np.cos(theta) eff_tr = abs(u_tr)**2 * beta_out.real / kz_inc eff_ref = abs(u_ref)**2 * beta_in.real / kz_inc return eff_ref, eff_tr, u_ref, u_tr # Testing the single layer mode solver N = 25 k0 = 2*np.pi/1.064 kx = 0.0 period = 3.0 width = 1.5 x = np.linspace(-period/2, period/2, 1001) x = x[0:-1] # leave the last one out because it's on the periodic boundary eps = np.where(abs(x) <= width/2, 4.0, 1.0) beta, phie = fmm1d_te_layer_modes(eps, period, k0, kx, N) i = np.argsort(beta.real) beta = beta[i] phie = phie[:,i] print('Three largest betas') print(beta[-3:]) u = fourier_to_real(phie, x, period) plt.rcParams['font.size'] = 12 plt.figure() plt.plot(x, u[:,-9:].real + u[:,-9:].imag) plt.legend([ f'beta = {beta[j]}' for j in range(-9,0)]) plt.vlines([-width/2, width/2], -2, 2, 'k') plt.show() # Testing the diffraction efficiency calculation wl0 = 1.064 theta = 0.0 period = 3.0 N = 20 x = np.linspace(0, period, 1001) x = x[0:-1] eps_in = 1.0 eps_low = 1.0 eps_out = 4.0 eps_high = 4.0 #eps_high = 1.0 Nlayers = 3 eps = np.zeros((Nlayers,len(x))) for il in range(Nlayers): eps[il,:] = np.where(x <= (il+1)*period/4, eps_high, eps_low) thick = 0.25*np.ones(Nlayers) k0 = 2*np.pi/wl0 kx = k0*np.sqrt(eps_in)*np.sin(theta) Nx = len(x) beta_in, phi_in = homogeneous_layer_modes(eps_in, period, k0, kx, N) beta_out, phi_out = homogeneous_layer_modes(eps_out, period, k0, kx, N) m = np.arange(2*N+1)-N kxw = kx + 2*np.pi/period*m effr, efft, ur, ut = fmm1d_te(wl0, theta, period, eps_in, eps_out, eps, thick, N) print(f'Sum of efficiencies: tr {efft.sum()}, ref {effr.sum()}, tot {efft.sum()+effr.sum()}') print(ut[N-1:N+2]) print(ur[N-1:N+2]) # The s-matrix implementation should give the same results as the transfer # matrix one. effr, efft, ur, ut = fmm1d_te_smatrix(wl0, theta, period, eps_in, eps_out, eps, thick, N) print(f'Sum of efficiencies: tr {efft.sum()}, ref {effr.sum()}, tot {efft.sum()+effr.sum()}') print(ut[N-1:N+2]) print(ur[N-1:N+2]) plt.figure() plt.plot(kxw, efft, 'o', kxw, effr, 'o') plt.legend(['T side', 'R side']) plt.vlines([k0*np.sqrt(eps_in), k0*np.sqrt(eps_out), -k0*np.sqrt(eps_in), -k0*np.sqrt(eps_out)], 0, 1, 'k') plt.xlabel('m') plt.ylabel('Efficiency') plt.plot() plt.show()