%% Visualisierung der Korrelationsfunktion - Revivals

% Vers. 2.8; 08.2013 
% Copyright W.Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

clear,close all,clc
%%
L = 2;       % Breite des unendlichen Potentialtopfs
lauf = 50;   % Anzahl der Entwicklungskoeffzienten
xs = 0.25;   %0.75;%0;0.1;%0;%0.5; % Verschiebung des Gausspakets in x-Richtung
gw=6;        % Breite des Gausspakets im Beispiel

fhcoeff = @(n,x) sin(pi*n/L*(x+L/2)).*exp(-gw*(x-xs).^2); % Eigenfunktion * Gauss zur Berechnung 
                                                          % der Entwicklungskoeffizienen
                                                          % a(k)
 
% nur die ungeraden n tragen fuer xs = 0 bei. 
% Grund gerade haben genauso hohen positiven wie negativen Anteil links und rechts von x=0
for k = 1:lauf
    n = k;           %2*k-1; % fuer xs=0
    a(k) = integral(@(x)fhcoeff(n,x),-L/2,L/2);
end
%% Wellenpaket
phi = @(n,x) sin(pi*n/L*(x + L/2));  % Eigenfunktionen unnormiert

%% Visualisierung
xv = linspace(-L/2,L/2);              % raeumliche Aufloesung (100 Schritte)
gauss = exp(-gw*(xv-xs).^2);          % dient nur dem Vergleich
wpaket = zeros(size(xv));
for k = 1:lauf
    wpaket = wpaket + a(k)*phi(k,xv); % Berechnung des Wellenpakets zur Zeit t=0
end
plot(xv,gauss,xv,wpaket,'o'),shg          % Vergleich mit Potential-freiem Gausspaket 
%%
En = @(n) (n*pi/L).^2/2;              % Energieeigenwerte
%% Berechnung des Korrelationskoeffzienten
t = linspace(0,2*pi*1.01,1000);       % Berechnungszeitpunkte
wp_corr =  0;
for k = 1:lauf
    wp_corr = wp_corr+a(k).^2*exp(i*En(k)*t);
end
wp_corr=wp_corr/wp_corr(1);
figure, plot(t,abs(wp_corr),'k'),shg  % Visualisierung
hold on
plot([t(1),t(end)],[1,1],'k--')
axis tight
xlabel('Zeit'), ylabel('Korr')
title('Korrelationskoeffizient');
ylim([0.4,1.05])
hold off
%% Animation
figure
tm = linspace(0,2*5.094,100);
for mk = 1:length(tm)
    xvisu = zeros(size(xv));
    for k = 1:lauf
        xvisu = xvisu + a(k).*phi(k,xv).*exp(i*En(k)*tm(mk));
    end
    plot(xv,abs(xvisu)),ylim([0,1.3])% ylim([0,1])
    F(mk) = getframe;
end
movie(F,1,3)
        