% ssfp_Pauly_equil.m;
% Calculates steady state magnetization using Cayley-Klein parameters from John Pauly paper.
% Apply N RF pulses preceded by a preparation sequence.
% Calculates time evolution along the RF train. Evolution for RF pulse train for an off resonance
% angle theta is displayed. M after the last RF pulse is compared to the theoretical steady-state signal.
% Mxy at TE = TR/2 is also shown.
% Two RF modes: chop, where the RF phase alternates by pi radians, and constant with constant RF phase.
% Preparation (prep) mode can be half_angle (alpha/2), inversion_half_angle (inversion + alpha/2), opt (optimal
% angle calculated from the steady-state Equations which is close to alpha/2) and inversion_opt (inversion + opt).
% 09-August-2003, Yuval Zur.

% Set parameters.
% =========================================================================================

% alp = 60;       % RF flip angle.
T1 = 300;
% T2 = 90;
% TR = 10;
T2 = 300;
TR = 2;
alp = 7.68;
alp = alp/180*pi;
M = 201;		    % Number of off resonance angles.
N = 5000;	    % Number of RF pulses.
time_vec = (0:N - 1)*TR + TR/2;
theta = pi/8;      % Off resonance angle where time evolution is calculated in radians.
th = (0:M - 1)*2*pi/(M - 1) - pi;

plot_Mxy = 0;   % Flag to plot Mxy for all RF pulses.

% Set prep and RF pulse parameters.
% RF_axis = 'x';
% RF_phase = pi/4;
RF_phase = 0;
% RF_mode = 'chop';   % 0 pi alternating RF phase.
RF_mode = 'constant';     % Constant RF phase.

if ~strcmp(RF_mode, 'chop') && ~strcmp(RF_mode, 'constant')
    error(' Illegal RF_mode. ');
end
% Preparation mode.
% prep = 'half_angle';        % alp/2.
% prep = 'inversion_half_angle';    % Inversion + alp/2.
% prep = 'opt';             % beta, where beta is calculated from the steady state Mxy/Mz.
% prep = 'inversion_opt';   % inversion + beta.
prep = 'noop';            % No peparation.

% ==========================================================================================

% Set TR1, the time between prep RF pulse and the start of the RF train.
if strcmp(RF_mode, 'constant')
    TR1 = TR;
else
    TR1 = TR/2;
end

% The off resonance angles. If alternating RF pulses shift by pi radians.
if strcmp(RF_mode, 'constant')
    Th = th;
else
    Th = th + pi;
end

[jj1, j1] = min(abs(th - theta));   % Find the index of theta.
[jj1, j2] = min(abs(Th - pi));      % Find the optimal index where beta is calculated from Mxy and Mz.

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

% Calculate signal from known Equations.
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;

% % S = sqrt(-1)*(a*exp(sqrt(-1)*Th) + b)./(c*cos(Th) + d);
% S = exp(sqrt(-1)*(RF_phase + pi/2))*(a*exp(-sqrt(-1)*Th) + b)./(c*cos(Th) + d);
S = exp(sqrt(-1)*(RF_phase + pi/2))*(a*exp(sqrt(-1)*Th) + b)./(c*cos(Th) + d);

MxMy = -E2*sin(Th)./(1 - E2*cos(Th));
% Theoretical Mz.
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*sin(alp)*E2 + (1 - E1)*cos(alp))/(1 - E1*cos(alp));

beta = atan(abs(S(j2))/Zm(j2));    % Steady state M angle. Prep RF pulse should push M to beta.

% Set alpha and beta of RF pulses.
% if ischar(RF_axis) && strcmp(RF_axis, 'x');
%     bRF = i*sin(alp/2);
%     bRF_half = i*sin(alp/4);
%     bRF_opt = i*sin(beta/2);
% elseif ischar(RF_axis) && strcmp(RF_axis, 'y');
%     bRF = -sin(alp/2);
%     bRF_half = -sin(alp/4);
%     bRF_opt = -sin(beta/2);
% else error(' Unknown RF pulse axis. ');
% end

