%% wpcayleyBos.m
% Beispiel: Zwei Bosonen mit bis auf das Vorzeichen gleichen k-Werten im
% Kasten; Potential gegen unendlich (aus numerischen Gruenden endlicher Wert 
% Stop-Button in der Abbildung rechts unten


% Vers. 1.2; 10.2013 
% Copyright W.Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
%% Parameter
tic
clear,clc,close all
n=2000;                     % Zahl finiten Differenzelemente
xl = -pi/20;                % linke raeumliche Grenze
xr = pi/20;                 % rechte raeumliche Grenze
xp1 = -0.025/2;             % Parameter dienen der Berechnung der Potentialbreite
xp2 = 2*0.075/2;
xg =0;% -0.1;               % Zentrum des Wellenpaktes
x=linspace(xl,xr,n);
dx = x(2)-x(1);             % raeumliche Aufloesung
dt = 5*dx^2;                % zeitliche Aufloesung s. Kommentare wpcayley.m
t=0;
lambda = 2 * dx^2/dt    
k=1000;%50;%10;%15;         % Wellenzahl und damit Geschwindigkeit des Pakets
%% Potential
V = zeros(size(x));
Vsc = 400000000;
V(3:7) = Vsc;               % Kastenwand
V(end-7:end-3) = Vsc;
%% Figure mit PushButton zum Beenden
fh = figure;
ph = uicontrol(fh, 'Style','PushButton','Units','normalized',...
    'Position',[0.93,0.0,0.08,0.1],'String','Stop',...
    'BackgroundColor', [0.901961 0.901961 0.901961],...
    'ToolTipString','Beendet Programm','CallBack','tend=0');
plot(x,V/Vsc)
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(-i*k.*x)).*exp(-a*(x-xg).^2)/2;
psi = psi.';
skal = max(abs(psi));
hpl=plot(x,abs(psi),xg,1.015,'pm');
%% 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 geeignet umskaliert
renn = true;
xort = xg;
tend = 0.001;            % Maximale Laufzeit
vg = 2*k;
ts = 0;
area(x,V/max(abs(V)),'FaceColor',[0.96 0.96 0.96]),hold on
while renn
%hold off
psi = propa*psi;
%area(x,V/max(abs(V)),'FaceColor',[0.96 0.96 0.96]),hold on
%drawnow
t=t+dt;
xort = xg + vg*(t-ts);       % vg = 2*k Gruppengeschwindigkeit
                         % zum Vergleich Bewegung eines klassischen
                         % Teilchens
if ((xort - x(end-8) > 0) && xort > 0)
    vg = -vg;
    xg = xort;
    ts = t;
elseif ((x(7) - xort > 0) && xort < 0) 
    vg = -vg;
    xg = xort;
    ts = t;
end

abspsi=abs(psi);
%plot(x,abspsi/skal,xort,1.015,'p'), 
hpl(1).XData = x;
hpl(1).YData = abspsi/skal;
hpl(2).XData = xort;
ylim([0,1.025])
title(['t = ' num2str(t)]),shg
drawnow
renn = (max(abs([psi(1),psi(end)])) < 0.001) && t<tend ;
%toc
end
toc