# A Phased Array Animation

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)

â€” Jacopo Bertolotti (@j_bertolotti) April 21, 2020

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

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)
```