function [FA, fssfp, res, Nm] = find_flip(band_num, M, Nshift, T1, T2, TR, filt_struct, RF_phase, plot_data)
% input: M = number of points in filter. M = 300 is a good number.
% plot data = 1 will plot the optimal filters and the least squares
% difference between actual and optimal filters vs. flip angle FA.
% ====================================================================
if nargin < 7
    filt_struct.Nfilt = M/2;
    filt_struct.method = 'time';
end
if nargin < 8
    RF_phase = 0;
end
if nargin < 9
    plot_data = 0;
end

FA_vec = linspace(0.1, 20, 150);
Nfilt = filt_struct.Nfilt;
Nfilt = min(M, Nfilt) - 1 + mod(Nfilt, 2);     % Nfilt odd and <= M.

[filt, Y1, P1] = calc_filter(band_num, Nfilt, M);
Nshift_vec = (0:Nshift - 1) - floor(Nshift/2);
Nfilt_vec = (0:Nfilt - 1) - floor(Nfilt/2);
Phssfp = ((0:M - 1) - floor(M/2))*2*pi/M;
shifted_Phss = Phssfp(:) + Nshift_vec*2*pi/Nshift;
res = zeros(1, length(FA_vec));
Nm = zeros(1, length(FA_vec));
for j = 1:length(FA_vec) + 1
    if j > length(FA_vec)
        [~, K] = min(res);
        [~, J] = min(Nm);
        FA = FA_vec(K);
        FA_beta = FA_vec(J);
    else
        FA = FA_vec(j);
        FA_beta = FA_vec(j);
    end

    % compute alpha.
    % ==================
    ssfp = ssfp_signal(T1, T2, TR, FA, Phssfp, RF_phase);
    Aecho = 1/M*fftshift(fft(ifftshift(ssfp)));     % similar to Am in the paper.
    Dssfp = reshape(ssfp_signal(T1, T2, TR, FA, shifted_Phss(:), RF_phase), M, []);
    L = crop(Aecho, Nfilt, 1).*exp(1j*(Nfilt_vec.'*Nshift_vec)*2*pi/Nshift);
    %     lambda = 0;
    lambda = 1e-3;
    % lambda = 1e-2;
    if lambda > 0
        s = max(svd(L));
        alpha = pinv(L'*L + lambda*s*eye(size(L, 2)))*(L'*filt);
    else
        %     alpha = L\filt;
        alpha = pinv(L)*filt;
    end

    % compute beta
    % =========================
    lambda1 = 1e-3;
    %     lambda1 = 0;
    %     lambda1 = 1e-2;
    if lambda1 > 0
        s = max(svd(Dssfp));
        beta = pinv(Dssfp'*Dssfp + lambda1*s*eye(Nshift))*(Dssfp'*Y1);
    else
        beta = Dssfp\Y1;
        %     beta = pinv(Dssfp)*Y1;
    end

    % Compute SSFP filters
    % ===========================
    dc = floor(size(alpha, 2)/2) + 1;
    switch filt_struct.method
        case 'time'
            fssfp = Dssfp*alpha;
            alpha = alpha/fssfp(M/2 + 1, dc);
            fssfp = Dssfp*alpha;
            Normval = norm(alpha(:, dc));
        case 'phi'
            fssfp = Dssfp*beta;
            beta = beta/fssfp(M/2 + 1, dc);
            fssfp = Dssfp*beta;
            Normval = norm(beta(:, dc));
        otherwise
            FA = []; fssfp = []; res = [];
            errordlg(' find_flip: filter method is not defined.')
            return
    end
    if j <= length(FA_vec)
        %         dc = floor(size(alpha, 2)/2) + 1;
        %         res(j) = mean(abs(fssfp(:, dc) - Y1(:, dc)));
        res(j) = mean(abs(fssfp(:, dc) - Y1(:, dc)));
        Nm(j) = Normval;       % compute norm beta.
        %         res(j) = mean(abs(fssfp(:, dc) - zpad(ones(M/band_num, 1), M, 1)));
    end
end
if plot_data > 0
    h = figure;
    plot(FA_vec, res/max(abs(res)));
    grid on;
    hold on;
    plot(FA_vec, Nm/max(abs(Nm)), 'm');
    xlabel(' FA, deg ');
    ylabel(' res ');
    %     title(' ls difference between optimal and actual filters. ');
    title([' FA\_ opt = ', num2str(FA, 3), ' deg. FA\_beta = ', num2str(FA_beta, 3), ' deg. '])
    legend(' || C - L\beta || ', ' || \beta || ')
    set(h, 'Name', ' FA_ opt = optimal filter. FA_beta = best snr. ');
    %     plot(Phssfp(:), abs(fssfp));
    %     grid on;
    %     title([' optimal FA = ', num2str(FA, 3), ' deg. ']);
end
