%## Autor: Peter Bokes  <peter.bokes zo stuba.sk>
%## podporené z projektu KEGA 006STU-4/2022
%## 
%## Skript pre demonštrovanie časového vývoja stavu častice v nekonečne hlbokej potenciálovej jame
%## z počiatočnej podmienky v tvare lineárne lomemenej funkcie. 
%## Vývoj je realizovaný rozvojom do konečného počtu N stacionárnych stavov.
%##
function Propagovanie_v_jame(L=2.0, a=0.5, T=4.0, Nt = 250, N=3)
%##
%## 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
%##
%
% šírka potenciálovej jamy 
% L = 2.0 
%
% Pozícia maxima počiatočnej podmienky
% a = 0.5
%
% diskretizácia osi x
Nx = 100
% 
% T = 4.0
% Nt = 500
%
% počet bázových funkcií 
% N = 50

n_vals = linspace(1, N, N );

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

for k = 1:Nt
   t = T/Nt*k
   c = coef(t,n_vals,a,L);

   % výpočet pravdepodobnosti
   for j = 1:Nx
      p_vals(j) = prob(x_vals(j),n_vals,c,L);
      p0_vals(j) = prob0(x_vals(j),a,L);
   end

   %grafy
   plot(x_vals,p_vals);
   hold on
   plot(x_vals,p0_vals);
   xlabel("x");
   ylabel("P(x)");
   grid on
   hold off

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

% koeficienty rozvoja
function c = coef(t,n,a,L)
   kn = pi.*n./L ;
   en = 0.5.*kn.^2;
   c = exp( -i.*en.*t)*(sqrt(6.0)/(a*L))*(1.0 + a/(L-a)).*sin(kn.*a)./(kn.*kn) ;
end

% vyskladaná pravdepodobnosť z bázových funkcií
function p = prob(x, n, c, L)
   kn = pi.*n./L  ;
   psi = sqrt(2./L).*sin(kn.*x) ;
   p = abs(dot(c,psi))^2 ;
end

% pravdepodobnosť priamo z počiatočnej podmienky
function p = prob0(x,a,L)
  if ( x < a ) 
     p = (3./a^2/L)*x^2;
  else
     p = (3./(L-a)^2/L)*(L-x)^2;
  endif 
end

end
