%## Autor: Peter Bokes  <peter.bokes zo stuba.sk>
%## podporené z projektu KEGA 006STU-4/2022
%## 
%## Skript pre demonštrovanie rozptylu úzkeho vlnového balíka na 
%## na potenciálovej bariére v tvare Diracovej delta-funkcie U(x) = alpha \delta(x)
%##
function Uzky_vlnovy_balik_na_bariere(L=100.0, k0=100, alpha = 200, dk=1, T=.4, Nt = 50)
%##
%## manuál pre octave: https://docs.octave.org/interpreter/index.html#SEC_Contents
%## rychly tutorial: https://sites.ualberta.ca/~kirchner/513/Matlab_Octave%20Tutorial.html
%##
% diskretizácia osi x
Nx = 400
% 
% T = 4.0
% Nt = 500

% pole x-hodnôt v ktorých počítame pravdepodobnosť
x_vals = linspace(-L/2, L/2, Nx);
% iniciovanie pravdepodobnosti, budú spočítané po jednotlivých hodnotách x
p_vals = linspace(-L/2, L/2, Nx);

for k = -Nt:Nt
   t = T/Nt*k;

   % výpočet pravdepodobnosti
   for j = 1:Nx
      p_vals(j) = prob(x_vals(j),t);
   end
 
   %grafy
   plot(x_vals,p_vals);
   hold on
   xlabel("x");
   ylabel("P(x,t)");
   axis([-Inf,Inf,0,0.2])  %% nastaví pevné hranice pre y-os
   grid on
   hold off

   pause(.01);
end
%##
%##  zadefunovanie funkcií
%## 

% vlnova funkcia uzkeho baliku na bariere
function p = Psi(x,t)
  eps = 1e-10;             %%  malé čílo na numerické odstránenie zdanivej singularity 
  %% koeficient transmisie a reflexie pri m = hbar = 1
  trans = i*k0/(alpha + i*k0);
  reflx = -1/(1 + i*k0/alpha);
  %%  3x použijeme na prednáške odvodený vzťah pre úzky vlnový balík, ale e^{ikx} je pri x<0 nahradené 
  %%  s e^{ikx} +  r(k0) e^{-ikx} a pri x>0 s t(k0) e^{ikx}
  if ( x < 0 ) 
     p = exp(i.*(k0.*x - 0.5.*(k0.^2).*t)).*2.*sin((x-k0.*t).*dk.*0.5)./(x .- k0.*t .+ eps)/(sqrt(2*pi*dk));
     p += reflx.*exp(i.*(-k0.*x - 0.5.*(k0.^2).*t)).*2.*sin((x .+ k0.*t).*dk.*0.5)./(x .+ k0.*t .+ eps)/(sqrt(2*pi*dk));
  else
     p = trans.*exp(i.*(k0.*x - 0.5.*(k0.^2).*t)).*2.*sin((x-k0.*t).*dk.*0.5)./(x .- k0.*t .+ eps)/(sqrt(2*pi*dk));
  endif 
end

% pravdepodobnost
function p = prob(x, t)
      p = abs(Psi(x,t)).^2;
end

end
