A Python program for showing plots of the thickness of bubbles based on their color patterns.

bubbles.py

# Computes bubble thickness refraction colormaps according to the paper
# <https://www.physics.mun.ca/~cdeacon/publications/Soap%20Bubbles%20-%20AJP%20Oct%202011.pdf>.
# License: MIT
# For this to run, matplotlib and numpy must be installed.

import matplotlib.pyplot as plt
import numpy as np
from math import pi

# Good enough; we only have one instantiation anyways
class Specs:
	pass

# List of bubble specifications.
specs = Specs()

specs.diameter = 14
specs.height = 7
specs.camera_dist = 25
specs.index_of_refrac = 1.3
specs.num_colors = 70
specs.max_thick = 1000
specs.num_thick = 40
specs.red_wl = 650
specs.green_wl = 575
specs.blue_wl = 475

specs.radius = radius = (specs.diameter ** 2 / 4 + specs.height ** 2) / (2 * specs.height)

def compute_colormap(specs, wavelength):
	""" Compute a numpy 2D array of the R_lambda color values. """
	# Base array of x-distances to work off of. `step` is added to the second range since
	# np.arange is [start, end) and we want [start, end].
	distances = np.arange(0, specs.diameter / 2, specs.diameter / 2 / specs.num_colors)

	# Angle calculations for alpha, beta, and theta, respectively
	center_to_dist = np.arcsin(distances / specs.radius)
	camera_to_dist = np.arctan(distances /
			(specs.radius - specs.radius * np.cos(center_to_dist) + specs.height))

	refrac_angle = center_to_dist + camera_to_dist

	# Calculate phi, the phase difference, for a variety of thicknesses. The thicknesses are
	# repeated into a 2D array for multiplication purposes.
	thicknesses = np.arange(0, specs.max_thick, specs.max_thick / specs.num_thick)
	thicknesses = thicknesses.reshape(1, specs.num_thick).repeat(specs.num_colors, axis=0)

	sqrt_part = np.sqrt(specs.index_of_refrac ** 2 - np.sin(refrac_angle) ** 2)

	phase_diff = thicknesses * (4 * pi / wavelength) * sqrt_part.reshape(specs.num_colors, 1)

	# Calculate R_s and R_p.
	amp_refrac_perp = (np.cos(refrac_angle) - sqrt_part) / (np.cos(refrac_angle) + sqrt_part)
	amp_refrac_par = ((specs.index_of_refrac ** 2 * np.cos(refrac_angle) - sqrt_part) /
			(specs.index_of_refrac ** 2 * np.cos(refrac_angle) + sqrt_part))

	# Finally, compute R_lambda for the final result.
	amp_refrac_perp = amp_refrac_perp.reshape(specs.num_colors, 1)
	amp_refrac_par  = amp_refrac_par .reshape(specs.num_colors, 1)

	cos_phase_diff = np.cos(phase_diff)
	frac = lambda amp: amp ** 2 * (1 - cos_phase_diff) / (1 + amp ** 4 - 2 * amp ** 2 * cos_phase_diff)
	return frac(amp_refrac_perp) + frac(amp_refrac_par)

def set_up_plot(plot, title, data, show_labels = False):
	plot.imshow(data, interpolation="none")
	plot.set_title(title)

	if show_labels:
		plot.set_xlabel("Thickness (µm)")
		plot.set_ylabel("x Distance (cm)")

	plot.set_xticks(np.arange(0, specs.num_thick, specs.num_thick / 2))
	plot.set_xticklabels(np.arange(0, specs.max_thick, specs.max_thick / 2))

	plot.set_yticks(np.arange(0, specs.num_colors, specs.num_colors / (specs.diameter / 2)))
	plot.set_yticklabels(np.arange(0, specs.diameter / 2, 1))

def display_colormap(red, green, blue, specs):
	""" Display the created colormap with matplotlib """
	empty = np.zeros(red.shape)

	figure, (rgb_plot, red_plot, green_plot, blue_plot) = plt.subplots(1, 4, figsize = (12, 6))
	figure.suptitle("Bubble Thickness Based on Color Matching")

	set_up_plot(rgb_plot, "Full RGB", np.stack((red, green, blue), axis=2), True)
	set_up_plot(red_plot, "Red", np.stack((red, empty, empty), axis=2))
	set_up_plot(green_plot, "Green", np.stack((empty, green, empty), axis=2))
	set_up_plot(blue_plot, "Blue", np.stack((empty, empty, blue), axis=2))

	figure.tight_layout()
	plt.show()

display_colormap(compute_colormap(specs, specs.red_wl), compute_colormap(specs, specs.green_wl),
		compute_colormap(specs, specs.blue_wl), specs)