% First set up the measuring device N = 16; L = 2*pi; m = 2.0; % Maximum reading of the dial dN = 1/N; theta = L*linspace(-0.5+dN,0.5,N)'; k = 2*pi*(1:N)'./L; % Define the fourier transform on the discrete circle F = exp(-1i*k*theta')/sqrt(N); % The operator which rotates by pi radians R = pi*F'*diag(k)*F; % Initial state of the apparatus is the dial pointing to 0. A0 = +(theta==0); % Plot probability clf plot(theta*m/pi,abs(A0).^2,'-sk'); shg pause plot(theta*m/pi,abs(expm(-1i*R)*A0).^2,'-sk'); pause % Example, measuring Sz for a spin-1/2 degree of freedom: Sx = 0.5*[0 1; 1 0]; Sy = 0.5 *[0 -1i; 1i 0]; Sz = 0.5*[1 0; 0 -1]; % Prepare the spin 1/2 in some random state psi0 = rand(2,1); psi0 = psi0/sqrt(psi0'*psi0); state = kron(A0,psi0); % Initial state of composite system: Hm = kron(R,Sz)/m; % Measurement hamiltonian U = @(t) expm(-1i*Hm*t); for t = 0:0.02:1; state_t = U(t)*state; prob_t = abs(state_t).^2; % Plot the spin up part of the state subplot(2,1,1); plot(theta*m/pi,prob_t(1:2:end),'-sk'); hold on; avgdialposition = theta'*prob_t(1:2:end)*m/(pi*sum(prob_t(1:2:end))); plot(avgdialposition,0,'-sr'); hold off; axis([-m m 0 1]); % Plot the spin down part of the state subplot(2,1,2); plot(theta*m/pi,prob_t(2:2:end),'-sk');hold on avgdialposition = theta'*prob_t(2:2:end)*m/(pi*sum(prob_t(2:2:end))); plot(avgdialposition,0,'-sr'); hold off; axis([-m m 0 1]); shg end