
% simulation of ssfp signal evolution to steady-state.
% Applied with lc_ssfp.
% Simulates magnetization after a set of Nshift RF shifts in [-pi, pi].
% ======================================================================

% Initialize parameters.
% =================================
flip_ang = 5;       % RF flip angle deg.
T1 = 1000;
T2 = 1000;
% T1 = 10000;
% T2 = 8000;
TR = 6.7;     % in msec.
TR = 20;

% keep the catalyse option for later versions.
% ===============================================
catalyse = 0;     % catalyze to improve approach to steady state.

E1 = exp(-TR/T1); E2 = exp(-TR/T2);

% N = round(10*T1/TR);	    % Number of RF pulses.
% N = 5000;	    % Number of RF pulses.
N = 400;
M = 201;     % number of offres angles from -pi to pi.
th = (0:M - 1)*2*pi/(M - 1) - pi;       % phi in radians.
Nshift = 20;       % No of shifted ssfp sequences.
Nshift_vec = (0:Nshift - 1) - floor(Nshift/2);

% phase increment vector
% ==========================
phase_inc = Nshift_vec/Nshift*360;     % phase increment in degrees.

% RF phase parameters.
% ======================
RF_phase = 0;     % 1st RF phase, degrees.
phase_init = RF_phase;

% off resonance Phi in radians. Phi E [-pi, pi].
% ===================================================
Phi = -pi/2;
Phi = angle(exp(1j*Phi));       % limit to [-pi, pi].
[~, K] = min(abs(th - Phi));       % index of Phi in th.

% Theoretical steady-state Mx, My and Mz
% ==========================================
Th = (0:M - 1)*2*pi/(M - 1) - pi;
a = -(1 - E1)*E2*sind(flip_ang);
b = (1 - E1)*sind(flip_ang);
c = E2*(E1 - 1)*(1 + cosd(flip_ang));
d = 1 - E1*cosd(flip_ang) - (E1 - cosd(flip_ang))*E2^2;
% Th = circshift(Th, [0, -round(phase_inc(loop)/360*(M - 1))]);
ssfp = exp(1j*(RF_phase*pi/180 + pi/2))*(a*exp(1j*Th) + b)./(c*cos(Th) + d);
ssfp1 = exp(1j*(RF_phase*pi/180 + pi/2))*(a*exp(1j*phase_inc*pi/180) + b)./(c*cosd(phase_inc) + d);
ssfp_shift = th(K)/(2*pi)*Nshift;
ssfp1 = circshift(ssfp1, -round(ssfp_shift));       % theoretical ssfp with Nshift values shifted by th(K).

% Theoretical Mz.
% ====================
MxMy = -E2*sin(Th)./(1 - E2*cos(Th));
MyZeroPhase = (a*cos(Th) + b)./(c*cos(Th) + d);
% Zm = (-(sin(Th).*MxMy + cos(Th)).*imag(S)*sin(alp)*E2 + (1 - E1)*cos(alp))/(1 - E1*cos(alp));
Zm = (-(sin(Th).*MxMy + cos(Th)).*MyZeroPhase*sind(flip_ang)*E2 + (1 - E1)*cosd(flip_ang))/(1 - E1*cosd(flip_ang));
% h = figure;
% set(h, 'Name', [' T1 ', num2str(T1), ', T2 ', num2str(T2), ', TR ', num2str(TR), ', flip deg ', num2str(flip_ang), '. '])
% plot(th, real(ssfp), th, imag(ssfp), 'm', th, Zm, 'g'); grid;
% legend(' Mx ', ' My ', ' Mz ');
% title(' Theoretical steady state signal ');

% initialize variables before running the do loop
% ========================================================
v0 = [0 0 1];			% start from z magnetization.
vin = repmat(v0(:), 1, M);     % initial M = Mz vs. offres.
Mout = zeros(M, N, Nshift, 3);
Mdemod = Mout;
wbar = waitbar(0, [' Running calculation with ', num2str(N), ' RF pulses... ']);

