Why is my IFFT plot showing different values from expected.

I am plotting total Electric field(incident and reflected) at some distance "zs" from the interface after it gets reflected from the dielectric interface . I am sweeping my frequency from 400 MHz to 2 GHz with stepping of 100 MHz. The time value calculated from IFFT should match with the formula "t = (2*zs)/speed of light(c)". But the results obtained from IFFT plot don't match.

close all;
clc;
clear all;
df=1e8;
f_low=0.4e9;
f_high=2e9;
freq = f_low:df:f_high;
Fs =6e9;
Nf =(Fs/df)+1;
epsilon_medium_real = 4;% for sand
epsilon_0 = 8.8544187817e-12; % permitivity of free space
epsilon_air_relv = ones*epsilon_0;
epsilon_medium_real_relv = epsilon_medium_real;
sigma =1e-5; % conductivity of dry sand
mu0 = pi*4e-7;
E0 = 1;
zs =60e-2; % source distance
zm = -5e-2;
w = 2*pi*freq;
epsilon_imag = sigma./(w*epsilon_0);
epsilon_medium= epsilon_0*((epsilon_medium_real_relv) -(1j*epsilon_imag));
beta_air =  w.*(sqrt(mu0*epsilon_air_relv)); % propagation constant in air
beta_medium = w.*(sqrt(mu0*epsilon_medium)); % propagation constant in medium
R = (sqrt(epsilon_air_relv)-sqrt(epsilon_medium))./(sqrt(epsilon_air_relv)+sqrt(epsilon_medium)); % reflection coefficient
T = (2*sqrt(epsilon_air_relv))./(sqrt(epsilon_air_relv)+sqrt(epsilon_medium));% transmission coefficient
Ey_air  = E0*((exp(-1j*beta_air*zs))+(R.*exp(1j*beta_air*zs))); % total field in air/ (Ey_i + Ey_r)
Ey_medium = E0*T.*exp(-1j*beta_medium*zm); % transmitted field in medium
Nf_half = (Nf-1)/2 + 1;
Evec_half = zeros(1,Nf_half);
fvec_half = linspace(0,Fs/2,Nf_half);
n1 = round(f_low/df);
n2 = round(f_high/df);
Evec_half(1,n1:n2) = Ey_air;
Evec(1,:) = [fliplr(conj(Evec_half)), Evec_half(1,2:Nf_half)];
fvec = linspace(-Fs/2,Fs/2,Nf);
subplot(211);
stem(freq,(Ey_air));
title('magnitude plot of electric field in Air');
xlabel('frequency');
ylabel('magnitude');
subplot(212);
stem(freq,imag(Ey_air));
title('phase plot of electric field in Air');
xlabel('frequency');
ylabel('phase');
figure(2);
subplot(211); stem(fvec,real(Evec));
subplot(212); stem(fvec,imag(Evec));
figure(3);
subplot(211); stem(fvec,abs(Evec));
subplot(212); stem(fvec,angle(Evec));
et = ifft(Evec,'symmetric');
dt = 1/Fs;
tmax = Nf*dt;
tvec = linspace(0,tmax,Nf);
figure(4); plot(tvec,real(et));