%% Visualsierung eines Wellenpakets im Kastenpotential
% wpcayley.m
% Beispiel: Rechteckpotential

% Vers. 2.8; 08.2013 
% Copyright W.Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016
%% Parameter
tic
clear,clc,close all
n=2000;                  % Zahl finiten Differenzelemente
xl = -1;                 % linke raeumliche Grenze
xr = 1;                  % rechte raeumliche Grenze
xp1 = -0.025/2;          % Grenzen des Potentials
xp2 = 2*0.075/2;
%xp2 = xp1+2*pi
xg = -0.1;               % Zentrum des Wellenpaktes
x=linspace(xl,xr,n);
dx = x(2)-x(1);          % raeumliche Aufloesung dx=epsilon
dt = 2*dx^2;             % zeitliche Aufloesung delta korrekt 2, fuer kleinere k-Werte auch bis 10 konvergent
t=0;
lambda = 2 * dx^2/dt     % vergl. Anmerkungen zur Konvergenz
k=200;%50;%10;%15;                    % Wellenzahl und damit Geschwindigkeit des Pakets
%% Potential
%V = 150*ones(size(x));
V = 80000*ones(size(x));
V(x<xp1) = 0;
V(x>xp2) = 0;
plot(x,V)
hold on
%% Start Wellenpaket
a = 4*4/(xp2-xp1)^2;                 % Bestimmt die Breite des Wellenpakets
psi = sqrt(sqrt((2*a)/pi))*exp(i*k.*x).*exp(-a*(x-xg).^2);
psi = psi.';
plot(x,abs(psi))
skal = max(abs(psi));                % Skalierung auf 1 zur Visualisierung
%% Erzeugen der Matrixstrukturen
xt = ones(1,n-1);
xt0 = zeros(1,n-1);
xt1 = zeros(n,1);
xm = diag(xt);
xm = [xt0;xm];
xm = [xm,xt1];
alpha = (i*lambda -2 - dx^2 * V);    % vgl. Buch
A = diag(alpha)+xm+xm'; A = sparse(A);            
beta = (i*lambda + 2 + dx^2 * V);
B = diag(beta)-xm-xm';  B = sparse(B);
propa = A\B;
propa = full(propa);
%% Zeitintegration
%  wird ausgefuehrt bis die Welle an die Raumgrenzen stoesst
%  Fuer die Visualisierung wird das Potential auf "1" umskaliert
renn = true;
xort = xg;
while renn
hold off
psi = propa*psi;
area(x,abs(V)/max(abs(V)),'FaceColor',[0.96 0.96 0.96]),hold on
t=t+dt;
xort = xg + 2*k*t;       % vg = 2*k Gruppengeschwindigkeit
                         % zum Vergleich Bewegung eines klassischen
                         % Teilchens. Bei Konvergenz bewegt sich das
                         % Wellenpaket und das klassische Teilchen mit
                         % derselben Geschwindigkeit.
abspsi=abs(psi);
plot(x,abspsi/skal,xort,1.015,'p'), ylim([0,1.025])
title(['t = ' num2str(t)]),shg
renn = (max(abs([psi(1),psi(end)])) < 0.001) && t<0.05 ;
drawnow;
end