from __future__ import division, print_function, absolute_import from numpy import cos, inf, zeros, array, exp, conj, nan, isnan, pi, sin import numpy as np import scipy as sp import sys EPSILON = sys.float_info.epsilon # typical floating-point calculation error def make_2x2_array(a, b, c, d, dtype=float): """ Makes a 2x2 numpy array of [[a,b],[c,d]] Same as "numpy.array([[a,b],[c,d]], dtype=float)", but ten times faster """ my_array = np.empty((2,2), dtype=dtype) my_array[0,0] = a my_array[0,1] = b my_array[1,0] = c my_array[1,1] = d return my_array def is_forward_angle(n, theta): """ if a wave is traveling at angle theta from normal in a medium with index n, calculate whether or not this is the forward-traveling wave (i.e., the one going from front to back of the stack, like the incoming or outgoing waves, but unlike the reflected wave). For real n & theta, the criterion is simply -pi/2 < theta < pi/2, but for complex n & theta, it's more complicated. See https://arxiv.org/abs/1603.02720 appendix D. If theta is the forward angle, then (pi-theta) is the backward angle and vice-versa. """ assert n.real * n.imag >= 0, ("For materials with gain, it's ambiguous which " "beam is incoming vs outgoing. See " "https://arxiv.org/abs/1603.02720 Appendix C.\n" "n: " + str(n) + " angle: " + str(theta)) ncostheta = n * np.cos(theta) if abs(ncostheta.imag) > 100 * EPSILON: # Either evanescent decay or lossy medium. Either way, the one that # decays is the forward-moving wave answer = (ncostheta.imag > 0) else: # Forward is the one with positive Poynting vector # Poynting vector is Re[n cos(theta)] for s-polarization or # Re[n cos(theta*)] for p-polarization, but it turns out they're consistent # so I'll just assume s then check both below answer = (ncostheta.real > 0) # convert from numpy boolean to the normal Python boolean answer = bool(answer) # double-check the answer ... can't be too careful! # error_string = ("It's not clear which beam is incoming vs outgoing. Weird" # " index maybe?\n" # "n: " + str(n) + " angle: " + str(theta)) # if answer is True: # assert ncostheta.imag > -100 * EPSILON, error_string # assert ncostheta.real > -100 * EPSILON, error_string # assert (n * cos(theta.conjugate())).real > -100 * EPSILON, error_string # else: # assert ncostheta.imag < 100 * EPSILON, error_string # assert ncostheta.real < 100 * EPSILON, error_string # assert (n * cos(theta.conjugate())).real < 100 * EPSILON, error_string return answer def snell(n_1, n_2, th_1): """ return angle theta in layer 2 with refractive index n_2, assuming it has angle th_1 in layer with refractive index n_1. Use Snell's law. Note that "angles" may be complex!! """ # Important that the arcsin here is scipy.arcsin, not numpy.arcsin! (They # give different results e.g. for arcsin(2).) th_2_guess = sp.arcsin(n_1*np.sin(th_1) / n_2) if is_forward_angle(n_2, th_2_guess): return th_2_guess else: return pi - th_2_guess def list_snell(n_list, th_0): """ return list of angle theta in each layer based on angle th_0 in layer 0, using Snell's law. n_list is index of refraction of each layer. Note that "angles" may be complex!! """ # Important that the arcsin here is scipy.arcsin, not numpy.arcsin! (They # give different results e.g. for arcsin(2).) angles = sp.arcsin(n_list[0]*np.sin(th_0) / n_list) # The first and last entry need to be the forward angle (the intermediate # layers don't matter, see https://arxiv.org/abs/1603.02720 Section 5) if not is_forward_angle(n_list[0], angles[0]): angles[0] = pi - angles[0] if not is_forward_angle(n_list[-1], angles[-1]): angles[-1] = pi - angles[-1] return angles def interface_r(polarization, n_i, n_f, th_i, th_f): """ reflection amplitude (from Fresnel equations) polarization is either "s" or "p" for polarization n_i, n_f are (complex) refractive index for incident and final th_i, th_f are (complex) propegation angle for incident and final (in radians, where 0=normal). "th" stands for "theta". """ if polarization == 's': return ((n_i * np.cos(th_i) - n_f * np.cos(th_f)) / (n_i * np.cos(th_i) + n_f * np.cos(th_f))) else: return ((n_f * np.cos(th_i) - n_i * np.cos(th_f)) / (n_f * np.cos(th_i) + n_i * np.cos(th_f))) def interface_t(polarization, n_i, n_f, th_i, th_f): """ transmission amplitude (frem Fresnel equations) polarization is either "s" or "p" for polarization n_i, n_f are (complex) refractive index for incident and final th_i, th_f are (complex) propegation angle for incident and final (in radians, where 0=normal). "th" stands for "theta". """ if polarization == 's': return 2 * n_i * np.cos(th_i) / (n_i * np.cos(th_i) + n_f * np.cos(th_f)) else: return 2 * n_i * np.cos(th_i) / (n_f * np.cos(th_i) + n_i * np.cos(th_f)) def coh_tmm(pol, n_list, d_list, float th_0, float lam_vac): """ Main "coherent transfer matrix method" calc. Given parameters of a stack, calculates everything you could ever want to know about how light propagates in it. (If performance is an issue, you can delete some of the calculations without affecting the rest.) pol is light polarization, "s" or "p". n_list is the list of refractive indices, in the order that the light would pass through them. The 0'th element of the list should be the semi-infinite medium from which the light enters, the last element should be the semi- infinite medium to which the light exits (if any exits). th_0 is the angle of incidence: 0 for normal, pi/2 for glancing. Remember, for a dissipative incoming medium (n_list[0] is not real), th_0 should be complex so that n0 sin(th0) is real (intensity is constant as a function of lateral position). d_list is the list of layer thicknesses (front to back). Should correspond one-to-one with elements of n_list. First and last elements should be "inf". lam_vac is vacuum wavelength of the light. * R--reflected wave power (as fraction of incident) """ # Convert lists to numpy arrays if they're not already. n_list = array(n_list) d_list = array(d_list, dtype=float) cdef int num_layers num_layers = n_list.size th_list = list_snell(n_list, th_0) # kz is the z-component of (complex) angular wavevector for forward-moving # wave. Positive imaginary part means decaying. kz_list = 2 * np.pi * n_list * np.cos(th_list) / lam_vac # delta is the total phase accrued by traveling through a given layer. # Ignore warning about inf multiplication delta = kz_list * d_list # t_list[i,j] and r_list[i,j] are transmission and reflection amplitudes, # respectively, coming from i, going to j. Only need to calculate this when # j=i+1. (2D array is overkill but helps avoid confusion.) t_list = zeros((num_layers, num_layers), dtype=complex) r_list = zeros((num_layers, num_layers), dtype=complex) for i in range(num_layers-1): t_list[i,i+1] = interface_t(pol, n_list[i], n_list[i+1], th_list[i], th_list[i+1]) r_list[i,i+1] = interface_r(pol, n_list[i], n_list[i+1], th_list[i], th_list[i+1]) # My M is a bit different than Sernelius's, but Mtilde is the same. M_list = zeros((num_layers, 2, 2), dtype=complex) for i in range(1, num_layers-1): M_list[i] = (1/t_list[i,i+1]) * np.dot( make_2x2_array(np.exp(-1j*delta[i]), 0, 0, np.exp(1j*delta[i]), dtype=complex), make_2x2_array(1, r_list[i,i+1], r_list[i,i+1], 1, dtype=complex)) Mtilde = make_2x2_array(1, 0, 0, 1) for i in range(1, num_layers-1): Mtilde = np.dot(Mtilde, M_list[i]) Mtilde = np.dot(make_2x2_array(1, r_list[0,1], r_list[0,1], 1, dtype=complex)/t_list[0,1], Mtilde) # Net complex transmission and reflection amplitudes r = Mtilde[1,0]/Mtilde[0,0] # t = 1/Mtilde[0,0] # # vw_list[n] = [v_n, w_n]. v_0 and w_0 are undefined because the 0th medium # # has no left interface. # vw_list = zeros((num_layers, 2), dtype=complex) # vw = array([[t],[0]]) # vw_list[-1,:] = np.transpose(vw) # for i in range(num_layers-2, 0, -1): # vw = np.dot(M_list[i], vw) # vw_list[i,:] = np.transpose(vw) return np.abs(r)**2