From fdc7f32c10427a43e2cbbe6a7807455f4ca6c007 Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Sun, 9 Sep 2018 23:26:59 -0500 Subject: [PATCH] Added a coupled wave library containing code to implement Bryn's homogeneous layered model --- python/coupledwave.py | 324 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ stim/image/image.h | 77 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------- 2 files changed, 384 insertions(+), 17 deletions(-) create mode 100644 python/coupledwave.py diff --git a/python/coupledwave.py b/python/coupledwave.py new file mode 100644 index 0000000..3002667 --- /dev/null +++ b/python/coupledwave.py @@ -0,0 +1,324 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Sep 6 22:14:36 2018 + +@author: david +""" + +import numpy as np +import matplotlib.pyplot as plt + +#evaluate a plane wave for all points in R = [Rx, Ry, Rz] +# E0 [Ex, Ey, Ez] is the electric field vector +# k is the k vector +# n is the refractive index of the material +# R is a list of coordinates in x, y, and z +def planewave2field(E0, k, n, R): + knr = n * (k[0] * R[0] + k[1] * R[1] + k[2] + R[2]) + exp_2pknr = np.exp(1j * 2 * np.pi * knr) + Ex = E0[0] * exp_2pknr + Ey = E0[1] * exp_2pknr + Ez = E0[2] * exp_2pknr + + return [Ex, Ey, Ez] + +def orthogonalize(E, k): + s = np.cross(E, k) + return np.cross(k, s) + +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 idx(self, l, c, d): + i = l * 6 + d * 3 + c - 3 + print(i) + return i + + #create a matrix for a single plane wave specified by k and E + def solve(self, k, E): + + #calculate the wavenumber + kn = np.linalg.norm(k) + + #s is the plane wave direction scaled by the refractive index + s = k / kn * 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.idx(l, 0, 1)] = s[0] + M[ei, self.idx(l, 1, 1)] = s[1] + M[ei, self.idx(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.idx(l, 0, 0)] = s[0] + M[ei, self.idx(l, 1, 0)] = s[1] + M[ei, self.idx(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 * kn * sz0 * (self.z[l] - self.z[l-1])) + if l < L - 1: + B = np.exp(-1j * kn * sz1 * (self.z[l] - self.z[l+1])) + + + #if this is the second layer, use the simplified equations that account for the incident field + if l == 1: + M[ei, self.idx(0, 0, 1)] = 1 + M[ei, self.idx(1, 0, 0)] = -1 + if L > 2: + M[ei, self.idx(1, 0, 1)] = -B + + b[ei] = -A * E[0] + ei = ei + 1 + + M[ei, self.idx(0, 1, 1)] = 1 + M[ei, self.idx(1, 1, 0)] = -1 + if L > 2: + M[ei, self.idx(1, 1, 1)] = -B + b[ei] = -A * E[1] + ei = ei + 1 + + M[ei, self.idx(0, 2, 1)] = s[1] + M[ei, self.idx(0, 1, 1)] = sz0 + M[ei, self.idx(1, 2, 0)] = -s[1] + M[ei, self.idx(1, 1, 0)] = sz1 + if L > 2: + M[ei, self.idx(1, 2, 1)] = -B*s[1] + M[ei, self.idx(1, 1, 1)] = -B*sz1 + b[ei] = A * sz0 * E[1] - A * s[1]*E[2] + ei = ei + 1 + + M[ei, self.idx(0, 0, 1)] = -sz0 + M[ei, self.idx(0, 2, 1)] = -s[0] + M[ei, self.idx(1, 0, 0)] = -sz1 + M[ei, self.idx(1, 2, 0)] = s[0] + if L > 2: + M[ei, self.idx(1, 0, 1)] = B*sz1 + M[ei, self.idx(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.idx(l-1, 0, 0)] = A + M[ei, self.idx(l-1, 0, 1)] = 1 + M[ei, self.idx(l, 0, 0)] = -1 + ei = ei + 1 + + M[ei, self.idx(l-1, 1, 0)] = A + M[ei, self.idx(l-1, 1, 1)] = 1 + M[ei, self.idx(l, 1, 0)] = -1 + ei = ei + 1 + + M[ei, self.idx(l-1, 2, 0)] = A*s[1] + M[ei, self.idx(l-1, 1, 0)] = -A*sz0 + M[ei, self.idx(l-1, 2, 1)] = s[1] + M[ei, self.idx(l-1, 1, 1)] = sz0 + M[ei, self.idx(l, 2, 0)] = -s[1] + M[ei, self.idx(l, 1, 0)] = sz1 + ei = ei + 1 + + M[ei, self.idx(l-1, 0, 0)] = A*sz0 + M[ei, self.idx(l-1, 2, 0)] = -A*s[0] + M[ei, self.idx(l-1, 0, 1)] = -sz0 + M[ei, self.idx(l-1, 2, 1)] = -s[0] + M[ei, self.idx(l, 0, 0)] = -sz1 + M[ei, self.idx(l, 2, 0)] = s[0] + ei = ei + 1 + #otherwise use the full set of boundary conditions + else: + M[ei, self.idx(l-1, 0, 0)] = A + M[ei, self.idx(l-1, 0, 1)] = 1 + M[ei, self.idx(l, 0, 0)] = -1 + M[ei, self.idx(l, 0, 1)] = -B + ei = ei + 1 + + M[ei, self.idx(l-1, 1, 0)] = A + M[ei, self.idx(l-1, 1, 1)] = 1 + M[ei, self.idx(l, 1, 0)] = -1 + M[ei, self.idx(l, 1, 1)] = -B + ei = ei + 1 + + M[ei, self.idx(l-1, 2, 0)] = A*s[1] + M[ei, self.idx(l-1, 1, 0)] = -A*sz0 + M[ei, self.idx(l-1, 2, 1)] = s[1] + M[ei, self.idx(l-1, 1, 1)] = sz0 + M[ei, self.idx(l, 2, 0)] = -s[1] + M[ei, self.idx(l, 1, 0)] = sz1 + M[ei, self.idx(l, 2, 1)] = -B*s[1] + M[ei, self.idx(l, 1, 1)] = -B*sz1 + ei = ei + 1 + + M[ei, self.idx(l-1, 0, 0)] = A*sz0 + M[ei, self.idx(l-1, 2, 0)] = -A*s[0] + M[ei, self.idx(l-1, 0, 1)] = -sz0 + M[ei, self.idx(l-1, 2, 1)] = -s[0] + M[ei, self.idx(l, 0, 0)] = -sz1 + M[ei, self.idx(l, 2, 0)] = s[0] + M[ei, self.idx(l, 0, 1)] = B*sz1 + M[ei, self.idx(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.Ptrans = [] + self.Prefl = [] + for l in range(L): + if l == 0: + self.Ptrans.append((E[0], E[1], E[2])) + else: + px = P[self.idx(l, 0, 0)] + py = P[self.idx(l, 1, 0)] + pz = P[self.idx(l, 2, 0)] + self.Ptrans.append((px, py, pz)) + + if l == L-1: + self.Prefl.append((0, 0, 0)) + else: + px = P[self.idx(l, 0, 1)] + py = P[self.idx(l, 1, 1)] + pz = P[self.idx(l, 2, 1)] + self.Prefl.append((px, py, pz)) + + #this converts the coefficients into a NumPy array for convenience later on + self.Ptrans = np.array(self.Ptrans) + self.Prefl = np.array(self.Prefl) + + #store values required for evaluation + #store k + self.k = kn + + #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])) + 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.Ptrans[LI] * Ph_t[:, :, None] + Er = self.Prefl[LI] * Ph_r[:, :, None] + + #add everything together coherently + E = (Et + Er) * Ph_xy[:, :, None] + + #return the electric field + return E + + +#This sample code produces a field similar to Figure 2 in Davis et al., 2010 +#set the material properties +depths = [-50, -15, 15, 40] +n = [1.0, 1.4+1j*0.05, 1.4, 1.0] + +#create a layered sample +layers = layersample(n, depths) + +#set the input light parameters +d = np.array([0.2, 0, 1]) +d = d / np.linalg.norm(d) +l = 5.0 +k = 2 * np.pi / l * d +E0 = [0, 1, 0] +E0 = orthogonalize(E0, d) + +#solve for the substrate field +layers.solve(k, E0) + +#set the simulation domain +N = 512 +M = 1024 +D = [-50, 80] + +x = np.linspace(D[0], D[1], 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 = Er[..., 0] ** 2 + Er[..., 1] **2 + Er[..., 2] ** 2 +plt.figure() +plt.set_cmap("afmhot") +plt.imshow(Er[..., 1]) diff --git a/stim/image/image.h b/stim/image/image.h index 9c9e42c..1898d01 100644 --- a/stim/image/image.h +++ b/stim/image/image.h @@ -33,6 +33,7 @@ class image{ T* img; //pointer to the image data (interleaved RGB for color) size_t R[3]; + bool interleaved = true; //by default the data is interleaved inline size_t X() const { return R[1]; } inline size_t Y() const { return R[2]; } @@ -73,13 +74,24 @@ class image{ #ifdef USING_OPENCV int cv_type(){ - if(typeid(T) == typeid(unsigned char)) return CV_MAKETYPE(CV_8U, (int)C()); - if(typeid(T) == typeid(char)) return CV_MAKETYPE(CV_8S, (int)C()); - if(typeid(T) == typeid(unsigned short)) return CV_MAKETYPE(CV_16U, (int)C()); - if(typeid(T) == typeid(short)) return CV_MAKETYPE(CV_16S, (int)C()); - if(typeid(T) == typeid(int)) return CV_MAKETYPE(CV_32S, (int)C()); - if(typeid(T) == typeid(float)) return CV_MAKETYPE(CV_32F, (int)C()); - if(typeid(T) == typeid(double)) return CV_MAKETYPE(CV_64F, (int)C()); + if (C() == 1 || C() == 3) { //if the image has 1 or 3 channels, this will work + if (typeid(T) == typeid(unsigned char)) return CV_MAKETYPE(CV_8U, (int)C()); + if (typeid(T) == typeid(char)) return CV_MAKETYPE(CV_8S, (int)C()); + if (typeid(T) == typeid(unsigned short)) return CV_MAKETYPE(CV_16U, (int)C()); + if (typeid(T) == typeid(short)) return CV_MAKETYPE(CV_16S, (int)C()); + if (typeid(T) == typeid(int)) return CV_MAKETYPE(CV_32S, (int)C()); + if (typeid(T) == typeid(float)) return CV_MAKETYPE(CV_32F, (int)C()); + if (typeid(T) == typeid(double)) return CV_MAKETYPE(CV_64F, (int)C()); + } + else if (C() == 4) { //OpenCV only supports saving BGR images - 4-channel isn't supported + if (typeid(T) == typeid(unsigned char)) return CV_MAKETYPE(CV_8U, 3); + if (typeid(T) == typeid(char)) return CV_MAKETYPE(CV_8S, 3); + if (typeid(T) == typeid(unsigned short)) return CV_MAKETYPE(CV_16U, 3); + if (typeid(T) == typeid(short)) return CV_MAKETYPE(CV_16S, 3); + if (typeid(T) == typeid(int)) return CV_MAKETYPE(CV_32S, 3); + if (typeid(T) == typeid(float)) return CV_MAKETYPE(CV_32F, 3); + if (typeid(T) == typeid(double)) return CV_MAKETYPE(CV_64F, 3); + } std::cout<<"ERROR in stim::image::cv_type - no valid data type found"<