
function [c, rankH3, M] = fun_FFA(k)

    global N Mass
    
    ui = sqrt(-1);
    
    %%%%%%%%%%%%%%%

    H1 = [-1, ui; ui, 1]/2;
    H2 = [Mass-cos(k), -ui*sin(k); ui*sin(k), -Mass+cos(k)];
    H3 = [-1, -ui; -ui, 1]/2;

    rankH3 = rank(H3);
    
    I = eye(N);
    c = zeros(N+1,2*N+1);
    M = zeros(N,N,N,2*N-1);

    c(N+1,1) = 1;
    M(:,:,1,1) = I;

    c(N,1) = -trace(H1);
    c(N,2) = -trace(H2);
    c(N,3) = -trace(H3);

    for k=2:N   % w polynomial index

        % Precomputing M_k terms given C_N = 1

        % Sum in positive indexes j=i+k
        %   where i goes from i=-(k-1) to i=k-1 -> j=1 to j=2k-1

        % M_k = H(z)*M_{k-1} + c_{n-k+1}I = sum_{j=1}^{2*k-1} M_{k,j}*z^{j-k}

        % The term H(z)*M_{k-1} can be decomposed in the 3 H(z)
        % contributions times the sum over M_{k-1,i} from j=1 to j=2k-3
        % i.e. H_{-1}M_{k-1} = H_{-1}sum_{i=-(k-2)}^{i=k-2} M_{k,j}*z^{i-1}

        for i=1:2*k-1   % z polynomial index in M_k
            if i <= 2*k-3  
                    % Z Order comming from H_{-1}*M_{k-1,i}
                    %   min: i=-(k-2)-1 -> j=+1
                    %   max: i=(k-2)-1 -> j=2k-3
               AUX = H1*M(:,:,k-1,i); % i in M_{k-1} goes from i=1 to i=2k-3
               M(:,:,k,i) = M(:,:,k,i) + AUX;
            end
            if (i >= 2) && (i <= 2*k-2)
                    % Z Order comming from H_{0}*M_{k-1,i}
                    %   min: i=-(k-2) -> j=+2
                    %   max: i=(k-2) -> j=2k-2
               AUX = H2*M(:,:,k-1,i-1); % i in M_{k-1} goes from i=1 to i=2k-3
               M(:,:,k,i) = M(:,:,k,i) + AUX;
            end
            if  i >= 3
                    % Z Order comming from H_{1}*M_{k-1,i}
                    %   min: i=-(k-2)+1 -> j=+3
                    %   max: i=(k-2)+1 -> j=2k-1
               AUX = H3*M(:,:,k-1,i-2); % i in M_{k-1} goes from i=1 to i=2k-3
               M(:,:,k,i) = M(:,:,k,i) + AUX;
            end
            M(:,:,k,i) = M(:,:,k,i) + c(N-k+2,i)*I;
        end

        % Sum in positive indexes j=i+k+1
        %   where i goes from i=-k to i=k -> j=1 to j=2k+1

        % C_{N-k} = -1/k*Tr(H(z)*M_k) = sum_{j=1}^{2*k+1} C_{N-k,j}*z^{j-k}

        for i=1:2*k+1 % z polynomial index in C_{N-k}
            if i <= 2*k-1
                    % Z Order comming from H_{-1}*M_{k,i}
                    %   min: i=-(k-1)-1 -> j=+1
                    %   max: i=(k-1)-1 -> j=2k-1
               c(N-k+1,i) = c(N-k+1,i) - trace(H1*M(:,:,k,i))/k; % i in M_{k} goes from i=1 to i=2k-1
            end
            if (i>=2) && (i<=2*k)
                    % Z Order comming from H_{0}*M_{k,i}
                    %   min: i=-(k-1) -> j=+2
                    %   max: i=(k-1) -> j=2k
               c(N-k+1,i) = c(N-k+1,i) - trace(H2*M(:,:,k,i-1))/k; % i in M_{k} goes from i=1 to i=2k-1
            end
            if i >=3
                    % Z Order comming from H_{1}*M_{k,i}
                    %   min: i=-(k-1)+1 -> j=+3
                    %   max: i=(k-1)+1 -> j=2k+1
               c(N-k+1,i) = c(N-k+1,i) - trace(H3*M(:,:,k,i-2))/k; % i in M_{k} goes from i=1 to i=2k-1
            end  
        end

    end
    
end

