%close all
clear
Which_case=1;
%% User options
softening = 'wes';  % 'wet' (no softening) or 'wes' (with softening)
folder = ['./', softening, '/output_pwp_timePlot']; % './wet/results_00/output_pwp00'; % Directory of the text files %USE TO EXTRACT THE PWP AT T=10 & FINAL STEP
%folder = ['./', softening, '/output_pwp'];% './wet/output_pwp'; % Directory of the text files %USE TO EXTRACT THE PWP AT T=14 & 20
ru_or_lambda=1; %1=ru=pw/sigma_z    ¦OR¦    2=lambda=pw/(tr(sigma)/3)
Visualization_pwp_or_lambda_or_stress = 2; %1=pwp    ¦OR¦    2=ru/lambda(depending on variable ru_or_lambda)c3 = sigma_zz    ¦OR¦    4 = sigma_mean=(s11+s22+s33)/3    ¦OR¦    5 = sigma_mean=(s11'+s22'+s33')/3    ¦OR¦    6 = bed displacement    ¦OR¦    7=sigma_yy,sigma_zz,sigma_xx


threshold_erosion = -0.01; % Threshold displacement (m) for material to be eroded. Negative to count all particles
particle_size = 0.005; % [m]
bed_volume = (1.38*0.24);

%%
files = dir(fullfile(folder, '*.txt')); % Get a list of all text files in the folder

% Initialize properties of all files (size=2*simulations, because for each sim I have both bed and flow)
sim_type = {};
data_type = {};
Phi = [];
k_Dp = [];
psi = [];
material = [];
time = [];
simulation_vis_sim_index = {};
sim_data_bed_14 = cell(1, length(files));
sim_data_bed_20 = cell(1, length(files));

% Initialize arrays for each unique simulation:
sim_type_sim = {};
Phi_sim = [];
k_Dp_sim = [];
psi_sim = [];
time_sim = [];

% Initialize cell arrays to store data for each simulation:
sim_data = cell(1, length(files)); %all simulations*2 (bot bed and flow)
sim_index=0; %unique simulation index

% Post-processed Output:
%x_runout = cell(1, length(files)/8); %the runout of all particles for each simulation (flow: all particles are counted. bed: particle is counted only if eroded)
RunoutSim = zeros(1, length(files)); %unique runout distance for each simulation. Extracted from x_runout based on the defined x_percentile
avg_lambda_14 = zeros(1, length(files));
avg_lambda_20 = zeros(1, length(files));
Relaxation = zeros(1, length(files));
x_percentile = 99;
ErosionVolSim = zeros(1, length(files)); %unique erosion volume (m^3/m) for each simulation.

% Regular expression pattern to extract properties:
pattern = '(dry|wet|wes|wec)_(\d+)_k(-?\d+)_d(-?\d+)_(particlesX|particlesDispl|particlespwp|particlesStress)_m(\d+)_t(\d+)';%(particlesX|particlesDispl)

%% Extract the particle positions for each simulation, for both bed and flow:
for i = 1:length(files)
    % Extract the filename without extension
    [~, name, ~] = fileparts(files(i).name);
    % Use regex to extract properties
    tokens = regexp(name, pattern, 'tokens');
    if ~isempty(tokens)
        % Assign properties to specific arrays
        sim_type{end+1} = tokens{1}{1};
        Phi(end+1) = str2double(tokens{1}{2});
        k_Dp(end+1) = str2double(tokens{1}{3});
        psi(end+1) = str2double(tokens{1}{4});
        data_type{end+1} = tokens{1}{5};
        material(end+1) = str2double(tokens{1}{6});
        time(end+1) = str2double(tokens{1}{7});
        filename = fullfile(folder, files(i).name); % Full path to the file
        fileContent = fileread(filename); % Read the entire file content
        modifiedContent = strrep(fileContent, '  ', ' '); % Replace double spaces with a single space
        % Write the modified content to a temporary file
        tempFilename = fullfile(folder, 'temp.txt');
        fileID = fopen(tempFilename, 'w');
        fwrite(fileID, modifiedContent);
        fclose(fileID);

        if strcmp(data_type{i}, 'particlesX')
            dataMatrix = readmatrix(tempFilename, 'Delimiter', ' ', 'HeaderLines', 2);  %, 'HeaderLines', 2  -> number of header lines
            sim_data{i} = struct('particleID',dataMatrix(:, 4), 'x_pos', dataMatrix(:, 5), 'y_pos', dataMatrix(:, 6), 'z_pos', dataMatrix(:, 7));
        elseif strcmp(data_type{i}, 'particlesDispl')
            dataMatrix = readmatrix(tempFilename, 'Delimiter', ' ', 'HeaderLines', 3);  %, 'HeaderLines', 3  -> number of header lines
            sim_data{i} = struct('particleID',dataMatrix(:, 4), 'x_displ', dataMatrix(:, 5), 'y_displ', dataMatrix(:, 6), 'z_displ', dataMatrix(:, 7));
        elseif strcmp(data_type{i}, 'particlespwp')
            dataMatrix = readmatrix(tempFilename, 'Delimiter', ' ', 'HeaderLines', 2);  %, 'HeaderLines', 3  -> number of header lines
            sim_data{i} = struct('particleID',dataMatrix(:, 4), 'pwp',dataMatrix(:, 5));
        elseif strcmp(data_type{i}, 'particlesStress')
            dataMatrix = readmatrix(tempFilename, 'Delimiter', ' ', 'HeaderLines', 2);  %, 'HeaderLines', 3  -> number of header lines
            sim_data{i} = struct('particleID',dataMatrix(:, 4), 's11', dataMatrix(:, 5), 's12', dataMatrix(:, 6), 's13', dataMatrix(:, 7), 's21', dataMatrix(:, 8), 's22', dataMatrix(:, 9), 's23', dataMatrix(:, 10), 's31', dataMatrix(:, 11), 's32', dataMatrix(:, 12), 's33', dataMatrix(:, 13));
            %it is sigma' !!!!
        end

        % After reading, delete the temporary file
        delete(tempFilename);
        % Check if data was successfully read
        if isempty(dataMatrix)
            warning('No data extracted for file: %s', files(i).name);
            continue;
        end

    else
        warning('No properties extracted for file: %s. Skipping...', files(i).name);
    end
end

%% Append, for the same simulation, "particlesX", "particlespwp", "particlesStress" on to "particlesDispl", so that it becomes easier to work with one file only
for i = 1:length(sim_data)
    % Check if the current index is 'particlesDispl'
    if strcmp(data_type{i}, 'particlesDispl')        
        % Now find matching indices for 'particlesX', 'particlespwp', and 'particlesStress'
        for j = 1:length(sim_data)
            if i ~= j && time(i) == time(j) && strcmp(sim_type{i}, sim_type{j}) && Phi(i) == Phi(j) && k_Dp(i) == k_Dp(j) && psi(i) == psi(j) && material(i) == material(j)
                % Check for unique simulation, matching properties and different material (bed + flow), but having different data_type=particlesX OR particlesDispl
                % Check and append data from 'particlesX' to 'particlesDispl'
                if strcmp(data_type{j}, 'particlesX')
                    sim_data{i}.x_pos = sim_data{j}.x_pos;
                    sim_data{i}.y_pos = sim_data{j}.y_pos;
                end
                % Check and append data from 'particlespwp' to 'particlesDispl'
                if strcmp(data_type{j}, 'particlespwp')
                    sim_data{i}.pwp = sim_data{j}.pwp;% + 101325 - sim_data{i}.pressCC ; %- atm_pwp;
                    %sim_data{j}.pressCC;
                end
                % Check and append data from 'particlesStress' to 'particlesDispl'
                if strcmp(data_type{j}, 'particlesStress')
                    sim_data{i}.s11 = sim_data{j}.s11;
                    sim_data{i}.s22 = sim_data{j}.s22;
                    sim_data{i}.s33 = sim_data{j}.s33;
                end
            end
        end
    end
end

%% correct all the pwp:
for i = 1:length(sim_data)
    % Check if the current index is 'particlesDispl'  --> where I have stored all my data
    if strcmp(data_type{i}, 'particlesDispl')        
        % Now find matching indices for 'particlesX', 'particlespwp', and 'particlesStress'
        
        num_bins = 1; % Number of bins you want to divide your x-axis into
        x_min = min(sim_data{i}.x_pos);
        x_max = max(sim_data{i}.x_pos);
        bin_edges = linspace(x_min, x_max, num_bins+1);
        minpwp = zeros(1, num_bins); %pwp at ymax (top bed particle) of each bin
        
        %Alternative 2: subtract from press_CC, taken a bit above the highest particle
        for k = 1:num_bins
            % Find indices of x_pos that fall into the current bin
            bin_indices = sim_data{i}.x_pos >= bin_edges(k) & sim_data{i}.x_pos < bin_edges(k+1);

            if any(bin_indices)
                xvalue = round((0.75 + 0.25)*100); % specify the x value .75  1.4   ------  0.75  ------
                yvalue = round((0.2 + 0.35)*100); % specify the y value 0.3   .2    ------  0.2   ------

                pressCC_value = PressCCExtractor_new_time(Which_case,softening, sim_type{i},Phi(i),k_Dp(i),psi(i),time(i), xvalue, yvalue); %NOW it will return only the value of press_CC for that particular simulation, at the specific location xvalue,yvalue
                sim_data{i}.pwp(bin_indices) = sim_data{i}.pwp(bin_indices) + 101325 - pressCC_value; %finally correct the p.pressure for all the MPs within the bin
            end
        end
    end
end


%% Calculation erosion volume & lambda
% Define a structure to hold results for unique simulations
% Define a function to create a valid field name


% Define a structure to hold results for unique simulations
uniqueSimulations = struct();

for i = 1:length(sim_data)
    if strcmp(data_type{i}, 'particlesDispl')
        % Calculate the necessary values
        ErosionVolSim_14(i) = sum(sqrt(sim_data{i}.x_displ.^2 + sim_data{i}.y_displ.^2) > threshold_erosion) * particle_size^2 / bed_volume;
        displacementMagnitudes_14 = sqrt(sim_data{i}.x_displ.^2 + sim_data{i}.y_displ.^2);
        nonZeroDisplacementIndices_14 = displacementMagnitudes_14 > threshold_erosion;
        displ_bed_14(:,i) = displacementMagnitudes_14;
        
        switch ru_or_lambda
            case 1
                fullExpressionLambda14 = sim_data{i}.pwp(nonZeroDisplacementIndices_14) ./ (( - sim_data{i}.s22(nonZeroDisplacementIndices_14) ) + 1e-6 + sim_data{i}.pwp(nonZeroDisplacementIndices_14));
            case 2
                fullExpressionLambda14 = sim_data{i}.pwp(nonZeroDisplacementIndices_14) ./ ((-sim_data{i}.s11(nonZeroDisplacementIndices_14) - sim_data{i}.s22(nonZeroDisplacementIndices_14) - sim_data{i}.s33(nonZeroDisplacementIndices_14)) ./ 3 + 1e-6 + sim_data{i}.pwp(nonZeroDisplacementIndices_14));
        end
        
        if ~isempty(fullExpressionLambda14)
            Lambda_bed_points_14(:,i) = fullExpressionLambda14(:);
        else
            Lambda_bed_points_14(:,i) = 0; 
        end
        avg_lambda_14(i) = median(Lambda_bed_points_14(:,i));

        % Define the unique key for the current simulation
        uniqueKey = createValidFieldName(Phi(i), k_Dp(i), psi(i));

        % Initialize storage if this is the first occurrence of the key
        if ~isfield(uniqueSimulations, uniqueKey)
            uniqueSimulations.(uniqueKey).avg_lambda = [];
            uniqueSimulations.(uniqueKey).time = [];
        end
        
        % Append the current values to the corresponding fields
        uniqueSimulations.(uniqueKey).avg_lambda = [uniqueSimulations.(uniqueKey).avg_lambda, avg_lambda_14(i)];
        uniqueSimulations.(uniqueKey).time = [uniqueSimulations.(uniqueKey).time, (time(i)-10)/10];
    end
end

% Display the results for each unique simulation
simulationFields = fieldnames(uniqueSimulations);
for j = 1:length(simulationFields)
    key = simulationFields{j};
    fprintf('Simulation %s:\n', key);
    disp(uniqueSimulations.(key));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%if strcmp(data_type{i}, 'particlesDispl') && Phi(i) && k_Dp(i) == k_Dp(j) && psi(i) == psi(j) time(i)
%% Define the specific simulations to plot
simulationsToPlot = [
    30, -3, -10;
    30, -5, -10;
    30, -3, 10;
    30, -5, 10;
    % Add more simulations as needed
];

colors = {'#FF9B69', '#FF0000', '#6EEC77', '#31AC24'}; % Add more colors if needed

% Create a cell array for the unique keys of the simulations to plot
plotKeys = cell(size(simulationsToPlot, 1), 1);
for idx = 1:size(simulationsToPlot, 1)
    Phi = simulationsToPlot(idx, 1);
    k_Dp = simulationsToPlot(idx, 2);
    k_Dp_sim_mapped = 10^k_Dp;
    psi = simulationsToPlot(idx, 3);
    plotKeys{idx} = createValidFieldName(Phi, k_Dp, psi);
end

% Plot the lambda evolution with time for the selected simulations
f3=figure;
dpi=600; lszs=1; 
box on;
set(f3,'Units','centimeters'); set(f3,'PaperUnits','centimeters'); 
pos=[1 1 9.5 8]; 
set(f3,'InnerPosition',pos); ax3.PositionConstraint = "outerposition"; %prevents text to get cut :)
set(gca,'linewidth',.8*lszs)
FontSize=9;%/72*dpi  *2.54;%Use this!!!=) (pt * in/pt * px/in)
ax3=gca; 


