% A more efficient way for multi shot SSFP signal combination.
% Echos are acquired at TE = TR/2.
% S0 - conj(Sm1) act as a single echo rotating at Th/2.
% S1 - conj(Sm2) behave as a single echo rotating at 3*Th/2.
% F0, the steady state signal at phi = pi, is lower than S0 - conj(Sm1) by 1/(1 + E).
% The noise of S0 - conj(Sm1) is sqrt(2) higher.
n = 1000;
N = 2;
Th = (0:n - 1)*2*pi/n + eps;	% Off resonance angle.
phi = (0:N - 1)*2*pi/N;
alp_degrees = 40;      % RF flip angle degrees.
alp = alp_degrees/180*pi;
RF_phase_degrees = 90;     % degrees.
RF_phase = RF_phase_degrees/180*pi;
T1 = 1000;
T2 = 100;
TR = 4;
TE = TR/2;
% TE = 0;
sig = 0.05;
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);
W = abs(am1/a0);                                    % weight of am1;
% W = 1;
F0 = (a*exp(-sqrt(-1)*pi) + b)./(c*cos(pi) + d);    % steady state signal at phi = pi.
Scom0 = a0 - am1;                                   % combined signal.
Swcom0 = a0 - W*am1;                                % weighted sum.
Mp = zeros(size(Th));
Mz = zeros(size(Th));
v = [0 sin(alp)*(1 - E1) cos(alp)*(1 - E1)].';
S = zeros(N, n);
Signal = zeros(N, n);
Sn = S;
for j = 1:N;
    % S(j, :) = sqrt(-1)*(a*exp(-sqrt(-1)*(Th + phi(j))) + b)./(c*cos(Th + phi(j)) + d).*exp(sqrt(-1)*Th*TE/TR);
    Signal(j, :) = exp(sqrt(-1)*(RF_phase + pi/2))*(a*exp(-sqrt(-1)*(Th + phi(j))) + b)./(c*cos(Th + phi(j)) + d).*exp(sqrt(-1)*Th*TE/TR);
%     S(j, :) = exp(-sqrt(-1)*(RF_phase + pi/2))*Signal(j, :);
    S(j, :) = (a*exp(-sqrt(-1)*(Th + phi(j))) + b)./(c*cos(Th + phi(j)) + d).*exp(sqrt(-1)*Th*TE/TR);
    Sn(j, :) = sig*randn(size(Th)) + sqrt(-1)*sig*randn(size(Th)) + S(j, :);
end;
% S1 = sum(S)/N;
F = fft(Sn)/N;
S0 = F(1, :);
Sm1 = F(end, :);
Scom = S0 - conj(Sm1);
figure; plot(Th, real(Scom.*exp(-sqrt(-1)*Th/2)), Th, imag(Scom.*exp(-sqrt(-1)*Th/2)), 'm'); grid;
xlabel(' Th ');
title([' Scom: snr = ', num2str(Swcom0/(sig*sqrt(1 + W^2)/sqrt(N)), 4), '. W = ', num2str(W, 4), '. snr of F0 = ', num2str(abs(F0)/(sig/sqrt(N)), 4), '. ']);
S1 = F(2, :);
Sm2 = F(end - 1, :);
Scom1 = S1 - conj(Sm2);
figure; plot(Th, real(Scom1.*exp(-sqrt(-1)*3*Th/2)), Th, imag(Scom1.*exp(-sqrt(-1)*3*Th/2)), 'm'); grid;
xlabel(' Th ');
title([' Scom1:  snr = ', num2str(E*Swcom0/(sig*sqrt(1 + W^2)/sqrt(N)), 4), '. ']);
disp(' ');
disp([' snr of S0 = ', num2str(a0/(sig/sqrt(N)), 4), '. ']);
disp([' Scom = ', num2str(Scom0, 4), '. Scom/(1 + E) = ', num2str(Scom0/(1 + E), 4), '. F0 = ', num2str(abs(F0), 4), '. ']);