for loop = 1:Nshift
    % Run the RF pulse train with N RF pulses.
    % ===========================================
    % phase_init = RF_phase;
    if loop > 1
        vin = squeeze(M1out(:, end, :)).';
    end
    % Mout is M immediately after the RF pulse, Mecho at the echo (TR/2
    % after the RF pulse), and M1out is TR after the RF pulse, just before
    % the next RF pulse.
    % Mout is M by N by Nshift by 3 matrix. It is the magnetization
    % immediately after the RF pulse.
    [Mout(:, :, loop, :), Mecho, M1out] = calc_M(vin, flip_ang, phase_init, phase_inc(loop), N, th, E1, E2);
    phase_end = mod(phase_init + (N - 1)*phase_inc(loop), 360);

    % catalyse. Do not use. Keep for later versions.
    % =====================================================
    if catalyse > 0
        NRF = 5;
        for jj = 1:40
            % for jj = 1:10
            vin = squeeze(M1out(:, end, :)).';
            phase_inc = phase_inc + 0.125;
            %         phase_inc = phase_inc + 1.5;
            phase_init = mod(phase_end + phase_inc, 360);
            [Mout1, ~, M1out] = calc_M(vin, flip_ang, phase_init, phase_inc, NRF, th, E1, E2);
            phase_end = mod(phase_init + (NRF - 1)*phase_inc, 360);
        end
        % end

        % Run second RF pulse train with N RF pulses.
        % ===============================================
        vin = squeeze(M1out(:, end, :)).';
        % phase_inc = phase_inc + 5;
        phase_init = mod(phase_end + phase_inc, 360);
        [Mout, Mecho, M1out] = calc_M(vin, flip_ang, phase_init, phase_inc, N, th, E1, E2);
    end

    % demodulate Mout.
    % ====================
    % RF_phase_vec = (phase_init + (0:N - 1)*phase_inc(loop))*pi/180;     % in radians.
    RF_phase_vec = (0:N - 1)*phase_inc(loop)*pi/180;     % in radians.
    Mdemod(:, :, loop, 1) = Mout(:, :, loop, 1).*exp(-1j*RF_phase_vec);
    waitbar(loop/Nshift);
end

close(wbar);
Mdemod(:, :, :, 2) = conj(Mdemod(:, :, :, 1));
Mdemod(:, :, :, 3) = Mout(:, :, :, 3);

% average Mdemod over the last 51 RF pulses.
dinput = mean(squeeze(Mdemod(K, end - 50:end, :, 1)));
[ssfp2, opt_angle] = call_fmin_angle(dinput, ssfp1);

h = figure;
set(h, 'Name', [' T1 = ', num2str(T1), ', T2 = ', num2str(T2), ', TR = ', ...
    num2str(TR), ', flip = ', num2str(flip_ang), ' deg, ', num2str(N), ' RF pulses.']);
% scale Figure height and width for dual figure;
fig_factor = 2.0;       % increase figure height by fig_factor.
fig_scale  = get(h, 'Position');
fig_scale(2) = fig_scale(2) - (fig_factor - 1)*fig_scale(4) + 10;
fig_scale(4) = fig_factor*fig_scale(4);
set(h, 'Position', fig_scale);

subplot(2, 1, 1);
plot(phase_inc*pi/180, real(ssfp1)), hold on, plot(phase_inc*pi/180, imag(ssfp1), 'm');
grid on;
title([' theoretical ssfp. \theta (off res shift) = ', num2str(-th(K)), ' rad. ']);
legend(' real ', ' imag ', 'Location', 'southeast')
ax = gca;
ax.XTick = pi*linspace(-1, 1, 5);
ax.XMinorTick = 'on';
ax.XTickLabel = {'-\pi', '-0.5\pi', '0', '0.5\pi', '\pi'};
xlabel(' Nshift (rad) ');
subplot(2, 1, 2);
plot(phase_inc*pi/180, real(ssfp2)), hold on, plot(phase_inc*pi/180, imag(ssfp2), 'm');
grid on;
title([' measured ssfp. opt angle = ', num2str(opt_angle, 3), ' deg. '])
legend(' real ', ' imag ', 'Location', 'southeast')
ax = gca;
ax.XTick = pi*linspace(-1, 1, 5);
ax.XMinorTick = 'on';
ax.XTickLabel = {'-\pi', '-0.5\pi', '0', '0.5\pi', '\pi'};
xlabel(' Nshift (rad) ');
