diff --git a/Examples/P_q_vs_a_p.py b/Examples/P_q_vs_a_p.py
new file mode 100644
index 0000000..fd8f17e
--- /dev/null
+++ b/Examples/P_q_vs_a_p.py
@@ -0,0 +1,46 @@
+"""
+Plots Mie and Rayleigh form factors at varying particle radii.
+"""
+
+from scattering import *
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Default Scattering Parameters
+n_p = 2.0  # particle refractive index
+n_s = 1.332  # medium refractive index (water)
+lambda_vac = 685e-9  # wavelength of light in vacuum [meters]
+phi = 0.03  # particle volume fraction
+n_ang = 100  # number of sampled angles
+
+a_p_range = [50e-9, 100e-9, 250e-9, 500e-9, 1000e-9]  # particle radii [meters]
+p_q = np.empty(shape=(n_ang, np.size(a_p_range), 2))  # initialize result array
+
+# Collect Scattering Data
+for i in range(np.size(a_p_range)):
+    _, _, _, _, [q, p_q[:, i, 0], _, _, _, _, _, _] = \
+        mie_scattering(n_p, n_s, a_p_range[i], lambda_vac, phi, n_ang=n_ang, struct='PY')
+    _, _, _, _, [_, p_q[:, i, 1], _, _] = \
+        rayleigh_scattering(n_p, n_s, a_p_range[i], lambda_vac, phi, n_ang=n_ang)
+
+# Generate Colors
+colormap = plt.get_cmap('plasma')
+c = np.empty(shape=(np.size(a_p_range), 3))
+for i in range(np.size(a_p_range)):
+    c[i, :] = colormap.colors[round(256 * i / np.size(a_p_range))]
+
+# Plot
+fig, ax1 = plt.subplots()
+for i in range(np.size(a_p_range)):
+    if i == 0:
+        ax1.plot(q, p_q[:, i, 0] / p_q[0, i, 0], label=r'Mie, $a_p$ = %i nm' % (a_p_range[i] * 1e9), c=c[i])
+        ax1.plot(q, p_q[:, i, 1] / p_q[0, i, 1], linestyle='--',
+                 label=r'Rayleigh, $a_p$ = %i nm' % (a_p_range[i] * 1e9), c=c[i])
+    else:
+        ax1.plot(q, p_q[:, i, 0] / p_q[0, i, 0], label=r'$a_p$ = %i nm' % (a_p_range[i] * 1e9), c=c[i])
+        ax1.plot(q, p_q[:, i, 1] / p_q[0, i, 1], linestyle='--', c=c[i])
+ax1.set(xlabel='q', ylabel='Normalized P(q)', title=r'$n_p$ = %.3f, $n_s$ = %.3f, $\lambda$ = %i nm, $\phi$ = %.2f' % (n_p, n_s, lambda_vac*1e9, phi))
+ax1.legend(loc='upper right')
+ax1.set(ylim=[-0.05, 1.1])
+plt.show()
+
diff --git a/Examples/Results/Intensity Plot.png b/Examples/Results/Intensity Plot.png
new file mode 100644
index 0000000..3c85445
Binary files /dev/null and b/Examples/Results/Intensity Plot.png differ
diff --git a/Examples/Results/Mie vs Rayleigh l_star.png b/Examples/Results/Mie vs Rayleigh l_star.png
new file mode 100644
index 0000000..3ddcc5a
Binary files /dev/null and b/Examples/Results/Mie vs Rayleigh l_star.png differ
diff --git a/Examples/Results/P(q) vs a_p.png b/Examples/Results/P(q) vs a_p.png
new file mode 100644
index 0000000..1f0e54f
Binary files /dev/null and b/Examples/Results/P(q) vs a_p.png differ
diff --git a/Examples/Results/S(q) Models.png b/Examples/Results/S(q) Models.png
new file mode 100644
index 0000000..00e80c4
Binary files /dev/null and b/Examples/Results/S(q) Models.png differ
diff --git a/Examples/Results/Silica Fit.png b/Examples/Results/Silica Fit.png
new file mode 100644
index 0000000..b14b79d
Binary files /dev/null and b/Examples/Results/Silica Fit.png differ
diff --git a/Examples/S_q_models.py b/Examples/S_q_models.py
new file mode 100644
index 0000000..cad8b64
--- /dev/null
+++ b/Examples/S_q_models.py
@@ -0,0 +1,38 @@
+"""
+Plots various structure factor models.
+"""
+
+from scattering import *
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Default Scattering Parameters
+a_p = 200e-9  # meters
+n_p = 1.5  # particle refractive index
+n_s = 1.332  # medium refractive index (water)
+lambda_vac = 685e-9  # wavelength of light in vacuum [meters]
+phi = 0.03  # particle volume fraction
+n_ang = 100  # number of sampled angles
+
+s_q = np.empty(shape=(n_ang, 3))  # initialize result array
+
+# Collect Scattering Data
+_, _, _, _, [q, _, s_q[:, 0], _, _, _, _, _] = mie_scattering(n_p, n_s, a_p, lambda_vac, phi, n_ang=n_ang,
+                                                              struct='PY')
+_, _, _, _, [_, _, s_q[:, 1], _, _, _, _, _] = mie_scattering(n_p, n_s, a_p, lambda_vac, phi, n_ang=n_ang,
+                                                              struct='PYAnnulus', optional_params=[1.2*a_p])
+_, _, _, _, [_, _, s_q[:, 2], _, _, _, _, _] = mie_scattering(n_p, n_s, a_p, lambda_vac, phi, n_ang=n_ang,
+                                                              struct='SHS', optional_params=[0.2, 0.05])
+
+# Plot
+fig, ax1 = plt.subplots()
+ax1.plot(q, np.ones_like(s_q[:, 0]), 'k', label='No Interaction', linewidth=2)
+ax1.plot(q, s_q[:, 0], 'k-.', label='Hard Sphere', linewidth=2)
+ax1.plot(q, s_q[:, 1], 'k:', label='1.2$a_p$ Excluded Annulus', linewidth=2)
+ax1.plot(q, s_q[:, 2], 'k--', label='Sticky Hard Sphere', linewidth=2)
+ax1.set(xlabel=r'$q$', ylabel='S(q)',
+        title=r'$n_p$ = %.3f, $n_s$ = %.3f, $a_p$ = %i nm, $\lambda$ = %i nm, $\phi$ = %.2f'
+              % (n_p, n_s, a_p*1e9, lambda_vac*1e9, phi))
+ax1.legend()
+
+plt.show()
diff --git a/Examples/intensity_plot.py b/Examples/intensity_plot.py
new file mode 100644
index 0000000..6e4c7e4
--- /dev/null
+++ b/Examples/intensity_plot.py
@@ -0,0 +1,56 @@
+"""
+Plots Mie intensities and form factor versus scattering angle.
+"""
+
+from scattering import *
+import matplotlib.pyplot as plt
+
+# Default Scattering Parameters
+a_p = 500e-9  # # particle radius [meters]
+n_p = 1.5  # particle refractive index
+n_s = 1.0  # medium refractive index
+lambda_vac = 685e-9  # wavelength of light in vacuum [meters]
+phi = 0.03  # particle volume fraction
+n_ang = 100  # number of sampled angles
+
+# Collect Data
+theta, i1, i2, _, _ \
+    = mie_scattering(n_p, n_s, a_p, lambda_vac, phi)
+_, i1r, i2r, _, _ \
+    = rayleigh_scattering(n_p, n_s, a_p, lambda_vac, phi)
+theta = np.append(theta, np.pi + theta)
+i1 = np.append(i1, i1[::-1])
+i2 = np.append(i2, i2[::-1])
+i1r = np.append(i1r, i1r[::-1])
+i2r = np.append(i2r, i2r[::-1])
+
+# Generate Colors
+n_colors = 3
+colormap = plt.get_cmap('plasma')
+c = np.empty(shape=(n_colors, 3))
+for i in range(n_colors):
+    c[i, :] = colormap.colors[round(256 * i / n_colors)]
+
+# Plot
+fig, (ax1, ax2) = plt.subplots(1, 2, subplot_kw=dict(projection='polar'))
+ax1.plot(theta, i1, label='Perpendicular Mie', c=c[0])
+ax1.plot(theta, i2, label='Parallel Mie', c=c[2])
+ax1.plot(theta, 0.5*(i1+i2), label='Mie P(q)', c=c[1])
+# ax1.plot(theta, i1r, label='Perpendicular Rayleigh', c=c[0], linestyle='--')
+# ax1.plot(theta, i2r, label='Parallel Rayleigh', c=c[2], linestyle='--')
+# ax1.plot(theta, 0.5*(i1r+i2r), label='Rayleigh P(q)', c=c[1], linestyle='--')
+ax1.set_title("Intensity")
+ax1.legend()
+
+ax2.plot(theta, np.log10(i1 + 1), label='Perpendicular Mie', c=c[0])
+ax2.plot(theta, np.log10(i2 + 1), label='Parallel Mie', c=c[2])
+ax2.plot(theta, np.log10(0.5*(i1+i2) + 1), label='Mie P(q)', c=c[1])
+# ax2.plot(theta, np.log10(i1r + 1), label='Perpendicular Rayleigh', c=c[0], linestyle='--')
+# ax2.plot(theta, np.log10(i2r + 1), label='Parallel Rayleigh', c=c[2], linestyle='--')
+# ax2.plot(theta, np.log10(0.5*(i1r+i2r) + 1), label='Rayleigh P(q)', c=c[1], linestyle='--')
+ax2.set_title("Log Intensity")
+ax2.legend()
+
+plt.suptitle(r'$n_p$ = %.3f, $n_s$ = %.3f, $a_p$ = %i nm, $\lambda$ = %i nm, $\phi$ = %.2f'
+             % (n_p, n_s, a_p*1e9, lambda_vac*1e9, phi))
+plt.show()
diff --git a/Examples/mie_vs_rayleigh.py b/Examples/mie_vs_rayleigh.py
new file mode 100644
index 0000000..3c11c34
--- /dev/null
+++ b/Examples/mie_vs_rayleigh.py
@@ -0,0 +1,58 @@
+"""
+Plots l* for varying input parameters for both Mie and Rayleigh scattering.
+"""
+
+from scattering import *
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Default Scattering Parameters
+n_p = 1.8  # particle refractive index
+n_s = 1.332  # medium refractive index (water)
+a_p = 250e-9  # particle radii [meters]
+lambda_vac = 685e-9  # wavelength of light in vacuum [meters]
+phi = 0.03  # particle volume fraction
+
+# Generate Parameter Ranges
+n_range = 100  # number of data points between values
+phi_range = np.linspace(0.01, 0.15, n_range)
+a_p_range = np.linspace(150e-9/2, 1000e-9/2, n_range)
+n_p_range = np.linspace(1.6, 2, n_range)
+
+# Initialize Result Array
+l_star = np.empty(shape=(n_range, 3, 2))
+
+# Collect Scattering Data
+for i in range(n_range):
+    _, _, _, l_star[i, 0, 0], _ \
+        = mie_scattering(n_p, n_s, a_p, lambda_vac, phi_range[i])
+    _, _, _, l_star[i, 1, 0], _ \
+        = mie_scattering(n_p, n_s, a_p_range[i], lambda_vac, phi)
+    _, _, _, l_star[i, 2, 0], _ \
+        = mie_scattering(n_p_range[i], n_s, a_p, lambda_vac, phi)
+
+    _, _, _, l_star[i, 0, 1], _ \
+        = rayleigh_scattering(n_p, n_s, a_p, lambda_vac, phi_range[i])
+    _, _, _, l_star[i, 1, 1], _ \
+        = rayleigh_scattering(n_p, n_s, a_p_range[i], lambda_vac, phi)
+    _, _, _, l_star[i, 2, 1], _ \
+        = rayleigh_scattering(n_p_range[i], n_s, a_p, lambda_vac, phi)
+
+# Plot Data
+fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
+ax1.plot(phi_range, l_star[:, 0, 0]*1e6, 'k', label='Mie', linewidth=2)
+ax1.plot(phi_range, l_star[:, 0, 1]*1e6, 'k--', label='Rayleigh', linewidth=2)
+ax1.legend()
+ax1.set(xlabel=r'$\phi$', ylabel='l* [µm]')
+
+ax2.plot(a_p_range*1e9, l_star[:, 1, 0]*1e6, 'k', label='Mie', linewidth=2)
+ax2.plot(a_p_range*1e9, l_star[:, 1, 1]*1e6, 'k--', label='Rayleigh', linewidth=2)
+ax2.set(xlabel=r'$a_p$ [nm]')
+
+ax3.plot(n_p_range, l_star[:, 2, 0]*1e6, 'k', label='Mie', linewidth=2)
+ax3.plot(n_p_range, l_star[:, 2, 1]*1e6, 'k--', label='Rayleigh', linewidth=2)
+ax3.set(xlabel=r'$n_p$')
+
+# title with default parameters listed
+plt.suptitle(r'$n_p$ = %.3f, $n_s$ = %.3f, $a_p$ = %i nm, $\lambda$ = %i nm, $\phi$ = %.2f' % (n_p, n_s, a_p*1e9, lambda_vac*1e9, phi))
+plt.show()
diff --git a/Examples/silica_data_fit.py b/Examples/silica_data_fit.py
new file mode 100644
index 0000000..f6e267d
--- /dev/null
+++ b/Examples/silica_data_fit.py
@@ -0,0 +1,55 @@
+"""
+Generates l* values versus particle volume concentration for silica particles with varying form factor models.
+Plots model versus experimental data (courtesy of Qi Li).
+"""
+
+from scattering import *
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Scattering Parameters
+n_p = 1.446  # particle refractive index (silica)
+n_s = 1.332  # medium refractive index (water)
+a_p = 345e-9/2  # particle radii [meters]
+lambda_vac = 685e-9  # wavelength of light in vacuum [meters]
+
+# Generate Parameter Range
+n_range = 100  # number of data points between values
+phi_range = np.linspace(0.03, 0.15, n_range)
+
+# Initialize Result Array
+l_star_phi = np.empty(shape=(n_range, 4))
+
+# Silica Data (courtesy of Qi Li)
+silica_data = np.array([[0.065, 0.067, 0.069, 0.072, 0.075, 0.135],
+                        [0.000230, 0.000225, 0.000216, 0.000215, 0.000205, 0.000130]])
+
+# Collect Scattering Data
+for i in range(n_range):
+    _, _, _,  l_star_phi[i, 0], _ \
+        = mie_scattering(n_p, n_s, a_p, lambda_vac, phi_range[i])
+    _, _, _, l_star_phi[i, 1], _ \
+        = mie_scattering(n_p, n_s, a_p, lambda_vac, phi_range[i], struct='PY')
+    _, _, _, l_star_phi[i, 2], _ \
+        = mie_scattering(n_p, n_s, a_p, lambda_vac, phi_range[i], struct='SHS', optional_params=[0.2, 0.05])
+    _, _, _, l_star_phi[i, 3], _ \
+        = mie_scattering(n_p, n_s, a_p, lambda_vac, phi_range[i], struct='PYAnnulus', optional_params=[1.2*a_p])
+
+# Generate Colors
+n_colors = 4
+colormap = plt.get_cmap('plasma')
+c = np.empty(shape=(n_colors, 3))
+for i in range(n_colors):
+    c[i, :] = colormap.colors[round(256 * i / n_colors)]
+
+# Plot Data
+fig, ax1 = plt.subplots()
+ax1.plot(phi_range, l_star_phi[:, 0]*1e6, label='No Interactions', linewidth=2, c=c[0])
+ax1.plot(phi_range, l_star_phi[:, 1]*1e6, linestyle='--', label='Hard Sphere', linewidth=2, c=c[1])
+ax1.plot(phi_range, l_star_phi[:, 2]*1e6, linestyle=':', label='Sticky Hard Sphere', linewidth=2, c=c[2])
+ax1.plot(phi_range, l_star_phi[:, 3]*1e6, linestyle='-.', label='1.2$a_p$ Excluded Annulus', linewidth=2, c=c[3])
+ax1.plot(silica_data[0, :], silica_data[1, :]*1e6, 'k.', label='Silica Data', markersize=10)
+ax1.legend()
+ax1.set(xlabel=r'$\phi$', ylabel='l* [µm]')
+
+plt.show()
diff --git a/scattering.py b/scattering.py
new file mode 100644
index 0000000..50702e7
--- /dev/null
+++ b/scattering.py
@@ -0,0 +1,421 @@
+"""
+Function package for generating theoretical light scattering models.
+
+Originally composed by Nicholas Sbalbi (nsbalbi@umass.edu) during the 2021 CHARM REU Program for the Furst Group in the
+University of Delaware Department of Chemical and Biomolecular Engineering. Last updated 08/13/2021
+"""
+
+import numpy as np
+import scipy.special as sps
+import sys
+from scattering_test import *
+
+
+def mie_pi_tau(theta, alpha):
+    """
+    Calculate the functions pi and tau which appear in the Mie scattering calculation.
+
+    These functions match those described by van de Hulst in Light Scattering by Small Particles, 1981, and are
+    calculated using the method outlined by Kerker in The Scattering of Light and Other Electromagnetic Radiation, 1969.
+
+    Args:
+        theta: the desired scattering angles
+        alpha: the size parameter of the spheres (2 * pi * radius / wavelength)
+
+    Returns:
+        pi: the first function
+        tau: the second function
+    """
+    # Constants
+    cos_theta = np.cos(theta)  # pre-calculated for efficiency
+    sin_theta = np.sin(theta)  # "
+    n_ang = np.size(theta)  # number of scattering angles
+    n_terms = np.rint(alpha + 4.05*alpha**0.3333 + 2).astype(int)  # number of terms in series for coeff calculation
+
+    # Initialize arrays
+    pi = np.empty(shape=(n_terms, n_ang))
+    pi_prime = np.empty(shape=(n_terms, n_ang))
+    tau = np.empty(shape=(n_terms, n_ang))
+
+    # Caluclation of pi and tau via Kerker's method
+    for i in range(n_terms):  # loop over rows
+        pi[i, :] = -sps.lpmv(1, i + 1, cos_theta) / sin_theta  # generates 1st order Legendre polynomials
+        if i == 0:
+            pi_prime[i, :] = np.zeros(shape=n_ang)
+            tau[i, :] = cos_theta * pi[i, :]
+        elif i == 1:
+            pi_prime[i, :] = 3 * np.ones(shape=n_ang)
+            tau[i, :] = cos_theta * pi[i, :] - 3 * sin_theta ** 2
+        else:
+            for j in range(n_ang):  # loop over columns
+                pi_prime[i, j] = (2 * i + 1) * pi[i - 1, j] + pi_prime[i - 2, j]
+            tau[i, :] = cos_theta * pi[i, :] - sin_theta ** 2 * pi_prime[i, :]
+
+    return pi, tau
+
+
+def mie_a_b(alpha, m):
+    """
+    Calculate the Mie coefficients for Mie scattering by spheres.
+
+    The size of the arrays is based on the recommended by Barber and Hill in Light Scattering by Particles:
+    Computational Methods, 1990. The coefficients are calculated using the method outlined by Kerker in
+    The Scattering of Light and Other Electromagnetic Radiation, 1969.
+
+    Args:
+        alpha: the size parameter of the spheres (2 * pi * radius / wavelength)
+        m: the ratio between the refractive indices of the spheres and the medium
+
+    Returns:
+        a: first set of Mie coefficients
+        b: second set of Mie coefficients
+    """
+    # Constants
+    beta = m * alpha
+    n_terms = np.rint(alpha + 4.05*alpha**0.3333 + 2).astype(int)  # number of terms in series for coeff calculation
+    x = np.arange(n_terms + 2) + 0.5  # two extra values for low and high bounds
+
+    # Calculation of coefficients via Kerker's method
+    psi = np.empty(shape=(2, n_terms + 2))  # add two extra terms for left and right bounds
+    chi = np.empty(shape=(2, n_terms + 2))
+    psi[0, :] = (np.pi * alpha / 2) ** 0.5 * sps.jv(x, alpha)  # generates Bessel function values of the first kind
+    psi[1, :] = (np.pi * beta / 2) ** 0.5 * sps.jv(x, beta)  # generates Bessel function values of the second kind
+    chi[0, :] = -(np.pi * alpha / 2) ** 0.5 * sps.yv(x, alpha)
+    chi[1, :] = -(np.pi * beta / 2) ** 0.5 * sps.yv(x, beta)
+
+    zeta = psi + 1j * chi  # conversion to complex numbers
+
+    psi_prime = np.empty(shape=(2, n_terms))
+    zeta_prime = np.empty(shape=(2, n_terms), dtype=np.complex64)
+    psi_prime[0, :] = 0.5 * psi[0, :-2] - 0.5 * psi[0, 2:] + 0.5 / alpha * psi[0, 1:-1]
+    psi_prime[1, :] = 0.5 * psi[1, :-2] - 0.5 * psi[1, 2:] + 0.5 / beta * psi[1, 1:-1]
+    zeta_prime[0, :] = 0.5 * zeta[0, :-2] - 0.5 * zeta[0, 2:] + 0.5 / alpha * zeta[0, 1:-1]
+    zeta_prime[1, :] = 0.5 * zeta[1, :-2] - 0.5 * zeta[1, 2:] + 0.5 / beta * zeta[1, 1:-1]
+
+    a = (psi_prime[1, :] * psi[0, 1:-1] - m * psi_prime[0, :] * psi[1, 1:-1]) / \
+        (psi_prime[1, :] * zeta[0, 1:-1] - m * zeta_prime[0, :] * psi[1, 1:-1])
+    b = (m * psi_prime[1, :] * psi[0, 1:-1] - psi_prime[0, :] * psi[1, 1:-1]) / \
+        (m * psi_prime[1, :] * zeta[0, 1:-1] - zeta_prime[0, :] * psi[1, 1:-1])
+
+    return a, b
+
+
+def mie_s1_s2(theta, alpha, m):
+    """
+    Calculate the scattering amplitude functions for Mie scattering by spheres.
+
+    Args:
+        theta: the angles at which scattering amplitudes are calculated
+        alpha: the size parameter of the spheres (2 * pi * radius / wavelength)
+        m: the ratio between the refractive indices of the spheres and the medium
+
+    Returns:
+        s1: perpendicular amplitude function at corresponding theta values
+        s2: parallel amplitude function at corresponding theta values
+        a: first set of Mie coefficients
+        b: second set of Mie coefficients
+    """
+    # Constants
+    n_ang = np.size(theta)
+    n_terms = np.rint(alpha + 4.05*alpha**0.3333 + 2).astype(int)  # number of terms in series for coeff calculation
+
+    # Calculation of coefficients using other functions
+    pi, tau = mie_pi_tau(theta, alpha)
+    a, b = mie_a_b(alpha, m)
+
+    # Calculation of scattering amplitude functions
+    s1 = np.reshape(np.tile(a, n_ang), (n_terms, n_ang), order='F') * pi \
+        + np.reshape(np.tile(b, n_ang), (n_terms, n_ang), order='F') * tau  # perpendicular amplitude function
+    s2 = np.reshape(np.tile(a, n_ang), (n_terms, n_ang), order='F') * tau \
+        + np.reshape(np.tile(b, n_ang), (n_terms, n_ang), order='F') * pi  # parallel amplitude function
+    temp = np.arange(n_terms) + 1  # factor multiplied by each term (row)
+    temp = (2 * temp + 1) / (temp ** 2 + temp)
+    s1 = s1 * np.reshape(np.tile(temp, n_ang), (n_terms, n_ang), order='F')
+    s2 = s2 * np.reshape(np.tile(temp, n_ang), (n_terms, n_ang), order='F')
+    s1 = np.sum(s1, axis=0)  # sum along angle (column)
+    s2 = np.sum(s2, axis=0)
+
+    return s1, s2, a, b
+
+
+def percus_yevick(qa, phi):
+    """
+    Calculate the structure factor for hard spheres using the Percus-Yevick
+    closure of the Ornstein-Zernicke equation.
+
+    Args:
+        qa: the scattering vector scaled by the sphere radius at the desired scattering angles
+        phi: volume fraction of spheres in solution
+
+    Returns:
+        s_q: the structure factor at the corresponding scattering angles
+    """
+    qa2 = 2 * qa  # twice qa, for efficiency/conciseness
+    lambda_1 = (1 + 2 * phi) ** 2 / (1 - phi) ** 4
+    lambda_2 = -(1 + phi / 2) ** 2 / (1 - phi) ** 4
+    c_1 = lambda_1 / qa2 ** 3 * (np.sin(qa2) - qa2 * np.cos(qa2))
+    c_2 = -6 * phi * lambda_2 / qa2 ** 4 * (qa2 ** 2 * np.cos(qa2) - 2 * qa2 * np.sin(qa2) - 2 * np.cos(qa2) + 2)
+    c_3 = -phi * lambda_1 / 2 / qa2 ** 6 * \
+        (qa2 ** 4 * np.cos(qa2) - 4 * qa2 ** 3 * np.sin(qa2) - 12 * qa2 ** 2 * np.cos(qa2)
+         + 24 * qa2 * np.sin(qa2) + 24 * np.cos(qa2) - 24)
+
+    s_q = np.divide(1, 1 + 24 * phi * (c_1 + c_2 + c_3))
+    return s_q
+
+
+def percus_yevick_annulus(q, a_p, a_e, phi):
+    """
+    Calculate the structure factor for hard spheres with an exluded annulus
+    using the Percus-Yevick closure of the Ornstein-Zernicke equation.
+
+    Args:
+        q: scattering vector at desired scattering angles
+        a_p: sphere radius (excluding annulus)
+        a_e: effective sphere radius (including annulus)
+        phi: volume fraction of spheres in solution
+
+    Returns:
+        s_q: the structure factor at the corresponding scattering angles
+    """
+    phi_e = phi * (a_e / a_p) ** 3  # effective volume fraction
+
+    qa2 = 2 * q * a_e  # twice qa, for efficiency/conciseness
+    lambda_1 = (1 + 2 * phi_e) ** 2 / (1 - phi_e) ** 4
+    lambda_2 = -(1 + phi_e / 2) ** 2 / (1 - phi_e) ** 4
+    c_1 = lambda_1 / qa2 ** 3 * (np.sin(qa2) - qa2 * np.cos(qa2))
+    c_2 = -6 * phi_e * lambda_2 / qa2 ** 4 * (qa2 ** 2 * np.cos(qa2) - 2 * qa2 * np.sin(qa2) - 2 * np.cos(qa2) + 2)
+    c_3 = -phi_e * lambda_1 / 2 / qa2 ** 6 * \
+        (qa2 ** 4 * np.cos(qa2) - 4 * qa2 ** 3 * np.sin(qa2) - 12 * qa2 ** 2 * np.cos(qa2)
+         + 24 * qa2 * np.sin(qa2) + 24 * np.cos(qa2) - 24)
+
+    s_q = np.divide(1, 1 + 24 * phi_e * (c_1 + c_2 + c_3))
+    return s_q
+
+
+def sticky_hard_sphere(q, a_p, phi, tau=0.2, epsilon=0.05):
+    """
+    Calculate the structure factor for hard spheres with square-well interparticle
+    potential. Based on the method of Rao et. al. in "A new interpretation of the
+    sticky hard sphere model", J. Chem. Phys. 95(12), 9186-9190 (1991).
+
+    Args:
+        q: scattering vector at desired scattering angles
+        a_p: sphere radius (excluding annulus)
+        phi: volume fraction of spheres in solution
+        tau: float, "stickiness" parameter within {0, 1}
+        epsilon: float, perturbation parameter within {0, 0.1}
+
+    Returns:
+        s_q: the structure factor at the corresponding scattering angles
+    """
+    eta = phi / (1 - epsilon) ** 3
+    eta1 = 1 - eta
+    sigma = 2 * a_p
+    a = sigma / (1 - epsilon)
+    kappa = q * a
+
+    # Quadratic for Lambda
+    a_q = eta/12
+    b_q = -tau - eta / eta1
+    c_q = (1 + eta / 2) / eta1 ** 2
+
+    disc = b_q ** 2 - 4 * a_q * c_q  # discriminant
+    if disc < 0:  # error if no real roots
+        print('Error: No Real Roots within SHS Model')
+        return -1*np.ones_like(q)
+    lamb = (-b_q + np.sqrt(disc)) / (2 * a_q)
+    lamb2 = (-b_q - np.sqrt(disc)) / (2 * a_q)
+    if lamb2 < lamb:
+        lamb = lamb2  # take smaller root, larger is unphysical
+
+    # Alpha and Beta
+    mu = lamb * eta * eta1
+    if mu > 1 + 2 * eta:  # model condition set in paper
+        print('Error: Unphysical SHS Model Result')
+        return -1*np.ones_like(q)
+    alpha = (1 + 2 * eta - mu) / eta1 ** 2
+    beta = (-3 * eta + mu) / 2 / eta1 ** 2
+
+    # S_q
+    k2 = kappa ** 2
+    k3 = kappa ** 3
+    ksin = np.sin(kappa)
+    kcos = np.cos(kappa)
+    A = 1 + 12 * eta * (alpha * (ksin - kappa * kcos) / k3 + beta * (1 - kcos) / k2 - lamb / 12 * ksin / kappa)
+    B = 12 * eta * (alpha * (1 / 2 / kappa - ksin / k2 + (1 - kcos) / k3)
+                    + beta * (1 / kappa - ksin / k2) - lamb / 12 * (1 - kcos) / kappa)
+    s_q = 1 / (A ** 2 + B ** 2)
+
+    return s_q
+
+
+def mie_scattering(n_p, n_s, a_p, lambda_vac, phi, n_ang=100, struct='none', optional_params=None):
+    """
+    Calculate intensities and statistics for Mie scattering by spheres.
+
+    Methods and algorithms are those outlined by Kerker in
+    The Scattering of Light and Other Electromagnetic Radiation, 1969
+
+    Args:
+        n_p: float, refractive index of the spheres
+        n_s: float, refractive index of the medium
+        a_p: float, sphere radius [m]
+        lambda_vac: float, wavelength of light in vacuum [m]
+        phi: float, volume fraction of particles in solution
+        n_ang: int, number of angles calculated (Default = 100)
+        struct: str, choice of equation used for structure factor
+            - 'none': S(q) = 1
+            - 'PY': S(q) is calculated using the Percus-Yevick approximation
+            - 'PYAnnulus': S(q) is calculated using the Percus-Yevick approximation with an excluded annulus
+            - 'SHS': S(q) is calculated via a perturbative sticky hard sphere solution of the Percus-Yevick closure
+        optional_params: list of optional parameters required for some structure factor calculations
+            - 'PYAnnulus': [a_e]
+                - a_e: float, effective sphere radius (including excluded annulus)
+            - 'SHS': [tau, epsilon]
+                - tau: float, "stickiness" parameter within {0, 1}
+                - epsilon: float, perturbation parameter within {0, 0.1}
+
+    Returns:
+        theta: the angles at which scattering intensities are calculated
+        i1: perpendicular scattering intensity at corresponding scattering angles
+        i2: parallel scattering intensity
+        l_star: photon mean free path
+        stats: list of other output scattering statistics
+            - q: array of scattering vectors corresponding with theta
+            - p_q: form factor at corresponding theta/q values
+            - s_q: structure factor at corresponding theta/q values
+            - l: scattering mean free path
+            - q_scat: scattering efficiency factor
+            - q_ext: extinction efficiency factor
+            - c_scat: scattering cross section
+            - c_ext: extinction cross section
+    """
+    # Initial Constants
+    m = n_p / n_s  # refractive index ratio
+    lambda_med = lambda_vac / n_s  # wavelength in medium
+    k_0 = 2 * np.pi / lambda_med  # scattering wavevector
+    alpha = k_0 * a_p
+    rho = phi / 4 * 3 / np.pi / a_p ** 3  # number density of scatterers
+    n_terms = np.rint(alpha + 4.05*alpha**0.3333 + 2).astype(int)  # number of terms in series for coeff calculation
+
+    theta = np.linspace(0.001, 0.999 * np.pi, n_ang)  # angles (ends are not absolute to prevent zeroing out)
+    qa = a_p * 2 * k_0 * np.sin(theta/2)  # radius times scattering vector
+
+    # Intensity Calculation
+    s1, s2, a, b = mie_s1_s2(theta, alpha, m)
+
+    i1 = np.real(s1 * np.conjugate(s1))  # perpendicular intensity
+    i2 = np.real(s2 * np.conjugate(s2))  # parallel intensity
+
+    # Scattering Statistics Calculation
+    temp = 2 * np.arange(n_terms) + 3  # factor multiplied by each term (row)
+    # scattering efficiency factor
+    q_scat = 2 / alpha ** 2 * np.sum(temp * np.real(a * np.conjugate(a) + b * np.conjugate(b)))
+    q_ext = 2 / alpha ** 2 * np.sum(temp * np.real(a + b))  # extinction efficiency factor
+    c_scat = q_scat * 3.14158 * a_p ** 2  # scattering cross section
+    c_ext = q_ext * 3.14158 * a_p ** 2  # extinction cross section
+
+    # Structure and Form Factor Calculation
+    # structure factor
+    if struct == 'none':
+        s_q = np.ones_like(theta)
+    elif struct == 'PY':
+        s_q = percus_yevick(qa, phi)
+    elif struct == 'PYAnnulus':
+        s_q = percus_yevick_annulus(qa / a_p, a_p, optional_params[0], phi)
+    elif struct == 'SHS':
+        s_q = sticky_hard_sphere(qa / a_p, a_p, phi, optional_params[0], optional_params[1])
+    else:
+        sys.exit('Error: Unexpected Input for Choice of Structure Factor Equation')
+    p_q = (i1 + i2) / 2  # form factor (average intensity)
+
+    # l Calculation
+    integrand = qa * p_q * s_q
+    l = (k_0 ** 4 * a_p ** 2) / \
+        (2 * np.pi * rho * 0.5 * np.sum((integrand[:-1] + integrand[1:]) * (qa[1:] - qa[:-1])))
+
+    # l* Calculation
+    integrand = qa ** 3 * p_q * s_q
+    l_star = (k_0 ** 6 * a_p ** 4) / \
+        (np.pi * rho * 0.5 * np.sum((integrand[:-1] + integrand[1:]) * (qa[1:] - qa[:-1])))
+
+    # Gather Misc Statistics
+    q = 2 * k_0 * np.sin(theta / 2)  # scattering vector
+    stats = [q, p_q, s_q, l, q_scat, q_ext, c_scat, c_ext]
+
+    return theta, i1, i2, l_star, stats
+
+
+def rayleigh_scattering(n_p, n_s, a_p, lambda_vac, phi, n_ang=100, struct='none', optional_params=None):
+    """
+    Calculate intensities and statistics for Rayleigh scattering by spheres.
+
+    Args:
+        n_p: float, refractive index of the spheres
+        n_s: float, refractive index of the medium
+        a_p: float, sphere radius [m]
+        lambda_vac: float, wavelength of light in vacuum [m]
+        phi: float, volume fraction of particles in solution
+        n_ang: int, number of angles calculated (Default = 100)
+        struct: str, choice of equation used for structure factor
+            - 'none': S(q) = 1
+            - 'PY': S(q) is calculated using the Percus-Yevick approximation
+            - 'PYAnnulus': S(q) is calculated using the Percus-Yevick approximation with an excluded annulus
+        optional_params: list of optional parameters required for some structure factor calculations
+            - 'PYAnnulus': [a_e]
+                - a_e: float, effective sphere radius (including excluded annulus)
+
+    Returns:
+        theta: the angles at which scattering intensities are calculated
+        i1: perpendicular scattering intensity at corresponding scattering angles
+        i2: parallel scattering intensity
+        l_star: photon mean free path
+        stats: list of other output scattering statistics
+            - q: array of scattering vectors corresponding with theta
+            - p_q: form factor at corresponding theta/q values
+            - s_q: structure factor at corresponding theta/q values
+            - l: scattering mean free path
+    """
+    # Initial Constants
+    m = n_p / n_s  # refractive index ratio
+    lambda_med = lambda_vac / n_s  # wavelength in medium
+    k_0 = 2 * np.pi / lambda_med  # scattering wavevector
+    rho = phi / 4 * 3 / np.pi / a_p ** 3  # number density of scatterers
+    v = 4 / 3 * np.pi * a_p ** 3  # sphere volume
+
+    theta = np.linspace(0.001, 0.999 * np.pi, n_ang)  # angles (ends are not absolute to prevent zeroing out)
+    cos_theta = np.cos(theta)
+    qa = a_p * 2 * k_0 * np.sin(theta/2)  # radius times scattering vector
+
+    # Intensity Calculation
+    i1 = k_0 ** 6 * v ** 2 * ((m - 1) / 2 / np.pi) ** 2 \
+        * np.divide(3 * (np.sin(qa) - qa * np.cos(qa)), qa ** 3) ** 2  # perpendicular Rayleigh  intensity
+    i2 = i1 * cos_theta ** 2  # parallel Rayleigh  intensity
+
+    # Form and Structure Factor Calculation
+    # structure factor
+    if struct == 'none':
+        s_q = np.ones_like(theta)
+    elif struct == 'PY':
+        s_q = percus_yevick(qa, phi)
+    elif struct == 'PYAnnulus':
+        s_q = percus_yevick_annulus(qa / a_p, a_p, optional_params[0], phi)
+    else:
+        sys.exit('Error: Unexpected Input for Choice of Structure Factor Equation')
+    p_q = (i1 + i2) / 2  # form factor (average intensity)
+
+    # l Calculation
+    integrand = qa * p_q * s_q
+    l = (k_0 ** 4 * a_p ** 2) / \
+        (2 * np.pi * rho * 0.5 * np.sum((integrand[:-1] + integrand[1:]) * (qa[1:] - qa[:-1])))
+
+    # l* Calculation
+    integrand = qa ** 3 * p_q * s_q
+    l_star = (k_0 ** 6 * a_p ** 4) / \
+        (np.pi * rho * 0.5 * np.sum((integrand[:-1] + integrand[1:]) * (qa[1:] - qa[:-1])))
+
+    # Gather Misc Statistics
+    q = 2 * k_0 * np.sin(theta / 2)  # scattering vector
+    stats = [q, p_q, s_q, l]
+
+    return theta, i1, i2, l_star, stats