% ================================================================ % 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