bRF = sqrt(-1)*exp(sqrt(-1)*RF_phase)*sin(alp/2);
bRF_half = sqrt(-1)*exp(sqrt(-1)*RF_phase)*sin(alp/4);
bRF_opt = sqrt(-1)*exp(sqrt(-1)*RF_phase)*sin(beta/2);

aRF = cos(alp/2);
aRF_half = cos(alp/4);
aRF_opt = cos(beta/2);

% Define matrices.
mxy = zeros(3, M);
mecho = zeros(3, M);
Mfull = zeros(3, M, N);

% Set RF pulses matrices.
RF1 = Rot(aRF, bRF);
RF2 = Rot(aRF, -bRF);

% Relaxation matrices.
Rel = [E2 0 0; 0 E2 0; 0 0 E1];
Rel_prep = [E21 0 0; 0 E21 0; 0 0 E11];

v0 = [0 0 1].';			% start from z magnetization.

% Define and apply prep sequence.
if exist('prep', 'var') && ~strcmp(prep, 'noop')
    if ~strcmp(prep, 'inversion_opt') && ~strcmp(prep, 'opt') && ~strcmp(prep, 'half_angle') && ...
            ~strcmp(prep, 'inversion_half_angle')
        error(' Illegal prep option. ');
    end
    if strcmp(prep, 'inversion_opt') || strcmp(prep, 'opt')
        aRF_prep = aRF_opt; bRF_prep = bRF_opt;
    elseif  strcmp(prep, 'inversion_half_angle') || strcmp(prep, 'half_angle')
        aRF_prep = aRF_half; bRF_prep = bRF_half;
    end
    if strcmp(prep, 'inversion_half_angle') || strcmp(prep, 'inversion_opt')
        v0 = -v0;
    end
    if strcmp(RF_mode, 'constant')
        v0 = Rel_prep*Rot(exp(-sqrt(-1)*th(j1)*TR1/TR/2), 0)*Rot(aRF_prep, bRF_prep)*v0 + [0 0 (1 - E11)].';
    else
        v0 = Rel_prep*Rot(exp(-sqrt(-1)*th(j1)*TR1/TR/2), 0)*Rot(aRF_prep, -bRF_prep)*v0 + [0 0 (1 - E11)].';
    end
end

% Run the RF pulse train.
vi = RF1*v0;	% Apply first RF pulse on v0.
Mfull(:, :, 1) = repmat(vi, 1, M);

tic;
for jj = 1:M		% Theta loop.
%     afree = exp(-sqrt(-1)*th(jj)/2);
afree = exp(sqrt(-1)*th(jj)/2);
    free = Rel*Rot(afree, 0);
    v1 = vi;
    if strcmp(RF_mode, 'constant')
        for k = 2:N	% RF pulse loop.
            v1 = RF1*(free*v1 + [0 0 (1 - E1)].');
            Mfull(:, jj, k) = v1;   % Store M for all Th and RF pulses in the train.
        end
    elseif strcmp(RF_mode, 'chop')
        for k = 2:N     % RF pulse loop.
            if mod(k, 2)
                v1 = RF1*(free*v1 + [0 0 (1 - E1)].');
            else
                v1 = RF2*(free*v1 + [0 0 (1 - E1)].');
            end
            Mfull(:, jj, k) = v1;   % Store M for all Th and RF pulses in the train.
        end
    else
        error(' Illegal RF_mode. ');
    end
    mxy(:, jj) = v1;        % Store M of the last RF pulse in the train for all th.
    mecho(:, jj) = sqrt(E2)*Rot(exp(-1j*th(jj)/4), 0)*v1;		% Echo signal TE/2 after the RF pulse.
end
toc;

Mtime = squeeze(Mfull(:, j1, :));   % M for all RF pulses at th = theta.
if strcmp(RF_mode, 'chop')         % Multiply by -+1.
    pm = ones(1, N); pm(2:2:end) = -1;
    Mtime(1, :) = pm.*Mtime(1, :);
    Mtime(2, :) = pm.*Mtime(2, :);
    mxy(1:2, :) = mxy(1:2, :)*pm(end);
    mecho(1:2, :) = mecho(1:2, :)*pm(end);
end
Mxy = mxy(1, :);
Mz = mxy(3, :);
Mxyecho = mecho(1, :);

% Plot data.
figure; subplot(2, 1, 1);
plot(th, real(Mxy), th, imag(Mxy)); grid;
legend(' real ', ' imaginary ');
title(strrep([' Calculated Mxy: RF_mode = ' RF_mode '. T1 = ' num2str(T1) ', T2 = ' num2str(T2) ', TR = ' num2str(TR) '. '], '_', '\_'));
subplot(2, 1, 2);
plot(th, real(S), th, imag(S)); grid;
title(' Theoretical Mxy ');
legend(' real ', ' imaginary ');
xlabel(' Theta, radians ');
figure;
plot(th, real(Mxyecho), th, imag(Mxyecho)); grid;
xlabel(' Theta, radians ');
title(strrep([' Signal at TE = TR/2: RF_mode = ' RF_mode  '. T1 = ' num2str(T1) ', T2 = ' num2str(T2)...
        ', TR = ' num2str(TR) '. '], '_', '\_'));
figure; plot(th, Mz, th, Zm, '+'); grid;
legend(' Calculated Mz ', ' Theoretical Mz ');
xlabel(' Theta, radians ');
ylabel(' Mz ');
title(strrep([' Mz: RF_mode = ' RF_mode '. T1 = ' num2str(T1) ', T2 = ' num2str(T2) ', TR = ' num2str(TR) '. '], '_', '\_'));

% Plot M for all RF pulses at th = theta.
figure; plot(1:N, real(Mtime(1, :)), 'b', 1:N, imag(Mtime(1, :)), 'r', 1:N, Mtime(3, :), 'g'); grid;
legend(' Mx ', ' My ', ' Mz ');
% title(strrep([' Theta = ' num2str(th(j1)) ' rads. RF_mode = ' RF_mode '. RF_axis = ' RF_axis '. prep = ' prep '. '], '_', '\_'));
title(strrep([' \phi = ' num2str(th(j1)) ' rads. RF_mode = ' RF_mode '. RF phase = ' num2str(RF_phase/pi*180, 2) ' deg. prep = ' prep '. '], '_', '\_'));
xlabel(' RF pulse number ');

if exist('plot_Mxy', 'var') && (plot_Mxy)
    % Plot absolute Mxy for all RF pulses.
    figure; plot(1:N, abs(Mtime(1, :)), 'b', 1:N, Mtime(3, :), 'g'); grid;        % Magnitude display.
    legend(' Mxy', ' Mz ');
%     title(strrep([' Theta = ' num2str(th(j1)) ' rads. RF_mode = ' RF_mode '. RF_axis = ' RF_axis '. prep = ' prep '. '], '_', '\_'));
    title(strrep([' \phi = ' num2str(th(j1)) ' rads. RF_mode = ' RF_mode '. RF phase = ' num2str(RF_phase/pi*180, 2) ' deg. prep = ' prep '. '], '_', '\_'));
    xlabel(' RF pulse number ');
end

disp(' ');
disp(' Parameters ');
disp(' +++++++++++ ');
disp([' T1 = ' num2str(T1) '. ']);
disp([' T2 = ' num2str(T2) '. ']);
disp([' TR = ' num2str(TR) '. ']);
disp([' RF flip angle (degrees) = ' num2str(alp/pi*180) '. ']);
disp([' Optimal preparation angle = ' num2str(beta/pi*180) ' degrees. ']);
disp([' RF mode = ' RF_mode '. ']);
disp([' prep = ' prep '. ']);
disp([' Off resonance angle = ' num2str(theta/pi*180) ' degrees. ']);
disp(' ');
