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

% Copyright W.Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 05/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;              % Grenze 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    
k=200;%50;%10;%15;                    % Wellenzahl und damit Geschwindigkeit des Pakets
%% Potential
%V = 150*ones(size(x));
V = 60000*ones(size(x));
V(x<xp1) = 0;
%V = cumsum(V);
V(x>xp2) = 0;
%V = zeros(size(x));
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))
%% 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 ausgef?hrt bis die Welle an die Raumgrenzen st??t
%  F?r die Visualisierung wird das Potential geeignet umskaliert
renn = true;
xort = xg;
%hpl=plot(x,abs(psi),xort,max(abs(psi))*1.015,'p'); %ylim([0,1.1])
nplot=0;
fh1 = figure;
fh2 = figure('Position',[520,363,1010,435]);
nsub=0;
while renn
%tic
hold off
%tic
psi = propa*psi;
%toc
%tic
figure(fh1)
area(x,abs(V)/max(abs(V)),'FaceColor',[0.96 0.96 0.96]),hold on
%toc
%[wert,pos]=max(abs(psi))
t=t+dt;
nplot=nplot+1;
xort = xg + 2*k*t;       % vg = 2*k Gruppengeschwindigkeit
                         % zum Vergleich Bewegung eines klassischen
                         % Teilchens
%tic
abspsi=abs(psi);
%skal=max(abspsi);
plot(x,abspsi/skal,xort,1.015,'p'), ylim([0,1.025])
%set(hpl(1),'yData',abs(psi)), set(hpl(2),'xData',xort,'yData',max(abs(psi))*1.015)
title(['t = ' num2str(t),'  ',num2str(nplot)]),shg
%toc
renn = (max(abs([psi(1),psi(end)])) < 0.001) && t<0.05 ;
drawnow;
if nplot==1
    nsub=nsub+1;
    figure(fh2)
    subplot(2,2,nsub)
    area(x,abs(V)/max(abs(V)),'FaceColor',[0.96 0.96 0.96]),hold on
    plot(x,abspsi/skal,'k',xort,1.015,'pk','MarkerSize',8,'MarkerFaceColor',[0.1,0.1,0.1]), ylim([0,1.025])
    title(['t = ' num2str(t)]),shg
elseif nplot==90
    nsub=nsub+1;
    figure(fh2)
    subplot(2,2,nsub)
    area(x,abs(V)/max(abs(V)),'FaceColor',[0.96 0.96 0.96]),hold on
    plot(x,abspsi/skal,'k',xort,1.015,'pk','MarkerSize',8,'MarkerFaceColor',[0.1,0.1,0.1]), ylim([0,1.025])
    title(['t = ' num2str(t)]),shg
    
elseif nplot==275
    nsub=nsub+1;
    figure(fh2)
    subplot(2,2,nsub)
    area(x,abs(V)/max(abs(V)),'FaceColor',[0.96 0.96 0.96]),hold on
    plot(x,abspsi/skal,'k',xort,1.015,'pk','MarkerSize',8,'MarkerFaceColor',[0.1,0.1,0.1]), ylim([0,1.025])
    title(['t = ' num2str(t)]),shg
elseif nplot==550
    nsub=nsub+1;
    figure(fh2)
    subplot(2,2,nsub)
    area(x,abs(V)/max(abs(V)),'FaceColor',[0.96 0.96 0.96]),hold on
    plot(x,abspsi/skal,'k',xort,1.015,'pk','MarkerSize',8,'MarkerFaceColor',[0.1,0.1,0.1]), ylim([0,1.025])
    title(['t = ' num2str(t)]),shg
end
%toc
end
toc