I have attempted to code up the simple math required to plot the total field set up by an incident plane wave on a semi-infinite flat plate which can be found here.
To summarise:
$$\phi_s(r,\theta ) = e^{ikr}/2 \left( w\left[ e^{i\pi/4}\sqrt{2kr}\sin(\frac{\theta + \Theta}{2}) \right]+ w\left[ e^{i\pi/4}\sqrt{2kr}\sin(\frac{\theta - \Theta}{2}) \right]\right) \enspace,$$
s the scattered potential in polar coordinates according to the article. $\Theta$ denotes the angle of the incoming plane wave and $k$ is the wave number. Finally
$$w(z) = e^{-z^2}(1-\mathrm{erf}(z)) \enspace$$
If you run the code you will see that I do not get the expected wave field as in the pictures.
My question is as follows: Have I made a mistake in my implementation or is the formula given incorrect?
My code is as follows:
import numpy as np
from scipy.special import erf
import matplotlib.pyplot as plt
"""Plane wave scattering by a plane wave incident on a seminfinite plane
"""
global T
global k
T = np.pi/4.0 # Incident plane wave angle
k = 10 # Wave number
def phi_scattered(r,t):
W = lambda z: np.exp(-z**2)*(1-erf(z))
term1 = W(np.exp(1j*np.pi/4.0)*np.sqrt(2*k*r)*np.sin((t+T)/2.0))
term2 = W(np.exp(1j*np.pi/4.0)*np.sqrt(2*k*r)*np.sin((t-T)/2.0))
res = np.exp(1j*k*r)/2.0*(term1 + term2)
return res
def phi_incident(r,t):
return np.exp(1j*k*r*np.cos(t-T))
x1 = np.linspace(-50,50,100)
x2 = np.linspace(-50,50,100)
X1, X2 = np.meshgrid(x1,x2)
R = np.sqrt(X1**2+X2**2)
THETA = np.arctan2(X2,X1)
phi = phi_scattered(R,THETA) + phi_incident(R,THETA)
plt.figure()
plt.imshow(np.real(phi), vmin = np.min(np.real(phi)), vmax = np.max(np.real(phi)))
plt.colorbar()
plt.show()