% an illustration of a plane (spherical) wave

plane_wave = 1;
spherical_wave = 2;

wave_type = plane_wave;
%wave_type = spherical_wave;

if wave_type == plane_wave

   % window size
   Lx=0.4;
   Ly=1;

   wavenumber = 60;
   
   % blow up the image by this factor to display better
   factor = 30;

   % a small shift to be added below for graph. purposes.
   shift = 1; 
   
elseif wave_type == spherical_wave
      
      Lx = 0.5;
      Ly = Lx;
      wavenumber = 100; 
      
      factor = 1;
      shift = 3;
end

Mx = Lx/2;
Wy = Ly/2;

M=400;
N = floor(M*Ly/Lx);

[X, Y]=meshgrid(linspace(-Lx/2, Lx/2, M), linspace(-Ly/2, Ly/2, N));

T = 1;
nt = 10;
Time = linspace(0, T, nt);

for repeat = 1:1

   % go over one time period of the field
   for iter = 1:(nt-1) % nt is same as 1 due to peridicity
      
      t = Time(iter);
      
      if wave_type == plane_wave
         
         % plane wave
         Z = real(exp(i*wavenumber*Y)*exp(-i*2*pi*t));
         
      elseif wave_type == spherical_wave
         
         % spherical wave
         Z = exp(sqrt(-1)*wavenumber*sqrt(X.^2+Y.^2))*exp(-sqrt(-1)*2*pi*t)./sqrt(X.^2+Y.^2);
      end
      
      % plot the real part of the field Z
      
      figure(1); clf; hold on; axis equal; axis off;
      image(factor*(real(Z+shift))); % add shift to Z for graphing purposes
      colormap jet; shading interp;

      file = sprintf('Movie_frame%d.eps', 1000+iter);
      disp(file); %show the frame number we are at
      saveas(gcf, file, 'psc2') %save the current frame
      
      pause(0.1);
      
   end

end

% The following command was used to create the animated figure.
% convert -antialias -loop 10000  -delay 15 -compress LZW Movie_frame10* Spherical_wave.gif