# Simulation Doppelspalt 
# mit Crank-Nicolson und duennbesetzten Matrizen (sparse)

l = 60;

s = 2;
wnd = [zeros(1,l/2-3*s/2),ones(1,s),zeros(1,s),ones(1,s),zeros(1,l/2-3*s/2)];
dom = zeros(l, l);
dom(2:l-1, 2:l-1) = ones(l-2, l-2);
for j=1:l
  dom(j,l-j+1) = wnd(j);
endfor

idx = dom;

pnt = zeros(sum(sum(dom)),2);

n = 0;
for j = 1:l
  for k = 1:l
    if(dom(j,k)>0)
      n++;
      idx(j,k) = n;
      pnt(n,:) = [j,k];
    endif
  endfor
endfor

lap = -4*speye(n,n);

nb = [0,1; 0,-1; 1,0; -1,0];
for j=1:n
  for k=1:4
    nn = pnt(j,:) + nb(k,:);
    if(dom(nn(1),nn(2)) > 0)
      lap(j,idx(nn(1),nn(2))) = 1;
    endif
  endfor
endfor

elo0 = zeros(l, l);
r = 3;
px = 4;
py = 4;
[x,y] = meshgrid(linspace(-r,r,l/2),linspace(-r,r,l/2));
elo0(1:l/2, 1:l/2) = exp(-i*(px*x+py*y)).*exp(-(x.^2 + y.^2));

u0 = zeros(n,1);
for j=1:n
  u0(j) = elo0(pnt(j,1),pnt(j,2));
endfor

ut = u0;
dt = 0.25;
for t=0:dt:20
  tmp = (speye(n) - i*lap*dt*0.5)*ut;
  ut = (speye(n) + i*lap*dt*0.5)\tmp;
  for j=1:n
    elo(pnt(j,1),pnt(j,2)) = ut(j);
  endfor
  surf(abs(elo));
  replot
endfor

