function bpm
close all

SOLVER=0; % wyznaczanie modow ?


lambda=.5;
k0=2*pi/lambda;

nm=1.2, % sr. wsp. zalamania
km=k0*nm; % do definicji amplitudy zespolonej
dx=1*lambda; % siatka
dz=2*lambda;
if SOLVER,
  dz=1i*dz; % propagacja w urojonym kierunku?
end


Nx=300; % rozmiar obszaru symulacji
Nz=1000;
A=zeros(Nx,Nz);

x=dx*((1:Nx) - Nx/2);
z=dz*(1:Nz);


n=ones(Nx,1); % wsp. zalamania
n(abs(x)<10)=1.5;

sgm=.5; % poczatkowy rozkald pola
A(:,1)=(exp(-(x-2).^2./sgm^2).*sign(x)).';
nrm=zeros(1,Nz);
r=dz / (2i * km*dx.^2);
ai=2+2*r+r*(dx^2).*(km^2-k0^2.*n.^2);

v=-r*ones(1,Nx-1);
M=diag(ai)+diag(v,1)+diag(v,-1); % Macierz ukladu trojdiagonalnego
Mi=inv(M);


for j=2:Nz, % propagacja w kierunku z
  if SOLVER,
    Nz(j-1)=norm(A(:,j-1));
    A(:,j-1)=A(:,j-1)./Nz(j-1);
   end;
  bi=((A(:,j-1).*(2-2*r-r*dx.^2.*(km.^2-k0^2.*n.^2)))+r.*([A(2:end,j-1);0]+[0;A(1:end-1,j-1)]));
  A(:,j)=Mi*bi;
end
figure(1)

imshow(abs(A));
%contour(real(A).^2);

hold on
contour(n);
hold off
%hold on
%plot(Nz);
%hold off
log(Nz(end))/k0/abs(dz)+nm

end

