<p>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.</p>

<h1>The Greens functions (i.e., analytical solutions)</h1>

<p>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.</p>

<h2>The source type</h2>

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

<h2>The Greens function equation</h2>

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

<h2>The Greens functions</h2>

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

<h1>The MATLAB implentation</h1>

<p>I made the following changes w.r.t. your code:</p>

<ol>
<li>The frequency vector can be defined from <span class="math-container">$f_1\in(0,\frac{1}{\Delta t})$</span> rather than how you did it, which is <span class="math-container">$f_2\in(-\frac{1}{2\Delta t},\frac{1}{2\Delta t})$</span>. 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 <code>fftshift</code> operations.</li>
<li>I changed the Green's function to the 2-D case, i.e., I changed <code>const=-1i*pi</code> to <code>const=-1i/4</code>. </li>
<li>I changed the <code>nfft</code>, 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 <code>nfft</code>, to make the results relatively independent of the chosen <code>dt</code>.</li>
</ol>

<p>Thus, the final code, having made only small changes compared to yours, becomes:</p>

<pre><code>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 &lt;= f &lt;= 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')
</code></pre>

<p>Provided that the signal is sampled to the Nyquist frequency, this should give good results:
<a href="https://i.sstatic.net/ES4um.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ES4um.png" alt="Matlab figure results"></a></p>

CC BY-SA 4.0
2019-05-19T18:02:21.530Z
2019-05-19T18:02:21.530Z
Erik
<p>I have been trying to evaluate the analytical solution for a wave travelling in a <strong>homogeneous</strong>, <strong>infinite</strong> media. For a given source $S(t)$, the wave-field can be calculated at a distance $r$, for a given velocity $v$ </p>

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

<p>Here $F^{-1}$ represent inverse Fourier transform.</p>

<p>$S(\omega)$ is the frequency domain representation of $S(t)$</p>

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

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

<p>I followed this <a href="https://scicomp.stackexchange.com/questions/8288/analytic-solution-2d-scalar-wave-equation-in-cylindrical-coordinates-numerical-i"><strong>link</strong></a> 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&lt;f&lt;f_{max}/2$) but the code in link takes frequencies between $0&lt;f&lt;f_{max}$ for calculating Hankal function.</p>

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

<p><strong>Edit1</strong></p>

<p>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?</p>

<p><a href="https://i.sstatic.net/hOgJr.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/hOgJr.jpg" alt="enter image description here"></a></p>

<p><strong>Edit2</strong></p>

<p>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 <strong>dt</strong> (e.g. $dt \in 10^{-4}[1, 2, 3, 4,....,10]$) the shape and amplitude of the solution changes.
<a href="https://i.sstatic.net/BjvJ7.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/BjvJ7.jpg" alt="enter image description here"></a></p>

<p>Any help is appreciated.</p>

<pre><code>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 &lt; f &lt; +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')
</code></pre>

CC BY-SA 3.0
2017-01-16T02:34:58.860Z
have you tried physics or math stack exchange?
f.thorpe
CC BY-SA 3.0
2017-01-16T10:06:47.687Z
I have asked the question at computational science stack exchange but still waiting for a solution.
Amartya
CC BY-SA 3.0
2017-01-16T10:35:24.440Z
@Amartya Please [don't multi-post](http://scicomp.stackexchange.com/q/25952/3212).  If you have not received an answer after a week on CompSci and you believe you will have better luck here (which I would not assume), you can flag your question on CompSci and ask for migration.
gerrit
CC BY-SA 3.0
2017-01-18T22:04:53.823Z
Please don't discourage geophysics questions from this forum. I keep seeing comments like that every time someone asks a computational geophysics question.
Antonio
CC BY-SA 3.0
2017-02-07T09:08:34.423Z
Please supply a full example. We don't have the source file available so testing is impossible.
Jesper
CC BY-SA 3.0
2017-02-08T14:27:12.787Z
Thank you @WayoftheGeophysicist for mentioning the problem and organizing it in a better way. I have inserted the source part which was missing in the code.
Amartya
CC BY-SA 3.0
2017-02-08T15:37:45.033Z
Yo can ignore negative frecuencies, that's not an issue.
Santiago
CC BY-SA 3.0
2017-02-12T11:10:54.120Z
@Santiago after ignoring the negative frequencies, the pulse appears at 2.5 sec instead of  5 sec. It also does not have right type shape.
Amartya
<seismology><waves><seismic><matlab>
Analytical solution for seismic wavefield in MATLAB
Amartya
# Analytical solution for seismic wavefield in MATLAB

(Asked by: username_0 on 2017-01-15)

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**](https://scicomp.stackexchange.com/questions/8288/analytic-solution-2d-scalar-wave-equation-in-cylindrical-coordinates-numerical-i) 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?

[![enter image description here](https://i.sstatic.net/hOgJr.jpg)](https://i.sstatic.net/hOgJr.jpg)

**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.
[![enter image description here](https://i.sstatic.net/BjvJ7.jpg)](https://i.sstatic.net/BjvJ7.jpg)

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

```

**Tags:** seismology, waves, seismic, matlab

> **Comments:**
> **username_1**: have you tried physics or math stack exchange?  
> **username_0**: I have asked the question at computational science stack exchange but still waiting for a solution.  
> **username_2**: @username_0 Please [don't multi-post](http://scicomp.stackexchange.com/q/25952/3212).  If you have not received an answer after a week on CompSci and you believe you will have better luck here (which I would not assume), you can flag your question on CompSci and ask for migration.  
> **username_3**: Please don't discourage geophysics questions from this forum. I keep seeing comments like that every time someone asks a computational geophysics question.  
> **username_4**: Please supply a full example. We don't have the source file available so testing is impossible.  
> **username_0**: Thank you @WayoftheGeophysicist for mentioning the problem and organizing it in a better way. I have inserted the source part which was missing in the code.  
> **username_5**: Yo can ignore negative frecuencies, that's not an issue.  
> **username_0**: @username_5 after ignoring the negative frequencies, the pulse appears at 2.5 sec instead of  5 sec. It also does not have right type shape.  

---

## 1 answers

### Answer by username_6 on 2019-05-19. Score: 1 | Accepted answer

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:
[![Matlab figure results](https://i.sstatic.net/ES4um.png)](https://i.sstatic.net/ES4um.png)


---

