Analytical solution for seismic wavefield in MATLAB

Question

I have been trying to evaluate the analytical solution for a wave travelling in a homogeneous, infinite media. For a given source $S(t)$, the wave-field can be calculated at a distance $r$, for a given velocity $v$

$$W(r)=F^{-1} [-i \pi S(\omega)H_0^{(2)} (kr) ]$$

Here $F^{-1}$ represent inverse Fourier transform.

$S(\omega)$ is the frequency domain representation of $S(t)$

$H_0^{(2)} (kr)$ is the Hankel function of zero order and second kind.

and $k$ is the wave number $(=\frac{\omega}{v})$.

I followed this link which demonstrate code in Python, however, I couldn't implement same in MATLAB. I also had doubts how it works. Since the fourier transform have positive as well as negative frequencies (i.e. $-f_{max}/2<f<f_{max}/2$) but the code in link takes frequencies between $0<f<f_{max}$ for calculating Hankal function.

I wrote following MATLAB code but got the wave-field all entries as NaN..!

Edit1

As a little modification, I forced the first term of the Hankel function to be zero (since it was containing a NaN in the real part). Now it is working nicely but it produced two peaks. I thought it to be similar to frequency shift effect, as observed after FFT for positive and negative frequencies, but it is not. In this case amplitudes are reversed. Can anyone explain it?

Edit2

