function [Mout, Mecho, M1out] = calc_M(vin, flip, phase_init, phase_inc, 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.
% Mecho M TR/2 after the RF pulse, i.e. half way between the current and
% the next RF pulses.
% M1out M before the next RF pulse vs. RF pulse number and offres.
% ======================================================================
M = length(th);
Mout = zeros(M, N, 3);
Mecho = Mout;
M1out = Mout;
Rel = [E2 0 0; 0 E2 0; 0 0 E1];
aRF = cosd(flip/2);
RF_phase = (phase_init + (0:N - 1)*phase_inc)*pi/180;
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.
        bRF = 1j*exp(1j*RF_phase(k))*sind(flip/2);
        RF1 = Rot(aRF, bRF);
        Mout(jj, k, :) = RF1*v1;
        Mecho(jj, k, :) = sqrt(free)*RF1*v1 + [0 0 (1 - sqrt(E1))].';
        v1 = free*RF1*v1 + [0 0 (1 - E1)].';
        M1out(jj, k, :) = v1;
    end
end
