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