In its simplest form, we can use a Green's function approach and just evaluate the solution of a single line source. Based on the analytical solution for the 3D Helmholtz equation (as can be read in this reference), we can write
$$ f(r) = \frac{e^{ikr}}{4πr} $$We can write some straightforward code to evaluate this on a grid, that I choose to be the same than in the tweet.
import numpy as np
wavelength = 1.
n_lambda = 20.
n_points = 401
dx = n_lambda * wavelength / (n_points - 1)
k = 2 * np.pi / wavelength
print(f"number of points per wavelenght: {wavelength/dx}")
x = np.linspace(-n_lambda//2 * wavelength, n_lambda//2 * wavelength, num=n_points)
y = np.linspace(0, n_lambda * wavelength, num=n_points)
X, Y = np.meshgrid(x, y)
Now let's write a function that allows us to compute the field amplitude.
import numba as nb
@nb.jit()
def compute_amplitude_from_source_point(source_location, X, Y, k, phase=0):
x0, y0 = source_location
R = np.sqrt((X - x0)**2 + (Y - y0)**2)
amp = np.exp(1j * (k * R + phase))
return np.real(amp)
amp = compute_amplitude_from_source_point((0., 0.), X, Y, k)
import holoviews as hv
hv.extension('matplotlib', logo=False)
hv.output(holomap='scrubber')
FIG_OPTS = hv.opts(aspect=1.05, fig_inches=5)
def plot(amp):
return hv.Image(amp[::-1], bounds=[-10, 0, 10, 20]).opts(colorbar=True, cmap='seismic').redim.label(x='x (wavelengths)', y='y (wavelengths)', z='amplitude').opts(FIG_OPTS)
plot(amp)
To compute the field produced by a set of sources on a grid, we can just add up the resulting fields. Let's write some helper functions to do that.
def make_symmetric_point_source(n_one_side, dx, wavelength_solid = 2.5 / 1.5):
pos = np.arange(-n_one_side, n_one_side + 1) * dx
phases = 2 * np.pi / wavelength_solid * pos
return np.c_[pos, np.zeros_like(pos)], phases
line_source, phases = make_symmetric_point_source(5, wavelength / 5.)
line_source, np.cos(phases)
@nb.jit()
def compute_amplitude_from_several_points(source_locations, phases, X, Y, k):
field = np.zeros_like(X)
for i in range(len(source_locations)):
source_location, phase = source_locations[i], phases[i]
field += compute_amplitude_from_source_point(source_location, X, Y, k, phase)
return field
field = compute_amplitude_from_several_points(line_source, phases, X, Y, k)
def plot_field(field, line_source, phases):
SCATTER_OPTS = hv.opts(color='z', cmap='gray', padding=0.05, s=100)
return (plot(field) * hv.Scatter(np.c_[line_source, np.cos(phases)], vdims=['y', 'z']).opts(SCATTER_OPTS)).opts(FIG_OPTS)
plot_field(field, line_source, phases)
Let's see what this looks like with several source points.
viz_data = {}
for n_points in [5, 10, 15, 20, 25, 30, 35, 40, 45]:
line_source, phases = make_symmetric_point_source(n_points, wavelength / 5.)
field = compute_amplitude_from_several_points(line_source, phases, X, Y, k)
viz_data[n_points] = plot_field(field, line_source, phases)
hv.HoloMap(viz_data, kdims='number of points')
viz_data = {}
for c_phi in np.linspace(2, 3, num=10):
wavelength_solid = c_phi / 1.5
line_source, phases = make_symmetric_point_source(80, wavelength / 10., wavelength_solid=wavelength_solid)
field = compute_amplitude_from_several_points(line_source, phases, X, Y, k)
viz_data[c_phi] = plot_field(field, line_source, phases)
hv.HoloMap(viz_data, kdims='phase velocity')
A way to show how this works is to make a polar diagram as a function of phase velocity.
We can use one of the features of holoviews
, a DynamicMap
to explore the parameters interactively:
def plot_with_options(c_phi, n_points, normalize=False):
wavelength_solid = c_phi / 1.5
line_source, phases = make_symmetric_point_source(n_points, wavelength / 10., wavelength_solid=wavelength_solid)
field = compute_amplitude_from_several_points(line_source, phases, X, Y, k)
if normalize:
field /= field.max()
return plot_field(field, line_source, phases)
dmap = hv.DynamicMap(plot_with_options, kdims=['c_phi', 'n_points']).redim.range(c_phi=(1.5, 6),
n_points=(0, 100)).redim.default(n_points=15)
# uncomment this for interactive exploration
hv.output(dmap, holomap='widgets')
The above function returns a plot for each call specifying the focal distance, angle and number of source points. Using repeated calls to the function as well as some interpolation between frames, we can try to remake the animation from the initial tweet with a little guesswork as to what the exact parameters are:
from tqdm import tqdm
states = [(10, 0), # c_phi, n_points for first state
20, # number of frames to interpolate
(10, 100), # c_phi, n_points for second state
20, # ...
(1.5, 100),
20,
(1.5, 2),
5,
(2.5, 2),
20,
(2.5, 100),
20,
(10., 100),
10,
(10, 0.)
]
def interp(start, end, nframes, frame, isint=False):
val = frame / (nframes - 1) * (end - start) + start
if isint:
val = int(val)
return val
hmap_data = []
frame_index = 0
for start_state, nframes, end_state in zip(states[::2], states[1::2], states[2::2]):
for frame in tqdm(range(nframes)):
c_phi = interp(start_state[0], end_state[0], nframes, frame)
n_points = interp(start_state[1], end_state[1], nframes, frame, isint=True)
hmap_data.append([frame_index, plot_with_options(c_phi, n_points, normalize=True)])
frame_index += 1
# optional: save animation locally
hv.save(hv.HoloMap(hmap_data), r'c:\Temp\holomap.gif', fps=5)