function [a0, am1, Mmax, ssfp] = calc_SSFP_signal(T1, T2, TR, alp_degrees)
alp = alp_degrees/180*pi;

E1 = exp(-TR/T1); E2 = exp(-TR/T2);
a = -(1 - E1)*E2*sin(alp);
b = (1 - E1)*sin(alp);
c = E2*(E1 - 1)*(1 + cos(alp));
d = 1 - E1*cos(alp) - (E1 - cos(alp))*E2^2;
MM = -c./d;
E = (1 - sqrt(1 - MM.^2))./MM;
a0 = a./c + (b./d - a./c)./sqrt(1 - MM.^2);
am1 = b./c + (a./d - b./c)./sqrt(1 - MM.^2);
Mmax = (b - a)./(d - c);
M = 301;
Phssfp = ((0:M - 1) - floor(M/2))*2*pi/M;
ssfp = (a*exp(-1j*Phssfp.') + b)./(c*cos(Phssfp.') + d);
plot_data = 1;
if plot_data > 0
    h = figure;
    set(h, 'Name', [' T1 ', num2str(T1), ', T2 ', num2str(T2), ', TR ', num2str(TR), ', flip deg ', num2str(alp_degrees), '. '])
    plot(Phssfp, real(ssfp), Phssfp, imag(ssfp), 'm'); grid;
    legend(' real ', ' imag ');
    xlabel('\phi rad');
    title(' ssfp signal ');
end
