function [K_p] = TM_Analysis_function(data, n, lambda, criterion)
    
    global base;
    global p;        
     
    lambda = flip(lambda(end-n : end))';        
    bar3D_plot = [7:-1:0];
    lim_max = max(data(:));
    
    %%%%%%%%%%%%%%% ALOCATION OF THE MEMORY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    moment = zeros(n + 1, length(p));
    lin = zeros(2, length(p));
    K_p = zeros(length(p), 1);
    R2 = zeros(length(p), 1);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    %%%%%%%%%%%%%%%%%%%%%% UP-SCALING PROCEDURE %%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i = 1 : length(p)
        moment(end, i) = mean2(data.^p(i));
    end
    
    for j = n : -1 : 1        
        n_el = length(data);
        d1 = data(1 : base : n_el, 1 : base : n_el) + data(2 : base : n_el + 1, 1 : base : n_el);
        d2 = data(1 : base : n_el, 2 : base : n_el + 1) + data(2 : base : n_el + 1, 2 : base : n_el + 1);
        data = (d1 + d2)/4;
                
        %%%%%%%%%%%%%%%%%%%%%%% 3D BAR IMAGE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if isempty(find(bar3D_plot == j-1))==0;
            fig3D = figure;
            b = bar3(data, 1);
            set(b, 'LineWidth', 0.1);
            hold on;
            [X,Y] = meshgrid(0.5 : length(data)+0.5);
            F_mean = criterion.*ones(length(X));
            b1 = surf(X, Y, F_mean);
            set(b1, 'edgecolor', 'none');
            colormap('jet');
            colorbar;  
            for r = 1:length(b)
                zdata = get(b(r), 'ZData');
                set(b(r), 'CData', zdata);
                set(b(r), 'FaceColor', 'interp');
            end
            caxis([0 lim_max]);
            xlim([0.5 length(data)+0.5]);
            ylim([0.5 length(data)+0.5]);
            set(gca, 'DataAspectRatio', [repmat(max(diff(get(gca, 'XLim')), diff(get(gca, 'YLim'))), [1 2]) lim_max]) 
            zlim([0 lim_max]);
            view(-61, 50);
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
        for k = 1 : length(p)
            moment(j, k) = mean2(data.^p(k));
        end        
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%% DETERMINATION OF THE BEST FITTING LINEAR REGRESSION %%%%%%%%%%%%%
    for k = 1 : length(p)
       lin(:,k) = polyfit(log(lambda), log(moment(:,k)), 1); 
       r2 = corrcoef(polyval(lin(:,k), log(lambda)), log(moment(:,k)));
       R2(k) = r2(1,2);
       K_p(k) = lin(1,k);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%%%%%%%%%%%% PLOTING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    index = [find(p == 0.2); find(p == 1.5); find(p == 2); find(p == 3)];
    
    % Scaling of statistical moments of order p 
    figure;
    plot(log(lambda), log(moment(:, index(1))), 'o', 'LineWidth', 3, 'MarkerSize', 12);
    hold on;
    plot(log(lambda), log(moment(:, index(2))), '^', 'LineWidth', 3, 'MarkerSize', 12);
    hold on;
    plot(log(lambda), log(moment(:, index(3))), 's', 'LineWidth', 3, 'MarkerSize', 12);
    hold on;
    plot(log(lambda), log(moment(:, index(4))), 'd', 'LineWidth', 3, 'MarkerSize', 12);
    hold on;
    plot(log(lambda), lin(2, index(1)) + lin(1, index(1))*log(lambda), 'r--', 'LineWidth', 3);
    hold on;
    plot(log(lambda), lin(2, index(2)) + lin(1, index(2))*log(lambda), 'r--', 'LineWidth', 3);
    hold on;
    plot(log(lambda), lin(2, index(3)) + lin(1, index(3))*log(lambda), 'r--', 'LineWidth', 3);
    hold on;
    plot(log(lambda), lin(2, index(4)) + lin(1, index(4))*log(lambda), 'r--', 'LineWidth', 3);
    hold off;
    grid on;
    
    moment_index = log(moment(:,index));
    ylabel('log(<E_{\lambda}^p>)', 'FontSize', 16, 'FontWeight', 'bold');
    xlabel('log({\lambda})', 'FontSize', 16, 'FontWeight', 'bold');
    ylim([2*min(moment_index(:)) 1.1*max(moment_index(:))]);
    xlim([0 ceil(log(lambda(end)))]);
    
    % Emirical moment scaling function K(p)
    figure;
    plot(p, K_p, 'k-', 'LineWidth', 3);
    ylabel('K(p)', 'FontSize', 16, 'FontWeight', 'bold');
    xlabel('p', 'FontSize', 16, 'FontWeight', 'bold');
    ylim([1.5*min(K_p) 1.1*max(K_p)]);
    xlim([0 1.1*p(end)]);
    grid on;
end