%
%  spin 1/2 in static+rotating magnetic field
%
omega = 1.5;
B = 1;
alpha = 0.7;       % angle with Oz
B0 = B*cos(alpha);  % z-component
B1 = B*sin(alpha);  % xy-component

ep = sqrt( (B0-omega/2)^2 + B1^2 );
time = linspace(0,2*pi/ep,5000);

ap = B1;
bp = (B0-omega/2-ep);
np = 1/sqrt(ap^2+bp^2);
ap = ap*np;
bp = bp*np;

am = bp;
bm = -ap;

Ap = bm/(ap*bm-am*bp);
Am =-bp/(ap*bm-am*bp);

phi_1 = Ap*ap*exp(-i*( omega/2-ep)*time) + Am*am*exp(-i*( omega/2+ep)*time);
phi_2 = Ap*bp*exp(-i*(-omega/2-ep)*time) + Am*bm*exp(-i*(-omega/2+ep)*time);

sx = conj(phi_1).*phi_2 + phi_1.*conj(phi_2);
sy = -i*conj(phi_1).*phi_2 + i*phi_1.*conj(phi_2);
sz = abs(phi_1).^2 - abs(phi_2).^2;

energy = B1*(cos(omega*time).*sx+sin(omega*time).*sy)+B0*sz;
B = sqrt(B0^2+B1^2);

subplot(2,1,1)

plot(time,sx,'g',time,sy,'r',time,sz,'b')
legend('s_x','s_y','s_z',0)
xlabel('time')
ylabel('spin')
grid on

subplot(2,1,2)

plot(time,energy)
axis([0,max(time),-B,B])
xlabel('time')
ylabel('energy')
grid on

