function [Mout, vout, vRF] = calc_magn(vin, RF1, N, th, E1, E2)
% Evolution of ssfp signal starting from vin.
% vin: input magnetization vs. offres.
% RF1 is 3 by 3 RF matrix.
% N number of RF pulses.
% th offres angles in radians.
% E1, E2: exp(-T1,2/TR).
% vout output magnetization vs offres.
% Mout M immediately after the RF pulse vs. RF pulse number and offres.
M = length(th);
Mout = zeros(M, N, 3);
vout = zeros(3, M);
Rel = [E2 0 0; 0 E2 0; 0 0 E1];
for jj = 1:M		% Theta loop.
    afree = exp(1j*th(jj)/2);
    free = Rel*Rot(afree, 0);
    v1 = vin(:, jj);
    for k = 1:N	% RF pulse loop.
        Mout(jj, k, :) = RF1*v1;
        v1 = free*RF1*v1 + [0 0 (1 - E1)].';
    end
    vout(:, jj) = v1;
end
vRF = squeeze(Mout(:, end, :));
