% ================================================================
% John Kerl
% kerl at math dot arizona dot edu
% 2008-02-06
%
% Ideas here are taken from a talk given by Itai Seggev at the
% AMS/MAA joint meetings in San Diego in January 2008.
% ================================================================
% This software is released under the terms of the GNU GPL.
% Please see LICENSE.txt in the same directory as this file.
% ================================================================

% Closer masses give slower oscillations.
% Eigenvalue degeneracy when a == b ...
% theta = pi/4 promotes mixing; theta = 0 does not.

% I am taking hbar = 1 and E_0 = 1.

% ----------------------------------------------------------------
a     = 0.4;
b     = 0.5;
theta = 1.0 * pi/4;

dt    = 0.1;
niter = 100;

initv = 0;
pev   = 0;
pa    = 0;

% ----------------------------------------------------------------
H = [
	a^2*cos(theta)^2 + b^2*sin(theta)^2, (a^2-b^2)*cos(theta)*sin(theta);
	(a^2-b^2)*cos(theta)*sin(theta),      a^2*sin(theta)^2 + b^2*cos(theta)^2
	];

ket0 = [1.0; 0.0];
ket1 = [0.0; 1.0];

[V,D] = eig(H);

Udt = expm(-i * H * dt);

%psi0 = V(:,1); % Start with 1st eigenvector
%psi0 = V(:,2); % Start with 2nd eigenvector
psi0 = ket0;   % Start with 1st standard basis vector
%psi0 = ket1;   % Start with 2nd standard basis vector

% ----------------------------------------------------------------
disp('H:');                          disp(H);
disp('U(dt):');                      disp(Udt);
disp('Eigenvalues of H:');           disp(D);
disp('Eigenvectors of H(columns):'); disp(V);
disp('psi(0):');                     disp(psi0);
disp('Is H Hermitian?  Here are H^* and H:');
disp(H');
disp(H);
disp('Is U(dt) unitary?  Here are U(dt)^*, U(dt), and U^*(dt) U(dt):');
disp(Udt');
disp(Udt);
disp(Udt' * Udt);


psit = psi0;
ts     = [];
prob0s = [];
prob1s = [];
phz0s  = [];
phz1s  = [];
for iter = 1:niter
	t     = iter*dt;
	psit  = Udt * psit;
	prob  = conj(psit).*psit;
	phz   = phase(psit);

	%string = sprintf('iter %4d %11.7f %11.7f', iter, prob(1), prob(2));
	%disp(string);

	%disp(norm(psit));

	ts     = [ts,    t];
	prob0s = [prob0s, prob(1)];
	prob1s = [prob1s, prob(2)];
	phz0s  = [phz0s, phz(1)];
	phz1s  = [phz1s, phz(2)];

end

hold on;
plot(ts, prob0s);
plot(ts, prob1s);
plot(ts, phz0s);
plot(ts, phz1s);