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]:

251



# Analytical expression

\begin{eqnarray} u_s(\rho, \phi, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \{ -i \pi H_0^{(2)}\left(k | \sigma - \sigma_s| \right) F(\omega) e^{i\omega t} d\omega\} \end{eqnarray}


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)




[  1.27868651e-03 +0.00000000e+00j   1.52464234e-01 -6.71253384e-02j
4.39030278e-01 -4.84469780e-01j   4.48789990e-01 -1.36552460e+00j
-2.49380371e-01 -2.46401395e+00j  -1.84470941e+00 -3.22829388e+00j
-4.11023678e+00 -3.01865400e+00j  -6.40188916e+00 -1.40008044e+00j
-7.84564286e+00 +1.61329225e+00j  -7.64813035e+00 +5.47024689e+00j
-5.41263793e+00 +9.19856699e+00j  -1.33671129e+00 +1.16922068e+01j
3.78890009e+00 +1.20693195e+01j   8.78591523e+00 +9.97319629e+00j
1.24296346e+01 +5.70956933e+00j   1.38091674e+01 +1.74012765e-01j
1.25875267e+01 -5.40235656e+00j   9.08191674e+00 -9.80037511e+00j
4.14907145e+00 -1.21396572e+01j  -1.07514044e+00 -1.20864139e+01j
-5.49093105e+00 -9.89084538e+00j  -8.30646899e+00 -6.26131622e+00j
-9.19359380e+00 -2.13155265e+00j  -8.29952847e+00 +1.59848746e+00j
-6.13362815e+00 +4.27193999e+00j  -3.38147513e+00 +5.58567564e+00j
-7.10048633e-01 +5.59063961e+00j   1.38221801e+00 +4.60388129e+00j
2.64054164e+00 +3.07446184e+00j   3.05156118e+00 +1.44883970e+00j
2.78313540e+00 +7.07116368e-02j   2.09649332e+00 -8.67816767e-01j
1.26058279e+00 -1.32514491e+00j   4.89671457e-01 -1.37357916e+00j
-8.69827265e-02 -1.14694011e+00j  -4.25696593e-01 -7.90362580e-01j
-5.48576276e-01 -4.24093980e-01j  -5.15073977e-01 -1.25614559e-01j
-3.94915795e-01 +7.18186540e-02j  -2.48596748e-01 +1.69706589e-01j
-1.17463358e-01 +1.90279184e-01j  -2.21846914e-02 +1.63091022e-01j
3.32913594e-02 +1.15296943e-01j   5.54751987e-02 +6.65226180e-02j
5.52609009e-02 +2.76974306e-02j   4.34585547e-02 +2.46035793e-03j
2.83406095e-02 -1.03244553e-02j   1.49158266e-02 -1.40580999e-02j
5.32681441e-03 -1.25238507e-02j  -2.49139608e-04 -8.80920626e-03j
-2.65961877e-03 -4.92129829e-03j  -3.05234072e-03 -1.85721745e-03j
-2.44272237e-03 +1.12585706e-04j  -1.53670004e-03 +1.12662139e-03j
-7.23548047e-04 +1.47890587e-03j  -1.52153641e-04 +1.46132893e-03j
1.73052920e-04 +1.29184928e-03j   3.11983619e-04 +1.10104361e-03j
3.36904638e-04 +9.48673998e-04j   3.06675174e-04 +8.49276141e-04j
2.59296015e-04 +7.94908280e-04j   2.14314504e-04 +7.70600609e-04j
1.78737339e-04 +7.62565220e-04j   1.52789241e-04 +7.61180433e-04j
1.34031986e-04 +7.61035265e-04j   1.19692036e-04 +7.59789564e-04j
1.07638653e-04 +7.56888330e-04j   9.65687632e-05 +7.52581493e-04j
8.58403777e-05 +7.47336901e-04j   7.52237057e-05 +7.41566874e-04j
6.46943126e-05 +7.35547625e-04j   5.43016146e-05 +7.29430940e-04j
4.41024602e-05 +7.23284687e-04j   3.41368092e-05 +7.17131179e-04j
2.44249007e-05 +7.10973304e-04j   1.49723857e-05 +7.04808696e-04j
5.77647624e-06 +6.98635685e-04j  -3.16948413e-06 +6.92454810e-04j
-1.18733683e-05 +6.86268463e-04j  -2.03430386e-05 +6.80080096e-04j
-2.85859383e-05 +6.73893542e-04j  -3.66090569e-05 +6.67712605e-04j
-4.44190135e-05 +6.61540859e-04j  -5.20221479e-05 +6.55381587e-04j
-5.94245823e-05 +6.49237778e-04j  -6.66322555e-05 +6.43112148e-04j
-7.36509390e-05 +6.37007171e-04j  -8.04862422e-05 +6.30925095e-04j
-8.71436150e-05 +6.24867966e-04j  -9.36283490e-05 +6.18837644e-04j
-9.99455797e-05 +6.12835811e-04j  -1.06100288e-04 +6.06863990e-04j
-1.12097304e-04 +6.00923551e-04j  -1.17941309e-04 +5.95015723e-04j
-1.23636838e-04 +5.89141606e-04j  -1.29188287e-04 +5.83302176e-04j
-1.34599912e-04 +5.77498298e-04j  -1.39875834e-04 +5.71730728e-04j
-1.45020043e-04 +5.66000129e-04j  -1.50036404e-04 +5.60307070e-04j
-1.54928656e-04 +5.54652037e-04j  -1.59700416e-04 +5.49035440e-04j
-1.64355188e-04 +5.43457613e-04j  -1.68896361e-04 +5.37918827e-04j
-1.73327212e-04 +5.32419289e-04j  -1.77650917e-04 +5.26959151e-04j
-1.81870543e-04 +5.21538510e-04j  -1.85989061e-04 +5.16157417e-04j
-1.90009345e-04 +5.10815875e-04j  -1.93934174e-04 +5.05513848e-04j
-1.97766238e-04 +5.00251263e-04j  -2.01508138e-04 +4.95028010e-04j
-2.05162394e-04 +4.89843948e-04j  -2.08731441e-04 +4.84698907e-04j
-2.12217635e-04 +4.79592689e-04j  -2.15623258e-04 +4.74525073e-04j
-2.18950517e-04 +4.69495815e-04j  -2.22201547e-04 +4.64504650e-04j
-2.25378417e-04 +4.59551295e-04j  -2.28483126e-04 +4.54635451e-04j
-2.31517611e-04 +4.49756804e-04j  -2.34483747e-04 +4.44915024e-04j
-2.37383348e-04 +4.40109771e-04j  -2.40218172e-04 +4.35340693e-04j
-2.42989920e-04 +4.30607430e-04j  -2.45700239e-04 +4.25909611e-04j
-2.48350725e-04 +4.21246858e-04j  -2.50942921e-04 +4.16618788e-04j
-2.53478326e-04 +4.12025010e-04j  -2.55958388e-04 +4.07465129e-04j
-2.58384512e-04 +4.02938747e-04j  -2.60758059e-04 +3.98445460e-04j
-2.63080347e-04 +3.93984863e-04j  -2.65352653e-04 +3.89556548e-04j
-2.67576216e-04 +3.85160106e-04j  -2.69752236e-04 +3.80795127e-04j
-2.71881877e-04 +3.76461198e-04j  -2.73966265e-04 +3.72157908e-04j
-2.76006496e-04 +3.67884845e-04j  -2.78003628e-04 +3.63641600e-04j
-2.79958691e-04 +3.59427760e-04j  -2.81872681e-04 +3.55242918e-04j
-2.83746566e-04 +3.51086666e-04j  -2.85581284e-04 +3.46958597e-04j
-2.87377746e-04 +3.42858307e-04j  -2.89136834e-04 +3.38785394e-04j
-2.90859405e-04 +3.34739459e-04j  -2.92546292e-04 +3.30720103e-04j
-2.94198301e-04 +3.26726933e-04j  -2.95816215e-04 +3.22759557e-04j
-2.97400796e-04 +3.18817584e-04j  -2.98952781e-04 +3.14900630e-04j
-3.00472887e-04 +3.11008310e-04j  -3.01961809e-04 +3.07140246e-04j
-3.03420225e-04 +3.03296061e-04j  -3.04848790e-04 +2.99475380e-04j
-3.06248141e-04 +2.95677836e-04j  -3.07618897e-04 +2.91903060e-04j
-3.08961661e-04 +2.88150690e-04j  -3.10277015e-04 +2.84420367e-04j
-3.11565528e-04 +2.80711734e-04j  -3.12827752e-04 +2.77024438e-04j
-3.14064221e-04 +2.73358132e-04j  -3.15275456e-04 +2.69712469e-04j
-3.16461965e-04 +2.66087107e-04j  -3.17624237e-04 +2.62481709e-04j
-3.18762752e-04 +2.58895940e-04j  -3.19877974e-04 +2.55329467e-04j
-3.20970355e-04 +2.51781964e-04j  -3.22040332e-04 +2.48253107e-04j
-3.23088334e-04 +2.44742573e-04j  -3.24114774e-04 +2.41250047e-04j
-3.25120057e-04 +2.37775215e-04j  -3.26104574e-04 +2.34317765e-04j
-3.27068707e-04 +2.30877390e-04j  -3.28012826e-04 +2.27453787e-04j
-3.28937291e-04 +2.24046655e-04j  -3.29842454e-04 +2.20655697e-04j
-3.30728654e-04 +2.17280619e-04j  -3.31596224e-04 +2.13921129e-04j
-3.32445484e-04 +2.10576941e-04j  -3.33276750e-04 +2.07247768e-04j
-3.34090325e-04 +2.03933330e-04j  -3.34886505e-04 +2.00633348e-04j
-3.35665578e-04 +1.97347546e-04j  -3.36427826e-04 +1.94075651e-04j
-3.37173518e-04 +1.90817393e-04j  -3.37902922e-04 +1.87572505e-04j
-3.38616293e-04 +1.84340723e-04j  -3.39313882e-04 +1.81121786e-04j
-3.39995932e-04 +1.77915433e-04j  -3.40662680e-04 +1.74721409e-04j
-3.41314355e-04 +1.71539459e-04j  -3.41951180e-04 +1.68369334e-04j
-3.42573374e-04 +1.65210783e-04j  -3.43181146e-04 +1.62063560e-04j
-3.43774701e-04 +1.58927422e-04j  -3.44354238e-04 +1.55802126e-04j
-3.44919951e-04 +1.52687434e-04j  -3.45472027e-04 +1.49583109e-04j
-3.46010649e-04 +1.46488914e-04j  -3.46535993e-04 +1.43404618e-04j
-3.47048231e-04 +1.40329990e-04j  -3.47547530e-04 +1.37264801e-04j
-3.48034052e-04 +1.34208824e-04j  -3.48507952e-04 +1.31161836e-04j
-3.48969385e-04 +1.28123612e-04j  -3.49418495e-04 +1.25093932e-04j
-3.49855428e-04 +1.22072577e-04j  -3.50280320e-04 +1.19059330e-04j
-3.50693307e-04 +1.16053974e-04j  -3.51094517e-04 +1.13056296e-04j
-3.51484076e-04 +1.10066084e-04j  -3.51862106e-04 +1.07083127e-04j
-3.52228724e-04 +1.04107215e-04j  -3.52584044e-04 +1.01138141e-04j
-3.52928175e-04 +9.81756981e-05j  -3.53261223e-04 +9.52196822e-05j
-3.53583289e-04 +9.22698894e-05j  -3.53894472e-04 +8.93261176e-05j
-3.54194868e-04 +8.63881658e-05j  -3.54484566e-04 +8.34558344e-05j
-3.54763654e-04 +8.05289248e-05j  -3.55032217e-04 +7.76072397e-05j
-3.55290336e-04 +7.46905830e-05j  -3.55538087e-04 +7.17787593e-05j
-3.55775545e-04 +6.88715746e-05j  -3.56002780e-04 +6.59688356e-05j
-3.56219861e-04 +6.30703503e-05j  -3.56426852e-04 +6.01759271e-05j
-3.56623814e-04 +5.72853757e-05j  -3.56810805e-04 +5.43985064e-05j
-3.56987881e-04 +5.15151304e-05j  -3.57155093e-04 +4.86350595e-05j
-3.57312491e-04 +4.57581062e-05j  -3.57460121e-04 +4.28840838e-05j
-3.57598027e-04 +4.00128061e-05j  -3.57726247e-04 +3.71440875e-05j
-3.57844821e-04 +3.42777430e-05j  -3.57953783e-04 +3.14135881e-05j
-3.58053164e-04 +2.85514387e-05j  -3.58142994e-04 +2.56911111e-05j
-3.58223297e-04 +2.28324220e-05j  -3.58294099e-04 +1.99751885e-05j
-3.58355419e-04 +1.71192280e-05j  -3.58407275e-04 +1.42643581e-05j
-3.58449682e-04 +1.14103968e-05j  -3.58482653e-04 +8.55716205e-06j
-3.58506196e-04 +5.70447209e-06j  -3.58520320e-04 +2.85214527e-06j]
[   0.    1.    2.    3.    4.    5.    6.    7.    8.    9.   10.   11.
12.   13.   14.   15.   16.   17.   18.   19.   20.   21.   22.   23.
24.   25.   26.   27.   28.   29.   30.   31.   32.   33.   34.   35.
36.   37.   38.   39.   40.   41.   42.   43.   44.   45.   46.   47.
48.   49.   50.   51.   52.   53.   54.   55.   56.   57.   58.   59.
60.   61.   62.   63.   64.   65.   66.   67.   68.   69.   70.   71.
72.   73.   74.   75.   76.   77.   78.   79.   80.   81.   82.   83.
84.   85.   86.   87.   88.   89.   90.   91.   92.   93.   94.   95.
96.   97.   98.   99.  100.  101.  102.  103.  104.  105.  106.  107.
108.  109.  110.  111.  112.  113.  114.  115.  116.  117.  118.  119.
120.  121.  122.  123.  124.  125.  126.  127.  128.  129.  130.  131.
132.  133.  134.  135.  136.  137.  138.  139.  140.  141.  142.  143.
144.  145.  146.  147.  148.  149.  150.  151.  152.  153.  154.  155.
156.  157.  158.  159.  160.  161.  162.  163.  164.  165.  166.  167.
168.  169.  170.  171.  172.  173.  174.  175.  176.  177.  178.  179.
180.  181.  182.  183.  184.  185.  186.  187.  188.  189.  190.  191.
192.  193.  194.  195.  196.  197.  198.  199.  200.  201.  202.  203.
204.  205.  206.  207.  208.  209.  210.  211.  212.  213.  214.  215.
216.  217.  218.  219.  220.  221.  222.  223.  224.  225.  226.  227.
228.  229.  230.  231.  232.  233.  234.  235.  236.  237.  238.  239.
240.  241.  242.  243.  244.  245.  246.  247.  248.  249.  250.]




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