# -*- coding: utf-8 -*- """ Created on Fri Jun 22 16:22:17 2018 @author: david """ import numpy as np import scipy as sp import scipy.ndimage import matplotlib import math import matplotlib.pyplot as plt #adjust the components of E0 so that they are orthogonal to d #if d is 2D, the z component is calculated such that ||[dx, dy, dz]||=1 def orthogonalize(E, d): if d.shape[-1] == 2: dz = 1 - d[..., 0]**2 - d[..., 1]**2 d = np.append(d, dz) s = np.cross(E, d) dE = np.cross(d, s) return dE #return dE/np.linalg.norm(dE) * np.linalg.norm(E) def intensity(E): Econj = np.conj(E) I = np.sum(E*Econj, axis=-1) return np.real(I) class planewave: #implement all features of a plane wave # k, E, frequency (or wavelength in a vacuum)--l # try to enforce that E and k have to be orthogonal #initialization function that sets these parameters when a plane wave is created def __init__ (self, k, E): self.k = np.array(k) self.E = E if ( np.linalg.norm(k) > 1e-15 and np.linalg.norm(E) >1e-15): s = np.cross(k, E) #compute an orthogonal side vector s = s / np.linalg.norm(s) #normalize it Edir = np.cross(s, k) #compute new E vector which is orthogonal self.k = np.array(k) self.E = Edir / np.linalg.norm(Edir) * np.linalg.norm(E) def __str__(self): return str(self.k) + "\n" + str(self.E) #for verify field vectors use print command def kmag(self): return math.sqrt(self.k[0] ** 2 + self.k[1] ** 2 + self.k[2] ** 2) #function that renders the plane wave given a set of coordinates def evaluate(self, X, Y, Z): k_dot_r = self.k[0] * X + self.k[1] * Y + self.k[2] * Z #phase term k*r ex = np.exp(1j * k_dot_r) #E field equation = E0 * exp (i * (k * r)) here we simply set amplitude as 1 Ef = ex[..., None] * self.E return Ef class objective: #initialize an objective lens # sin_alpha is the sine of the angle subtended by the objective (numerical aperture in a refractive index of 1.0) # n is the real refractive index of the imaging media (default air) # N is the resolution of the 2D objective map (field at the input and output aperture) def __init__(self, sin_alpha, n, N): #set member variables self.sin_alpha = sin_alpha self.n = n #calculate the transfer function from the back aperture to the front aperture self.N = N self._CI(N) self._CR(N) #calculates the basis vectors describing the mapping between the back and front aperture def _calc_bases(self, N): #calculate the (sx, sy) coordinates (angular spectrum) corresponding to each position in the entrance/exit pupil s = np.linspace(-self.sin_alpha, self.sin_alpha, N) [SX, SY] = np.meshgrid(s, s) S_SQ = SX**2 + SY**2 #calculate sx^2 + sy^2 BAND_M = S_SQ <= self.sin_alpha**2 #find a mask for pixels within the bandpass of the objective self.BAND_M = BAND_M #pull out pixels that pass through the objective bandpass SX = SX[BAND_M] SY = SY[BAND_M] S_SQ = S_SQ[BAND_M] SZ = np.sqrt(1 - S_SQ) #calculate sz = sqrt(1 - sx^2 - sy^2) in normalized (-1, 1) coordinates #calculate the apodization functions for each polarization self.FS = 1.0 / np.sqrt(SZ) self.FP = self.FS #calculate the S and P basis vectors THETA = np.arctan2(SY, SX) SX_SXSY = np.cos(THETA) SY_SXSY = np.sin(THETA) v_size = (len(SX), 3) VS_VSP = np.zeros(v_size) VP = np.zeros(v_size) VPP = np.zeros(v_size) VS_VSP[:, 0] = -SY_SXSY VS_VSP[:, 1] = SX_SXSY VP[:, 0] = -SX_SXSY VP[:, 1] = -SY_SXSY SXSY = np.sqrt(S_SQ) VPPZ = np.zeros(SXSY.shape) DC_M = SXSY > 0 VPPZ[DC_M] = S_SQ[DC_M] / SXSY[DC_M] VPP[:, 0] = -SX_SXSY * SZ VPP[:, 1] = -SY_SXSY * SZ VPP[:, 2] = VPPZ return [VS_VSP, VP, VPP] #calculate the transfer function that maps the back aperture field to the front aperture field (Eq. 19 in Davis et al.) def _CI(self, N): [VS_VSP, VP, VPP] = self._calc_bases(N) VSP_VST = np.matmul(VS_VSP[:, :, None], VS_VSP[:, None, :]) VPP_VPT = np.matmul(VPP[:, :, None], VP[:, None, :]) self.CI = self.FS[:, None, None] * VSP_VST + self.FP[:, None, None] * VPP_VPT def _CR(self, N): [VS_VSP, VP, VPP] = self._calc_bases(N) VSP_VST = np.matmul(VS_VSP[:, :, None], VS_VSP[:, None, :]) VP_VPPT = np.matmul(VP[:, :, None], VPP[:, None, :]) self.CR = 1.0 / self.FS[:, None, None] * VSP_VST + 1.0 / self.FP[:, None, None] * VP_VPPT #condenses a field defined at the back aperture of the objective # E is a 3xNxN image of the field at the back aperture def back2front(self, E): Ein = E[self.BAND_M, :, None] Efront = np.zeros((self.N, self.N, 3), dtype=np.complex128) Efront[self.BAND_M, :] = np.squeeze(np.matmul(self.CI, Ein)) return Efront #propagates a field from the front aperture of the objective to the back aperture # E is a 3xNxN image of the field at the front aperture def front2back(self, E): Ein = E[self.BAND_M, :, None] Eback = np.zeros((self.N, self.N, 3), dtype=np.complex128) Eback[self.BAND_M, :] = np.squeeze(np.matmul(self.CR, Ein)) return Eback #calculates the image at the focus given the image at the back aperture # E is the field at the back aperture of the objective # k0 is the wavenumber in a vacuum (2pi/lambda_0) def focus(self, E, k0, RX, RY, RZ): #calculate the angular spectrum at the front aperture Ef = self.back2front(E) #generate an array describing the valid components of S s = np.linspace(-self.sin_alpha, self.sin_alpha, self.N) [SX, SY] = np.meshgrid(s, s) S = [SX[self.BAND_M], SY[self.BAND_M]] SXSY_SQ = S[0]**2 + S[1]**2 S.append(np.sqrt(1 - SXSY_SQ)) S = np.array(S) #calculate the k vector (accounting for both k0 and the refractive index of the immersion medium) K = k0 * self.n * np.transpose(S)[:, None, :] #calculate the dot product of K and R R = np.moveaxis(np.array([RX, RY, RZ]), 0, -1)[..., None, :, None] K_DOT_R = np.squeeze(np.matmul(K, R)) EX = np.exp(1j * K_DOT_R)[..., None] Ewaves = EX * Ef[self.BAND_M] Efocus = np.sum(Ewaves, -2) return Efocus class gaussian_beam: #generate a beam # (Ex, Ey) specifies the polarization # p specifies the focal point # l is the free-space wavelength # sigma_p is the standard deviation of the Gaussian beam at the focal point # amp is the amplitude of the DC component of the Gaussian beam def __init__(self, x_polar, y_polar, lambda_0=1, amp=1, sigma_p=None, sigma_k=None): self.E = [x_polar, y_polar] self.l = lambda_0 #calculate the standard deviation of the Gaussian function in k-space if(sigma_p is None and sigma_k is None): self.sigma = 2 * math.pi / lambda_0 elif(sigma_p is not None): self.sigma = 1.0 / sigma_p else: self.sigma = sigma_k #self.s = self.sigma_p2k(sigma_p) self.A = amp #calculates the maximum k-space frequency supported by this beam def kmax(self): return 2 * np.pi / self.l #calculate the standard deviation of the Gaussian function in k space given the std at the focal point # def sigma_p2k(self, sigma): # this code is based on: http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf # f(r) = 1 / (2 pi sigma^2) e^(-r^2 / (2 sigma^2)) # F(kx, ky) = f( |k_par| ) = e^( -2 pi^2 |k_par|^2 sigma^2 ) # Basically, the Fourier transform of a normally distributed Gaussian is # a non-scaled Gaussian with standard deviation 1/(2 pi sigma) #return 1.0/(2 * math.pi * sigma) #return 1.0 / sigma def kspace(self, Kx, Ky): #sigma = self.s * (np.pi / self.l) #G = 1 / (2 * np.pi * sigma **2) * np.exp(- (0.5/(sigma**2)) * (Kx ** 2 + Ky ** 2)) G = self.A * np.exp(- (Kx ** 2 + Ky ** 2) / (2 * self.sigma ** 2)) #calculate k-vectors outside of the band-limit B = ((Kx ** 2 + Ky **2) > self.kmax() ** 2) #calculate Kz (required for calculating the Ez component of E) #Kz = np.lib.scimath.sqrt(self.kmax() ** 2 - Kx**2 - Ky**2) Ex = self.E[0] * G Ey = self.E[1] * G #Ez = -Kx/Kz * Ex Ez = np.zeros(Ex.shape) #set all values outside of the band limit to zero Ex[B] = 0 Ey[B] = 0 Ez[B] = 0 return [Ex, Ey, Ez] # The PSF at the focal point is the convolution of a Bessel function and a Gaussian # This function calculates the parameters of the Gaussian function in terms of a*e^{-x^2 / (2*c^2)} # a (return) is the scale factor for the Gaussian # c (return) is the standard deviation def psf_gaussian(self): # this code is based on: http://mathworld.wolfram.com/FourierTransformGaussian.html #a = math.sqrt( 2 * math.pi * (self.s ** 2) ) #c = 1.0 / (2 * math.pi * self.s) a = math.sqrt( (self.sigma ** 2) ) c = 1.0 / (self.sigma) return [a, c] # The PSF at the focal point is the convolution of a Bessel function and a Gaussian # This function calculates the parameters of the bessel function in terms of a J_1(pi * a * x) / x def psf_bessel(self): return self.kmax() / math.pi # calculate the beam point spread function at the focal plane # D is the image extent for the square: (-D, D, -D, D) # N is the number of spatial samples along one dimension of the image (NxN) # order is the mathematical order of the calculation (in terms of the resolution of the FFT) def psf(self, D, N): [a, c] = self.psf_gaussian() alpha = self.psf_bessel() #calculate the sample positions along the x-axis x = np.linspace(-D, D, N) dx = x[1] - x[0] [X, Y] = np.meshgrid(x, x) R = np.sqrt(X**2 + Y**2) bessel = alpha * sp.special.j1(math.pi * alpha * R) / R return a * sp.ndimage.filters.gaussian_filter(bessel, c / dx) #return the k-space transform of the beam at p as an NxN array def kspace_image(self, N): #calculate the band limit [-2*pi/lambda, 2*pi/lambda] kx = np.linspace(-self.kmax(), self.kmax(), N) ky = kx #generate a meshgrid for the field [Kx, Ky] = np.meshgrid(kx, ky) #generate a Gaussian with a standard deviation of sigma * pi/lambda return self.kspace(Kx, Ky) # calculate the practical band limit of the beam such that # all frequencies higher than the returned values are less than epsilon def practical_band_limit(self, epsilon = 0.001): return min(self.sigma * math.sqrt(math.log(1.0/epsilon)), self.kmax()) #return a plane wave decomposition of the beam using N plane waves def decompose(self, N): PW = [] #create an empty list of plane waves # sample a cylinder and project those samples onto a unit sphere theta = np.random.uniform(0, 2*np.pi, N) # get the practical band limit to minimize sampling band_limit = self.practical_band_limit() / self.kmax() phi = np.arccos(np.random.uniform(1 - band_limit, 1, N)) kx = -self.kmax() * np.sin(phi) * np.cos(theta) ky = -self.kmax() * np.sin(phi) * np.sin(theta) mc_weight = 2 * math.pi / N E0 = self.kspace(kx, ky) kz = np.sqrt(self.kmax() ** 2 - kx**2 - ky**2) for i in range(0, N): k = [kx[i], ky[i], kz[i]] E = [E0[0][i] * mc_weight, E0[1][i] * mc_weight, E0[2][i] * mc_weight] w = planewave(k, E) PW.append(w) return PW class layersample: #generate a layered sample based on a list of layer refractive indices # layer_ri is a list of complex refractive indices for each layer # z is the list of boundary positions along z def __init__(self, layer_ri, z): self.n = np.array(layer_ri).astype(np.complex128) self.z = np.array(z) #calculate the index of the field component associated with a layer # l is the layer index [0, L) # c is the component (x=0, y=1, z=2) # d is the direction (0 = transmission, 1 = reflection) def i(self, l, c, d): i = l * 6 + d * 3 + c - 3 return i #create a matrix for a single plane wave specified by k and E # d = [dx, dy] are the x and y coordinates of the normalized direction of propagation # k0 is the free space wave number (2 pi / lambda0) # E is the electric field vector def solve1(self, d, k0, E): #s is the plane wave direction scaled by the refractive index s = np.array(d) * self.n[0] #allocate space for the matrix L = len(self.n) M = np.zeros((6*(L-1), 6*(L-1)), dtype=np.complex128) #allocate space for the RHS vector b = np.zeros(6*(L-1), dtype=np.complex128) #initialize a counter for the equation number ei = 0 #calculate the sz component for each layer self.sz = np.zeros(L, dtype=np.complex128) for l in range(L): self.sz[l] = np.sqrt(self.n[l]**2 - s[0]**2 - s[1]**2) #set constraints based on Gauss' law for l in range(0, L): #sz = np.sqrt(self.n[l]**2 - s[0]**2 - s[1]**2) #set the upward components for each layer # note that layer L-1 does not have a downward component # David I, Equation 7 if l != L-1: M[ei, self.i(l, 0, 1)] = s[0] M[ei, self.i(l, 1, 1)] = s[1] M[ei, self.i(l, 2, 1)] = -self.sz[l] ei = ei+1 #set the downward components for each layer # note that layer 0 does not have a downward component # Davis I, Equation 6 if l != 0: M[ei, self.i(l, 0, 0)] = s[0] M[ei, self.i(l, 1, 0)] = s[1] M[ei, self.i(l, 2, 0)] = self.sz[l] ei = ei+1 #enforce a continuous field across boundaries for l in range(1, L): sz0 = self.sz[l-1] sz1 = self.sz[l] A = np.exp(1j * k0 * sz0 * (self.z[l] - self.z[l-1])) if l < L - 1: dl = self.z[l] - self.z[l+1] arg = -1j * k0 * sz1 * dl B = np.exp(arg) #if this is the second layer, use the simplified equations that account for the incident field if l == 1: M[ei, self.i(0, 0, 1)] = 1 M[ei, self.i(1, 0, 0)] = -1 if L > 2: #print(-B, M[ei, self.i(1, 0, 1)]) M[ei, self.i(1, 0, 1)] = -B b[ei] = -A * E[0] ei = ei + 1 M[ei, self.i(0, 1, 1)] = 1 M[ei, self.i(1, 1, 0)] = -1 if L > 2: M[ei, self.i(1, 1, 1)] = -B b[ei] = -A * E[1] ei = ei + 1 M[ei, self.i(0, 2, 1)] = s[1] M[ei, self.i(0, 1, 1)] = sz0 M[ei, self.i(1, 2, 0)] = -s[1] M[ei, self.i(1, 1, 0)] = sz1 if L > 2: M[ei, self.i(1, 2, 1)] = -B*s[1] M[ei, self.i(1, 1, 1)] = -B*sz1 b[ei] = A * sz0 * E[1] - A * s[1]*E[2] ei = ei + 1 M[ei, self.i(0, 0, 1)] = -sz0 M[ei, self.i(0, 2, 1)] = -s[0] M[ei, self.i(1, 0, 0)] = -sz1 M[ei, self.i(1, 2, 0)] = s[0] if L > 2: M[ei, self.i(1, 0, 1)] = B*sz1 M[ei, self.i(1, 2, 1)] = B*s[0] b[ei] = A * s[0] * E[2] - A * sz0 * E[0] ei = ei + 1 #if this is the last layer, use the simplified equations that exclude reflections from the last layer elif l == L-1: M[ei, self.i(l-1, 0, 0)] = A M[ei, self.i(l-1, 0, 1)] = 1 M[ei, self.i(l, 0, 0)] = -1 ei = ei + 1 M[ei, self.i(l-1, 1, 0)] = A M[ei, self.i(l-1, 1, 1)] = 1 M[ei, self.i(l, 1, 0)] = -1 ei = ei + 1 M[ei, self.i(l-1, 2, 0)] = A*s[1] M[ei, self.i(l-1, 1, 0)] = -A*sz0 M[ei, self.i(l-1, 2, 1)] = s[1] M[ei, self.i(l-1, 1, 1)] = sz0 M[ei, self.i(l, 2, 0)] = -s[1] M[ei, self.i(l, 1, 0)] = sz1 ei = ei + 1 M[ei, self.i(l-1, 0, 0)] = A*sz0 M[ei, self.i(l-1, 2, 0)] = -A*s[0] M[ei, self.i(l-1, 0, 1)] = -sz0 M[ei, self.i(l-1, 2, 1)] = -s[0] M[ei, self.i(l, 0, 0)] = -sz1 M[ei, self.i(l, 2, 0)] = s[0] ei = ei + 1 #otherwise use the full set of boundary conditions else: M[ei, self.i(l-1, 0, 0)] = A M[ei, self.i(l-1, 0, 1)] = 1 M[ei, self.i(l, 0, 0)] = -1 M[ei, self.i(l, 0, 1)] = -B ei = ei + 1 M[ei, self.i(l-1, 1, 0)] = A M[ei, self.i(l-1, 1, 1)] = 1 M[ei, self.i(l, 1, 0)] = -1 M[ei, self.i(l, 1, 1)] = -B ei = ei + 1 M[ei, self.i(l-1, 2, 0)] = A*s[1] M[ei, self.i(l-1, 1, 0)] = -A*sz0 M[ei, self.i(l-1, 2, 1)] = s[1] M[ei, self.i(l-1, 1, 1)] = sz0 M[ei, self.i(l, 2, 0)] = -s[1] M[ei, self.i(l, 1, 0)] = sz1 M[ei, self.i(l, 2, 1)] = -B*s[1] M[ei, self.i(l, 1, 1)] = -B*sz1 ei = ei + 1 M[ei, self.i(l-1, 0, 0)] = A*sz0 M[ei, self.i(l-1, 2, 0)] = -A*s[0] M[ei, self.i(l-1, 0, 1)] = -sz0 M[ei, self.i(l-1, 2, 1)] = -s[0] M[ei, self.i(l, 0, 0)] = -sz1 M[ei, self.i(l, 2, 0)] = s[0] M[ei, self.i(l, 0, 1)] = B*sz1 M[ei, self.i(l, 2, 1)] = B*s[0] ei = ei + 1 #store the matrix and RHS vector (for debugging) self.M = M self.b = b #evaluate the linear system P = np.linalg.solve(M, b) #save the results (also for debugging) self.P = P #store the coefficients for each layer self.Pt = np.zeros((3, L), np.complex128) self.Pr = np.zeros((3, L), np.complex128) for l in range(L): if l == 0: self.Pt[:, 0] = [E[0], E[1], E[2]] else: px = P[self.i(l, 0, 0)] py = P[self.i(l, 1, 0)] pz = P[self.i(l, 2, 0)] self.Pt[:, l] = [px, py, pz] if l == L-1: self.Pr[:, L-1] = [0, 0, 0] else: px = P[self.i(l, 0, 1)] py = P[self.i(l, 1, 1)] pz = P[self.i(l, 2, 1)] self.Pr[:, l] = [px, py, pz] #store values required for evaluation #store k self.k = k0 #store sx and sy self.s = np.array([s[0], s[1]]) self.solved = True #evaluate a solved homogeneous substrate def evaluate(self, X, Y, Z): if not self.solved: print("ERROR: the layered substrate hasn't been solved") return #this code is a bit cumbersome and could probably be optimized # Basically, it vectorizes everything by creating an image # such that the value at each pixel is the corresponding layer # that the pixel resides in. That index is then used to calculate # the field within the layer #allocate space for layer indices LI = np.zeros(Z.shape, dtype=np.int) #find the layer index for each sample point L = len(self.z) LI[Z < self.z[0]] = 0 for l in range(L-1): idx = np.logical_and(Z > self.z[l], Z <= self.z[l+1]) LI[idx] = l LI[Z > self.z[-1]] = L - 1 #calculate the appropriate phase shift for the wave transmitted through the layer Ph_t = np.exp(1j * self.k * self.sz[LI] * (Z - self.z[LI])) #calculate the appropriate phase shift for the wave reflected off of the layer boundary LIp = LI + 1 LIp[LIp >= L] = 0 Ph_r = np.exp(-1j * self.k * self.sz[LI] * (Z - self.z[LIp])) # print(Ph_r) Ph_r[LI >= L-1] = 0 #calculate the phase shift based on the X and Y positions Ph_xy = np.exp(1j * self.k * (self.s[0] * X + self.s[1] * Y)) #apply the phase shifts Et = self.Pt[:, LI] * Ph_t[:, :] Er = self.Pr[:, LI] * Ph_r[:, :] #add everything together coherently E = (Et + Er) * Ph_xy[:, :] #return the electric field return np.moveaxis(E, 0, -1) def example_psf(): #specify the angle of accepted waves (NA if n = 1) sin_angle = 0.8 #refractive index of the imaging medium n = 1.0 #specify the size of the spectral domain (resolution of the front/back aperture) N = 100 #create a new objective, specifying the NA (angle * refractive index) and resolution of the transform O = objective(sin_angle, n, N) #specify the size of the spatial domain D = 10 #specify the resolution of the spatial domain M = 100 #free-space wavelength l0 = 1 #calculate the free-space wavenumber k0 = 2 * np.pi * l0 #specify the field incident at the aperture pap = np.array([1, 0, 0]) #polarization sap = np.array([0, 0, 1]) #incident angle of the plane wave (will be normalized) sap = sap / np.linalg.norm(sap) #calculate the normalized direction vector for the plane wave k0ap = sap * k0 #calculate the k-vector of the plane wave incident on the aperture #create a plane wave pw = planewave(k0ap, pap) #generate the domain for the simulation (where the focus will be simulated) x = np.linspace(-D, D, M) z = np.linspace(-D, D, M) [RX, RZ] = np.meshgrid(x, z) RY = np.zeros(RX.shape) #create an image of the back aperture (just an incident plane wave) xap = np.linspace(-D, D, N) [XAP, YAP] = np.meshgrid(xap, xap) ZAP = np.zeros(XAP.shape) #evaluate the plane wave at the back aperture to produce a 2D image of the field EAP = pw.evaluate(XAP, YAP, ZAP) #calculate the field at the focus (coordinates are specified in RX, RY, and RZ) Efocus = O.focus(EAP, k0, RX, RY, RZ) #display an image of the calculated field plt.figure() plt.imshow(intensity(Efocus)) plt.colorbar() return Efocus #This sample code produces a field similar to Figure 2 in Davis et al., 2010 def example_layer(): #set the material properties depths = [-50, -15, 15, 40] #specify the position (first boundary) of each layer (first boundary is where the field is defined) n = [1.0, 1.4+1j*0.05, 1.4, 1.0] #specify the refractive indices of each layer #create a layered sample layers = layersample(n, depths) #set the input light parameters d = np.array([0.5, 0]) #direction of propagation of the plane wave #d = d / np.linalg.norm(d) #normalize this direction vector l0 = 5 #specify the wavelength in free-space k0 = 2 * np.pi / l0 #calculate the wavenumber in free-space E0 = [1, 1, 0] #specify the E vector E0 = orthogonalize(E0, d) #make sure that both vectors are orthogonal #solve for the substrate field layers.solve1(d, k0, E0) #set the simulation domain N = 512 M = 1024 D = [-40, 80, 0, 60] x = np.linspace(D[2], D[3], N) z = np.linspace(D[0], D[1], M) [X, Z] = np.meshgrid(x, z) Y = np.zeros(X.shape) E = layers.evaluate(X, Y, Z) Er = np.real(E) I = intensity(E) plt.figure() plt.set_cmap("afmhot") matplotlib.rcParams.update({'font.size': 32}) plt.subplot(1, 4, 1) plt.imshow(Er[..., 0], extent=(D[3], D[2], D[1], D[0])) plt.title("Ex") plt.subplot(1, 4, 2) plt.imshow(Er[..., 1], extent=(D[3], D[2], D[1], D[0])) plt.title("Ey") plt.subplot(1, 4, 3) plt.imshow(Er[..., 2], extent=(D[3], D[2], D[1], D[0])) plt.title("Ez") plt.subplot(1, 4, 4) plt.imshow(I, extent=(D[3], D[2], D[1], D[0])) plt.title("I") def example_psflayer(): #specify the angle of accepted waves (NA if n = 1) sin_angle = 0.7 #refractive index of the imaging medium n = 1.0 #specify the size of the spectral domain (resolution of the front/back aperture) N = 31 #create a new objective, specifying the NA (angle * refractive index) and resolution of the transform O = objective(sin_angle, n, N) #specify the size of the spatial domain D = 5 #specify the resolution of the spatial domain M = 200 #free-space wavelength l0 = 1 #calculate the free-space wavenumber k0 = 2 * np.pi /l0 #specify the field incident at the aperture pap = np.array([1, 0, 0]) #polarization sap = np.array([0, 0, 1]) #incident angle of the plane wave (will be normalized) sap = sap / np.linalg.norm(sap) #calculate the normalized direction vector for the plane wave k0ap = sap * k0 #calculate the k-vector of the plane wave incident on the aperture #create a plane wave pw = planewave(k0ap, pap) #create an image of the back aperture (just an incident plane wave) xap = np.linspace(-D, D, N) [XAP, YAP] = np.meshgrid(xap, xap) ZAP = np.zeros(XAP.shape) #evaluate the plane wave at the back aperture to produce a 2D image of the field Eback = pw.evaluate(XAP, YAP, ZAP) Efront = O.back2front(Eback) #set the material properties depths = [0, -1, 0, 2] #specify the position (first boundary) of each layer (first boundary is where the field is defined) n = [1.0, 1.4+1j*0.05, 1.4, 1.0] #specify the refractive indices of each layer #create a layered sample L = layersample(n, depths) #generate the domain for the simulation (where the focus will be simulated) x = np.linspace(-D, D, M) z = np.linspace(-D, D, M) [RX, RZ] = np.meshgrid(x, z) RY = np.zeros(RX.shape) i = 0 ang = np.linspace(-sin_angle, sin_angle, N) for yi in range(N): dy = ang[yi] for xi in range(N): dx = ang[xi] dxdy_sq = dx**2 + dy**2 if dxdy_sq <= sin_angle**2: d = np.array([dx, dy]) e = Efront[xi, yi, :] L.solve1(d, k0, e) if i == 0: E = L.evaluate(RX, RY, RZ) i = i + 1 else: E = E + L.evaluate(RX, RY, RZ) plt.imshow(intensity(E)) def example_spp(): l = 0.647 k0 = 2 * np.pi / l n_oil = 1.51 n_mica = 1.56 n_gold = 0.16049 + 1j * 3.5726 n_water = 1.33 #all units are in um d_oil = 100 d_mica = 80 d_gold = 0.05 #specify the parameters for the layered sample zp = [0, d_oil, d_oil + d_mica, d_oil + d_mica + d_gold] n = [n_oil, n_mica, n_gold, n_water] #generate a layered sample L = layersample(n, zp) N = 10000 min_theta = -90 max_theta = 90 DEG = np.linspace(min_theta, max_theta, N) I = np.zeros(N) for i in range(N): theta = DEG[i] * np.pi / 180 #create a plane wave x = np.sin(theta) z = np.cos(theta) d = [x, 0] #calculate the E vector by finding a value orthogonal to d E0 = [-z, 0, x] #solve for the layered sample L.solve1(d, k0, E0) Er = L.Pr[:, 0] I[i] = intensity(Er) plt.plot(I)