% Calculates the optimal Flip Angle for spectrum with least squares
% calculation method.
% ======================================================================
% time and sampling parameters
T1 = 1000;
T2 = 500;
% T1 = 10000;
% T2 = 2500;
TR = 3;
% FA = 13;
bw_kHz = 3;
N1 = 0; N2 = 0;
Nshift = 30;     % number of shifts.
band_num = Nshift/2;     % number of bands.
% band_num = Nshift - 5;     % number of bands.
% bw_kHz = floor(TR*bw_kHz)/TR;      % integer sampled points in TR.
points_in_TR = round(TR*bw_kHz);       % no of points from t = 0 to t = TR.
% dt = 1/points_in_TR;
% samp_points = points_in_TR - N1 - N2;       % points actually sampled.
% norm_time = (N1:samp_points - 1 + N1)*dt;
norm_time = (0:points_in_TR - 1)/points_in_TR;
% freq_points = points_in_TR;       % freq bins: total frequency bins = freq_points*band_num.
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(:)*ones(1, Nshift) + ones(fvals, 1)*(Nshift_vec*2*pi/Nshift);       % size: M1 by Nshift.
shifted_Ph = Ph(:) + (Nshift_vec*2*pi/Nshift);       % size: M1 by Nshift.
A = exp(-1j*norm_time(:)*Ph);
FA = 2:0.5:20;
var_val = zeros(1, length(FA));
std_dev = 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.
    % A = exp(-1j*norm_time(:)*Ph);
    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 = inv(Gn'*Gn);       % inv(Gn'*Gn) is the known variance of the ls parameters.
    var_val(k) = mean(diag(TempM));
    Res_ls = Gn\m_data;
    std_dev(k) = std(Res_ls);
end
[~, J] = min(var_val);
figure; plot(FA, var_val, '-o'); grid;
title([' optimal FA = ', num2str(FA(J)), ' deg. '])