hold on;
legendEntries = cell(size(simulationsToPlot, 1), 1);
for idx = 1:length(plotKeys)
    key = plotKeys{idx};
    if isfield(uniqueSimulations, key)
        time_values = uniqueSimulations.(key).time;
        lambda_values = uniqueSimulations.(key).avg_lambda;

        % Plot line
        plot(time_values, lambda_values, '-', 'LineWidth', 1.2*lszs, 'Color', colors{idx});
        legendEntries{idx} = sprintf('\x03C8_{b0} = %.0f°, k_{b0} = 10^{%.0f} m^2', simulationsToPlot(idx, 3), 2*simulationsToPlot(idx, 2)-3);

    else
        warning('Simulation %s not found in the data.', key);
    end
end
% scatter(psi_sim(avg_lambda_14>-999), avg_lambda_14(avg_lambda_14>-999), 40, k_Dp_sim_mapped(avg_lambda_14>-999), 'filled', 'Marker', 'o', 'LineWidth', 0.8*lszs, 'MarkerEdgeColor', 'k');
% 
% my_colormap = [0.812, 0.894, 0.953
%                0.275, 0.600, 0.812
%                0.09803921568627451, 0.3607843137254902, 0.5725490196078431];
hold off;
xlabel('Time (s)');
ylabel('Bed pore pressure ratio, \lambda_b');
ylim([0 1])
xtickformat('%.1f'); ytickformat('%.1f');
legend(legendEntries, 'Position', [0.149547963660841,0.140604841741815,0.513940520446097,0.309602649006622],FontSize=8.5);
set(gca,'linewidth',.8*lszs)

