function plot_psic(input_data)

% Postprocessing script to get Power-law Sensitivity to Initial Condition (PSIC) curve
% Copyrighted by Jing Hu and Jianbo Gao, ECE department, University of Florida

% input_data: the data you want to apply PSIC to. Actually this is the
% output file from the computation of the Time-dependent exponent (TDE)
% curves. The C and Fortran codes for computing TDE can be downloaded from 
% http://www.gao.ece.ufl.edu/GCTH_Wileybook/programs/lambda_k_curves/
% The output file of TDE curves has the following format:
% the first row stores the number of pairs of vectors within each shell
% the rest rows store the information of \Lambda(k) curves
% each column corresponds to one shell used for computing TDE curves

% close all the open figure windows
%close all;

% get the number of pairs of vectors within each shell
num_of_pair = input_data(1, :);

% note even though here I labeled it as max shell, it actually corresponds
% to the smallest scale resolvable
max_shell = max(find(num_of_pair >0)); 

row_num = size(input_data, 1) - 1;

% plot the PSIC curves
plot(log(1:row_num), input_data(2:end, 1:max_shell), '.-');
set(gca, 'FontSize', 18);
xlabel('ln k', 'FontSize', 24);
ylabel('\Lambda(k)', 'FontSize', 24);
firstseries=get(gca,'Children');
firstseries=firstseries(1);
masxxx=get(firstseries,'XData');
masyyy=get(firstseries,'YData');
while (strcmp(questdlg('Провести линию наклона не по всем точкам?','','Да', 'Нет','Да'),'Да')==1)
        %axes(h2);
        [mouseX,mouseY]=ginput(2);
        xmin=min(mouseX(1),mouseX(2));
        xmax=max(mouseX(1),mouseX(2));
        nesessary_points_coords=find(masxxx>xmin);
        nesessary_pointsx=masxxx(nesessary_points_coords);
        nesessary_pointsy=masyyy(nesessary_points_coords);
        nesessary_points_coords2=find(nesessary_pointsx<xmax);
        nesessary_pointsx=nesessary_pointsx(nesessary_points_coords2);
        nesessary_pointsy=nesessary_pointsy(nesessary_points_coords2);
        coefs=polyfit(nesessary_pointsx,nesessary_pointsy,1);
        a1=coefs(1);
        a0=coefs(2);
        yyy_teor=a0+a1*masxxx;
        yyy_teor2=a0+a1*nesessary_pointsx;
        u2=nesessary_pointsy-yyy_teor2;
        R2=(var(nesessary_pointsy)-var(u2))/var(nesessary_pointsy);
        %axes(h2);
        hold on;
        plot(masxxx,yyy_teor,'r');
        title(['PSIC, H=' num2str(a1) ', R^2=' num2str(R2)]);
    end
