function M = solvebloch(teta,M0,N,T1,T2,U,V)
% 
% Purpose : Solve the bloch equation
%
% Arguments in : 
% T   : pulse duration
% M0  : the initial magnetization
% N   : truncation order for the series expansion;
%       the size of the truncated matrix is 3*(2N+1) 
% T1,T2 : relaxation times
%
% Arguments out :
% M  : the magnetisation vector
%
% ARAMIS - (c) S. Balac  - ENSSAT - Université de Rennes 1 - 2007 (version 2.0)
%

global B0; 
global B1; 
global gamma; 
w0 = gamma*B0;
global eta0;

T = teta/(gamma*B1);

tau1 = 1/T1;
tau2 = 1/T2;

u1 = U(1);
u2 = U(2);
v1 = V(1);
v2 = V(2);

I = sqrt(-1);
b = [0 0 eta0*tau1]';  % RHS of Bloch equation

wa0 = (1/2)*gamma*(u2 + v1);
wb0 = (1/2)*gamma*(2*B1-v2 + u1);
A0=[-tau2,    0,   -wa0;
	0, -tau2,   wb0;
      wa0,  -wb0, -tau1];
wa2 = - (1/4)*gamma*(I*(u1 + v2) + v1 - u2);
wb2 =   (1/4)*gamma*(I*(u2-v1) + v2 + u1);
A2=[  0,    0, -wa2;
      0,    0,  wb2;
    wa2, -wb2,   0];

Am2=conj(A2);

% The truncated matrix A

A = [A0+2*N*I*w0*eye(3,3), Am2; A2, A0+2*(N-1)*I*w0*eye(3,3)];
for k=-N+2:N
  [p,q]=size(A);
  A=[A,[zeros(p-3,3);Am2]; [zeros(3,q-3),A2], A0-2*k*I*w0*eye(3,3)];
end

[P,D] = eig(A);

B=[zeros(1,3*N), b', zeros(1,3*N)]';
V0=[zeros(1,3*N), M0', zeros(1,3*N)]';

for k=1:size(D,1)
     lambda=D(k,k);
     e1(k)=exp(lambda*T);
     if abs(lambda)<eps 
       e2(k) = T;
     else
       e2(k)=(1/lambda)*(e1(k)-1);
     end
end

invP = inv(P);

V=P*diag(e1)*invP*V0 + P*diag(e2)*invP*B;

% computation of the magnetization at time T

m=[0;0;0];
for k=-N:N
  m=m+V(3*(N+k)+1:3*(N+k)+3)*exp(2*I*k*w0*T);
end

M=[cos(w0*T), sin(w0*T), 0; -sin(w0*T), cos(w0*T), 0; 0 , 0, 1]*m;
M=real(M);

