% ================================================================
% 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 = 0.6 * pi/4;

dt    = 0.1;
niter = 1000;

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);

measurement_iter = round(4*niter/5);

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));

	% Measurement
	if (iter == measurement_iter)
		u = rand(1);
		if (u < prob(1))
			psit = [1;0];
		else
			psit = [0;1];
		end
	end

	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, 'k.', 'MarkerSize', 1);
plot(ts, prob1s, 'b.', 'MarkerSize', 1);
%plot(ts, phz0s, 'r');
%plot(ts, phz1s, 'g');

if 1
	set(gcf, 'PaperPositionMode', 'Manual');
	set(gcf, 'PaperPosition', [0 0 5.0 2.0])
	axis([0,niter*dt,-0.1,1.1]);
	print -depsc neutrinos
end