location=sprintf('./%s/Figures/pwp_evolution.png',softening);
saveas(f3,location)
location=sprintf('./%s/Figures/pwp_evolution.svg',softening);
saveas(f3,location)

% Plot the lambda evolution with time for all simulations
% Define a set of colors
% Define a set of colors using the hsv colormap
% colors = hsv(length(simulationFields));
% 
% % Define a set of markers to use
% markers = {'o', 's', 'd', '^', 'v', '>', '<', 'p', 'h', 'x', '*', '+', '.', '|', '_', 'square', 'diamond', 'pentagram'};
% 
% figure;
% hold on;
% for j = 1:length(simulationFields)
%     key = simulationFields{j};
%     time_values = uniqueSimulations.(key).time;
%     lambda_values = uniqueSimulations.(key).avg_lambda;
% 
%     % Select color and marker
%     color = colors(j, :);
%     marker = markers{mod(j-1, length(markers)) + 1};
% 
%     % Plot with different colors and markers
%     plot(time_values, lambda_values, '-', 'DisplayName', key, 'Color', color, 'Marker', marker);
% end
% hold off;
% xlabel('Time');
% ylabel('Average Lambda');
% title('Lambda Evolution with Time for All Simulations');
% legend('show');
% grid on;



%%
function validKey = createValidFieldName(Phi, k_Dp, psi)
    validKey = sprintf('Phi_%g_kDp_%g_psi_%g', Phi, k_Dp, psi);
    validKey = strrep(validKey, '-', 'n'); % Replace '-' with 'n'
end


% Function to determine marker face color based on k_Dp value
function color = getMarkerFaceColor(k_Dp)
    % Define colors based on k_Dp values
    colorMap = [
        0.812, 0.894, 0.953;  % k_Dp = -3
        0.275, 0.600, 0.812   % k_Dp = -5
        % Add more colors as needed
    ];
    
    % Determine index based on k_Dp value
    if k_Dp == -3
        color = colorMap(1, :);
    elseif k_Dp == -5
        color = colorMap(2, :);
    else
        color = [0.5, 0.5, 0.5];  % Default color if not found
    end
end