In [1]:
from scipy.fftpack import rfft
from scipy.special import hankel2
import numpy as np
import matplotlib.pyplot as plt
# Analytical solution for the constant medium acoustic wave equation
# Given for a nx by ny domain of velocity v on a h grid with source at position xsrc,ysrc , 2second propagation at dt=1ms
T = 1.0
dt=0.002 #time increment dt = .5 * hstep /maxv;
nt = int(T/dt)
nf = int(nt/2+1)
f0 = 15.0
fnyq = 1. / (2*(dt))
df = 1.0/T
faxis = df*np.arange(nf)
nf
Out[1]:
In [2]:
# Ricker wavelet time domain
def ricker(f, T,dt,t0):
t = np.linspace(-t0,T-t0, T/dt)
y = (1.0 - 2.0*(np.pi**2)*(f**2)*(t**2)) * np.exp(-(np.pi**2)*(f**2)*(t**2))
return y
rick=ricker(f0,T,dt,1.0/f0)
# Ricker wavelet in frew domain
R= np.fft.fft(rick)
R=R[0:nf-1]
nf=len(R)
nt=len(rick)
In [3]:
nx=200
ny=200
h=5 #space increment d = minv/(10*f0);
xsrc=100
ysrc=100
v=1500
U_a=np.zeros((nx,ny,nf),dtype=complex)
for a in range(1,nf-1):
k = 2*np.pi*faxis[a]/v
for m in range(0,nx):
for n in range(0,ny):
tmp = k*np.sqrt((h*(m - xsrc))**2 + (h*(n - ysrc))**2)
U_a[m,n,a] = -1j*np.pi* hankel2(0,tmp)*R[a]
In [4]:
# Do inverse fft on 0:dt:T (third dimension on U) and you have analytical solution o
U_t=np.zeros((nx,ny,nt))
for m in range(0,nx):
for n in range(0,ny):
U_t[m,n,:]=np.real(np.fft.ifft(U_a[m,n,:],nt))
In [5]:
from matplotlib import animation
fig = plt.figure()
plts = [] # get ready to populate this list the Line artists to be plotted
plt.hold("off")
for i in range(1,nf):
r = plt.imshow(U_t[:,:,i], vmin=-.1, vmax=.1) # this is how you'd plot a single line...
plts.append( [r] )
ani = animation.ArtistAnimation(fig, plts, interval=50, repeat = False) # run the animation
plt.show()