If I assume that the only the first peak in the solution is correct one and second arises because of some effect (which I don' know). Problem comes when I change the dt (e.g. $dt \in 10^{-4}[1, 2, 3, 4,....,10]$) the shape and amplitude of the solution changes.

Any help is appreciated.

clc;    clear all;      close all;

%% Parameter Setting
vel=2000;   % velocity of medium
r=1000;     % distance at which wavefield will be observed
% the final wavelet should appear at time, t=r/vel=0.5sec

T=1;        % Total time
dt=.0001;   % time step
f0=25;      % central freq

%% Create Source Wavelet
N = round(T/dt);        % No of samples
t0 = 5/(sqrt(2)*pi*f0);     % Zero time shift
t = dt*(0:N-1);             % time vector
tau = t-t0;     % time vector shifted by t0
src = (1 -2* tau.*tau * f0^2 * pi^2).*exp(-tau.^2 * pi^2 * f0^2);  %source     wavelet

%%FT of signal
nfft=2^nextpow2(N);     % no of fft points
fmax=1/dt;              % max frequency 
freq= fmax*(-nfft/2:nfft/2-1)/nfft;    % Freq vector:  -Fmax/2 < f < +Fmax/2 

sw=fft(src,nfft)/nfft;                % fourier transform of wavelet, S(w)
sw_shift=fftshift(sw);      % shifting to make S(w) amplitude corrosponding to above "freq" vector 
amp= abs(sw_shift);         % amplitude of the S(w)

%% Green Function 
%(I have taken only positive freq. since hankel func defined only for positive numbers.
% Suggestions please..May be I am wrong..!)                  
freq_new= fmax*(0:nfft/2-1)/nfft;
w=2*pi*freq_new;
const = -1i*pi;
H02 = besselh(0,2,w*r/vel);
H02(1)=eps;
GF = const*H02;

%% Calculating Wavefield 
%I have made Green Function (GF) symmetric and then multipled with source spectrum.
GF_new=[fliplr(GF),GF]; 
WF= ifft(GF_new.*sw);
t_new=dt*(0:length(WF)-1);

%% Plot
figure(); 

% Wavelet
subplot(2,2,1);     
plot(t,src);     title('\bf{Source wavelet}'); 
xlabel('Time(sec)');    ylabel('Amplitude ')

% Frequencies
subplot(2,2,2);   
plot(freq,amp);    title('\bf{Source freq spec}');  
xlabel('Frequencies(Hz)');  ylabel('Amplitude ')

% Gneens Fkt
subplot(2,2,3); 
plot(freq_new,abs(GF));   title('\bf{Green Function}'); 
xlabel('Frequencies(Hz)'); ylabel('Amplitude ')

% Wafefield
subplot(2,2,4);      
plot(t_new,real(WF));    title('\bf{Trace}');   
xlabel('Time(sec)');  ylabel('Amplitude')

Answer 1

So, I know this question was asked a long time ago. But in the spirit of Stack Exchange, I will post my answer for future users perhaps.

The Greens functions (i.e., analytical solutions)

The first thing to be aware of are the various solutions to the wave equations change for the (1) source type and (2) the dimension.

The source type

Assuming acoustic waves in a homogeneous domain, that means you are either solving: $$ \frac{\partial^2 p(x,t)}{\partial t^2} = c^2(x)\left(\nabla^2 p(x,t) +f(x,t) \right) , \tag{1}$$ (where $p$ is pressure, $c$ the compressional velocity and $f$ the source function), or you are solving: $$ \left\{ \begin{array}{rl} \frac{\partial p(x,t)}{\partial t} & = -c^2(x)\bigg( \nabla\cdot v(x,t) + g(x,t) \bigg), \\ \frac{\partial v(x,t)}{\partial t} & = -\nabla p(x,t). \end{array} \right. \tag{2} $$ The difference between (1) and (2) is a time differentiation: $$ f(x,t) = \frac{\partial}{\partial t}g(x,t). $$ This is an important difference in source functions!

The Greens function equation

It is standard to restate equation (1) in the following form: $$ \left( \frac{1}{c^2(x)}\frac{\partial^2}{\partial t^2} - \nabla^2 \right)p(x,t) = f(x,t). \tag{3a}$$ Transforming this equation to the frequency domain, for which $\partial^2/\partial t^2 \to -\omega^2$, and where $\omega^2/c^2=k^2$ we obtain: $$ \left( -k^2 - \nabla^2 \right)p(x,\omega) = f(x,\omega), \tag{3b}$$ or $$ \left( k^2 + \nabla^2 \right)p(x,\omega) = -f(x,\omega), \tag{3c}$$ We can consult, e.g., https://www3.nd.edu/~atassi/Teaching/AME%2060633/Notes/greens.pdf to find the Greens functions (for which $f(x,\omega)=\delta(x-x_s)$ with $x_s$ the source location).

The Greens functions

I'll just copy the results. For equations of type (1) we have the Green's functions $G$: $$ \begin{array}{lrl} \text{1-D:} & G & = −\frac{i}{2k}e^{−ik\lvert x−x_s \rvert}, \\ \text{2-D:} & G & = −\frac{i}{4}H_0^{(2)}(k \lvert x−x_s \rvert),\\ \text{3-D:} & G & = \frac{1}{4\pi}\frac{e^{−ik\lvert x−x_s \rvert} }{\lvert x−x_s \rvert}. \end{array} $$ Indeed, $H_0^{(2)}$ is the Hankel function, in MATLAB we simply write that with besselh(0,2, k*r).

The MATLAB implentation

I made the following changes w.r.t. your code:

1. The frequency vector can be defined from $f_1\in(0,\frac{1}{\Delta t})$ rather than how you did it, which is $f_2\in(-\frac{1}{2\Delta t},\frac{1}{2\Delta t})$. In this way, the Hankel function only receives positive frequencies. Moreover, it is exactly the order in which the FFT produces its output. I thus also removed the fftshift operations.
2. I changed the Green's function to the 2-D case, i.e., I changed const=-1i*pi to const=-1i/4.
3. I changed the nfft, to just using the FFT of the size of the wavelet or time vector. This made it easier to obtain results of similar amplitudes... Also, I don't divide the source function by nfft, to make the results relatively independent of the chosen dt.

Thus, the final code, having made only small changes compared to yours, becomes:

clc;    clear all;

%% Parameter Setting
vel=2000;   % velocity of medium
r=1000;     % distance at which wavefield will be observed
% the final wavelet should appear at time, t=r/vel=0.5sec

T=1;        % Total time
dt=.0025;   % time step
f0=25;      % central freq

%% Create Source Wavelet
N = round(T/dt);        % No of samples
t0 = 5/(sqrt(2)*pi*f0);     % Zero time shift
t = dt*(0:N-1);             % time vector
tau = t-t0+1*dt;     % time vector shifted by t0
src = (1 -2* tau.*tau * f0^2 * pi^2).*exp(-tau.^2 * pi^2 * f0^2);  % source wavelet

%%FT of signal
fmax=1/dt; % max frequency 
freq= linspace(0,fmax,N);    % Freq vector:   0 <= f <= Fmax
sw=fft(src);      % fourier transform of wavelet, S(w)

%% 2-D Greens Function 
w=2*pi*freq;
const = -1i/4;
H02 = besselh(0,2,w/vel*r);
H02(1)=0;
GF = const*H02;

%% Calculating the Wavefield  
WF= ifft(GF.*sw);

%% Plot
% Wavelet
subplot(2,2,1);     
plot(t,src);     title('\bf{Source wavelet}'); 
xlabel('Time(sec)');    ylabel('Amplitude ')

% Frequencies
subplot(2,2,2);   
plot(freq,abs(sw));    title('\bf{Source freq spec}');  
xlabel('Frequencies(Hz)');  ylabel('Amplitude ')
xlim([0 fmax/2])

% Gneens Fkt
subplot(2,2,3); 
plot(freq,abs(GF));   title('\bf{Green Function}'); 
xlabel('Frequencies(Hz)'); ylabel('Amplitude ')
xlim([0 fmax/2])

% Wafefield
subplot(2,2,4);      
hold on
plot(t,real(WF));    title('\bf{Trace}');   
xlabel('Time(sec)');  ylabel('Amplitude')

Provided that the signal is sampled to the Nyquist frequency, this should give good results: