% Calculates the optimal Flip Angle for least squares calculation method.
function FA_opt = Find_FAopt(T1, T2, TR, Nshift)
bw_kHz = 3;
% Nshift = 40;     % number of shifts.
band_num = Nshift/2;     % number of bands.
% band_num = Nshift - 5;     % number of bands.
points_in_TR = round(TR*bw_kHz);       % no of points from t = 0 to t = TR.
norm_time = (0:points_in_TR - 1)/points_in_TR;
fvals = points_in_TR*band_num;
Ph = ((0:fvals - 1) - floor(fvals/2))*2*points_in_TR*pi/fvals;     % phase range [-pi, pi]*points_in_TR with fvals points.
Nshift_vec = (0:Nshift - 1) - floor(Nshift/2);
shifted_Ph = Ph(:) + (Nshift_vec*2*pi/Nshift);       % size: M1 by Nshift.
A = exp(-1j*norm_time(:)*Ph);
FA = 2:20;
var_val = zeros(1, length(FA));
for k = 1:length(FA)
    D = reshape(ssfp_signal(T1, T2, TR, FA(k), shifted_Ph(:), 0), fvals, []);       % size(D) = M1 by Nshift. Use T1, T2.
    % sigma = 2e-5;
    % m_data = A*D;       % measured data.
    % m_data1 = m_data(:);
    % m_data = m_data1 + sigma/sqrt(2)*complex(randn(size(A, 1)*Nshift, 1), randn(size(A, 1)*Nshift, 1));

    Gn = zeros(points_in_TR*Nshift, points_in_TR*band_num);

    for jj = 1:fvals
        tmp = A(:, jj).*D(jj, :);
        Gn(:, jj) = tmp(:);
    end

    TempM = pinv(Gn'*Gn);       % inv(Gn'*Gn) is the known variance of the ls parameters.
    var_val(k) = mean(diag(TempM));
end
var_val = abs(var_val);
% [~, J] = min(var_val);

FAint = linspace(FA(1), FA(end), 1000);
y = interp1(FA, var_val, FAint, 'spline');
FA_opt = FAint(y == min(y));
figure; plot(FA, var_val, 'o', FAint, y, 'm'); xlabel(' FA '); grid;
title([' FA\_opt = ', num2str(FA_opt, 3), ' deg. '])
