% demo waves in 2d, time domain, arbitrary dirichlet structures, PML
%
% barnett 11/11/03

nx = 70;
dt = 0.04;
ndamp = 4;

% build dof mask
dx = 1/nx;
xs = dx*(0:nx);
dof = zeros(nx+1, nx+1); % dofs for edges at x=0, x=1.
dmsk = [];
is = 1:nx+1;
n = 0;
for i = is
  x = dx*(i-1);
  for j = is
    y = dx*(j-1);
    if ~(abs(x-0.5)<0.02 & abs(y-0.6)>0.05); % single   
      %~(abs(x-0.3)<0.01 & (abs(y-0.5)>0.2 | abs(y-0.5)<0.1)) %
      %double slit
      n = n+1;
      dof(i,j) = n;
      id(n) = i;
      jd(n) = j;
      % damping mask (0 or 1), d is grid dist into PML
      dmx = max(0, abs(i-1-nx/2)-nx/2+ndamp);
      dmy = max(0, abs(j-1-nx/2)-nx/2+ndamp);
      dmsk(n) = (max(dmx, dmy)/ndamp)^2; %max(dmx, dmy)
    end
  end
end
if 0
  figure;
  surf(xs, xs, dof');
  figure;
  for i=is
    for j=is
      if dof(i,j)>0
	psi(i,j) = dmsk(dof(i,j));
      end
    end
  end
  surf(xs, xs, psi');
end

if 1
  disp('building H...');
  H = -4 * speye(n); % diag
  nnz = 0;
  ii = [];
  jj = [];
  for m=1:n
    i = id(m);
    j = jd(m);
    if (i<nx+1 & dof(i+1,j)>0)
      nnz = nnz+1;
      ii(nnz) = dof(i+1,j); % connect up
      jj(nnz) = m;
    end
    if (j<nx+1 & dof(i,j+1)>0)
      nnz = nnz+1;
      ii(nnz) = dof(i,j+1); % connect right
      jj(nnz) = m;
    end
  end
  nnz
  Hodu = sparse(ii, jj, ones(nnz,1), n, n); % off diag 'upper'
  H = H + Hodu + Hodu';
  H = H; % rescale
end
  
% initial condition
x = zeros(n, 1);
v = zeros(n, 1);

figure; 
sn = dof(nx*0.1,nx*0.3);
sn2 = dof(nx*0.1,nx*0.7);
t = 0;
for f=1:100  % display frame, or seconds
  psi = zeros(nx+1, nx+1);
  for i=is
    for j=is
      if dof(i,j)>0
	psi(i,j) = x(dof(i,j));
      end
    end
  end
  imagesc(xs, xs, psi'); set(gca, 'ydir', 'normal'); caxis([-1 1]);
  axis equal; title(sprintf('t = %f', t)); drawnow
  for s=1:(1/dt) % timesteps
    dv = dt*10*sin(pi*t/20)^2*(t<20)*cos(t);   %(t<1);
    v(sn) = v(sn) + dv;
    v(sn2) = v(sn2) + dv;
    %30*cos(10*t)*(t<10);
    x = x + dt*v;
    v = v.*(1 - dt*dmsk') + dt*H*x;  % including damping
    t = t + dt;
  end
end