%
%  code06.m -- wave radiation (djm, 26 sept 01)
%

%  parameters
circ = 0:pi/40:2*pi;  r0 = 1.0;  IO = 1;

%  coordinates
L = 4*pi;  x = -L:L/20:L;  y = -L:L/20:L;
T = 8*pi;  t = 0:T/32:T;

%  2D coordinates (rectangular & polar)
[xg,yg] = meshgrid(x,y);
zg = xg + i*yg;   rg = abs(zg);   tg = angle(zg);
re = max(rg,r0/2);

%  colormap (cold,blue -> hot,red)
z16 = zeros(1,16);  o32 = 32*ones(1,32);
djm  = [[z16,(0:2:32),o32];[(0:1:32),(31:-1:0)];[o32,(32:-2:0),z16]]'/32;

%  2D plot: time sequence

figure(1);  clf;  colormap(djm)
for j=1:size(t,2)
  v = besselh(0,IO,re)./besselh(0,IO,r0);
  u = v*exp(-i*t(j));

  pcolor(x,y,real(u));  shading interp;  hold on
  caxis([-1 1]);  colorbar
  axis(L*[-1 1 -1 1]);  axis square

  fill(cos(circ),sin(circ),'k');  hold off;  
  title(['\bf wave radiation:  real part @ t = ' num2str(t(j))])
  pause(0.1)
end

%  radial plot: time sequence

xp = r0:L/50:4*L+r0;

figure(2);  clf;
for j=1:size(t,2)
  v = besselh(0,IO,xp)./besselh(0,IO,r0);
  u = v*exp(-i*t(j));

  plot(xp,real(u),'r');  hold on
  axis([r0 4*L -1 1]);

  title(['\bf wave radiation:  real part @ t = ' num2str(t(j))])
  xlabel('\bf r-axis');  ylabel('\bf real(u)')
  pause(0.5)
  plot(xp,real(u),'g');  hold on
end

plot(xp,real(u),'r');  hold on