function GaussPacket(t,params)

%  Function GaussPacket(t)
%  Plot Gaussian wave packet
%    params:
%     a
%     k0
%     xwidth
%     x
%     plotOmegaQuadraticFit

a = params.a;
k0 = params.k0;
xwidth = params.xwidth;
x = params.x;
plotOmegaQuadraticFit = params.plotOmegaQuadraticFit;

% Plot free surface at time t

xmin = min(x);
xmax = max(x);
xmid = 0.5*(xmin+xmax);
xp = x(1:(end-1));
np = length(xp);

L = xmax - xmin;
k = (2*pi/L)*[0:(np/2 - 1) 0 -(np/2-1):-1];
om = sqrt(abs(k)).*sign(k);

eta0 = a*exp(-(xp-xmid).^2/(2*xwidth^2)).*cos(k0*xp);
etahat0 = fft(eta0);
etahat = etahat0.*exp(-1i*om*t);
eta = ifft(etahat);

% Wave packet in physical space

subplot(2,1,1)
plot(x,eta([1:np 1]),'-')
ylabel('z')
zmax = 1.3*a;
zmin = -2.5*a;
text(0.98*xmin + 0.02*xmax, 0.95*zmin,sprintf('t = %5.2f',t));
axis([xmin xmax zmin zmax]);
title(['Gaussian deep water wave packet, k_0 = ' num2str(k0)])

if(plotOmegaQuadraticFit)
  om0 = k0^0.5;
  cg0 = 0.5*k0^-0.5;     % dom/dk(k0)
  beta0 = -0.25*k0^-1.5; % d2om/dk2(k0)
  omQF = (om0 + cg0*(abs(k)-k0) + 0.5*beta0*(abs(k)-k0).^2).*sign(k);
  etaQF = ifft(etahat0.*exp(-1i*omQF*t));

  subplot(2,1,2)
  plot(x,etaQF([1:np 1]),'-')
  ylabel('z')
  zmax = 1.3*a;
  zmin = -2.5*a;
  text(0.98*xmin + 0.02*xmax, 0.95*zmin,sprintf('t = %5.2f',t));
  axis([xmin xmax zmin zmax]);
  title('Using quadratic fit to dispersion reln around k_0')
end
    
