Frolian's blog standing on the shoulders of PyGiants

• About
• Archive

A Phased Array Animation

Wed, 22 Apr 2020

Waves Physics Animations

This morning, I saw a cool looking tweet on sci-twitter. Since it's about waves and it looks very nice, I want to try to reproduce it and share the results.

That's the tweet:

#PhysicsFactlet (220)
Wavefront shaping: If you can control the relative phase of a number of point emitters, you can control the shape of your propagating wave.
(Shown: plane wave, focus, and Airy beam) pic.twitter.com/5QN5e1xZn2

— Jacopo Bertolotti (@j_bertolotti) April 21, 2020

The way I understand what's going on, we need the following building blocks to get to make an animation that looks like the one in the tweet:

• single cylindrical wave source on a grid
• line of cylindrical wave sources on a grid
• take into account phase delays to allow for focusing and beam deviation

Let's get started.

Cylindrical waves on a grid

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)

number of points per wavelenght: 20.0

Now let's write a function that allows us to compute the field amplitude.

def compute_amplitude_from_source_point(source_location, X, Y, k):
    x0, y0 = source_location
    R = np.sqrt((X - x0)**2 + (Y - y0)**2)
    amp = np.exp(1j * k * R) 
    amp[np.isinf(amp)] = 0.
    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=6)

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)

A line of cylindrical sources

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):
    pos = np.arange(-n_one_side, n_one_side + 1) * dx
    return np.c_[pos, np.zeros_like(pos)]

line_source = make_symmetric_point_source(5, wavelength / 2.)
line_source

array([[-2.5,  0. ],
       [-2. ,  0. ],
       [-1.5,  0. ],
       [-1. ,  0. ],
       [-0.5,  0. ],
       [ 0. ,  0. ],
       [ 0.5,  0. ],
       [ 1. ,  0. ],
       [ 1.5,  0. ],
       [ 2. ,  0. ],
       [ 2.5,  0. ]])

def compute_amplitude_from_several_points(source_locations, X, Y, k):
    field = np.zeros_like(X)
    for source_location in source_locations:
        field += compute_amplitude_from_source_point(source_location, X, Y, k)
    return field

field = compute_amplitude_from_several_points(line_source, X, Y, k)

def plot_field(field, line_source):
    return (plot(field).opts(FIG_OPTS) * hv.Points(line_source).opts(color='black', padding=0.05)).opts(FIG_OPTS)

plot_field(field, line_source)

Let's see what this looks like with several source points.

viz_data = {}
for n_points in [0, 1, 2, 3, 4, 5, 6, 7]:
    line_source = make_symmetric_point_source(n_points, wavelength / 2.)
    field = compute_amplitude_from_several_points(line_source, X, Y, k)
    viz_data[n_points] = plot_field(field, line_source)
    
hv.HoloMap(viz_data, kdims='number of points')