Function PostProcessing

    %% PostProcessing
    clear
    close all
    
    %% Data information
    
    % The name of the pictures folder
    FolderName = 'save';
    
    % The name of the simulation (the output of the script Displacement 2D.m)
    % This is where the data are saved
    FileName = 'Test.mat';
    
    % The iteration limits
    % To have a faster script, you can focus on iterations where the shear band
    % changes (see pictures generated)
    ts = 15; % start
    te = 23; % stop
    
    % Set a threshold value
    % values lower or greater than this will be considered as error
    % Need to be set for every postprocessing
    min_max_eps = 0.5;
    
    %% Load the data
    
    s = load(strcat(FolderName,'/',... %Folder name
                    FileName)); % File name
    
    %% Read data and Initialization
    
    % Read data
    TimeStep = s.TimeStep ; % the time step used
    Maps = s.Displacement.Maps; % the mesh used
    Times = s.Displacement.Times; % the list of the time you used
    M1 = s.Displacement.Discretisation_X(1,:); % The x coordinates of the mesh
    M2 = s.Displacement.Discretisation_Y(1,:); % The x coordinates of the mesh
    
    % Initialization
    ListWidthX = zeros(4,te-ts+1); % 3 slices + 1 mean
    ListWidthY = zeros(4,te-ts+1); % 3 slices + 1 mean
    ListWidth = zeros(4,te-ts+1); % 3 slices + 1 mean
    
    %% Iteration on time
    
    for t = ts:te
    
        % Read strain maps
        e11 = s.Displacement.e11(:,:,t);
        e12 = s.Displacement.e12(:,:,t);
        e22 = s.Displacement.e22(:,:,t);
    
        % Build the map of X and Y
        X = zeros(size(Maps));
        Y = zeros(size(Maps));
        for l = 1 : size(Maps,1)
            for c = 1 : size(Maps,2)
                indice = (l-1)*size(Maps,2) + c ;
                X(l,c)= M1(indice);
                Y(l,c)= M2(indice);
            end
        end
    
        % Adapt strain maps with the threshold value
        for l = 1 : size(Maps,1)
            for c = 1 : size(Maps,2)
              if e11(l,c) < 0 || min_max_eps < e11(l,c)
                e11(l,c) = 0;
              end
              if e12(l,c) < 0 || min_max_eps < e12(l,c)
                e12(l,c) = 0;
              end
              if e22(l,c) < 0 || min_max_eps < e22(l,c)
                e22(l,c) = 0;
              end
            end
        end
    
        % Replot the figure
        NameFigure = ['Iteration ' int2str(t-ts+1)];
        figure('Name', NameFigure);
    
        % Plot e11
        subplot(131)
        surf(X,Y,e11,'EdgeColor', 'None', 'facecolor', 'interp')
        set(gca,'DataAspectRatio',[1,1,1])
        xlabel('x')
        ylabel('y')
        TitleName = ['\epsilon_{11} at time ' int2str(t-ts+1)];
        title(TitleName)
        colorbar
        view(2)
    
        % Plot e22
        subplot(132)
        surf(X,Y,e22,'EdgeColor', 'None', 'facecolor', 'interp')
        set(gca,'DataAspectRatio',[1,1,1])
        xlabel('x')
        ylabel('y')
        TitleName = ['\epsilon_{22} at time ' int2str(t-ts+1)];
        title(TitleName)
        colorbar
        view(2)
    
        % Plot e12
        subplot(133)
        surf(X,Y,e12,'EdgeColor', 'None', 'facecolor', 'interp')
        set(gca,'DataAspectRatio',[1,1,1])
        xlabel('x')
        ylabel('y')
        TitleName = ['\epsilon_{12} at time ' int2str(t-ts+1)];
        title(TitleName)
        colorbar
        view(2)
    
        saveas(gcf,strcat('png/pp_t_',int2str(t-ts+1),'.png'))
        close gcf
    
        % Build the list of x/y axis for the slices
        ListCutY = X(1,:); % list of x
        StepX = abs(ListCutY(2)-ListCutY(1));
        ListCutX = Y(:,1)'; % list of y
        StepY = abs(ListCutX(2)-ListCutX(3));
    
    %% Definition of the slides
    % At the first iteration, it is asked for the user to select 3 points
    % Those points define the slices where the algorithm is applied
    
        if t==ts
            % Here the overview is on e11
            % It can be on other strain map, just change
            % The same slices are used for the other strain maps
    
            % Create the figure
            ViewE11 = figure('Name','Overlook on \epsilon11');
            surf(X,Y,e11,'EdgeColor', 'None', 'facecolor', 'interp')
            set(gca,'DataAspectRatio',[1,1,1])
            xlabel('x')
            ylabel('y')
            TitleName = ['\epsilon_{11} at initial time'];
            title(TitleName)
            colorbar
            view(2)
            % wait for the input of the user
            uiwait(msgbox('You have to select 3 slides of \epsilon_{11}'))
    
            % Receive the information from the GUI
            [Cut1_x,Cut1_y] = ginput(1);
            [Cut2_x,Cut2_y] = ginput(1);
            [Cut3_x,Cut3_y] = ginput(1);
    
            close(ViewE11)
    
            % Look for the nearest lines and column
            indiceCutX1 = find(abs(ListCutY-Cut1_x)