%
% generate a Gaussian wave packet in coordinate form
%
% arguments:
%   n --       dimension of matrix.  (should be a power of 2 because of FFT.)
%   rx,ry --   components of initial avg position.
%   kx,ky --   components of initial avg wave vector.
%   delta --   discrete space step.
%   eps --     half-width of Gaussian (in coordinate form).
%   ncoef --   normalization constant.
% Last six aren't long double because of compiler bugs.
%

function [pack] = generate_packet (n, rx, ry, kx, ky, delta, epp, ncoef)

nhalf = n/2;

pack = 0 * ([1:n]'*[1:n]);
 
% This one is not symmetric...
  for i = (-nhalf:(nhalf-1)),
    ii = rem(i+n,n)+1;
    for j = (-nhalf:(nhalf-1)),
% ncoef * exp(i*k0 \cdot r) * exp(-1/4*(r-r0)^2/\mu^2)
      ikr = delta*(i*kx + j*ky)*1i;
      r2 = -.25*((delta*i-rx).^2 + (delta*j-ry).^2)/(epp.^2);
      pack(ii,rem(j+n,n)+1) = exp(ikr) * exp(r2) * ncoef;
    end
  end

