function [FA, fssfp, res, Nm] = find_flip_new(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

