%twoper
% Clear all variables and close all figures
clear; close all; clc;

% Parameters
headerSize = 11;  % Number of header lines in CSV file
filename = 'HV7 250 mmHg 3_1.csv';  % Name of the CSV file
%filename = '2% offset correct - use code_1.csv';  % Name of the CSV file

crossSectionalArea = 28.27433;  % Cross-sectional area in mm^2
initialHeight = 12;  % Initial height in mm
toeSlopeThreshold = 2000;  % Threshold for the slope to end the toe region
linearSlopeThreshold = 300;  % Threshold for the deviation to end the linear region

% Read the data from CSV file
opts = detectImportOptions(filename);
opts.DataLines = [headerSize+1, Inf]; % Adjust for header size
data = readtable(filename, opts);

% Assuming the data has columns: Displacement (mm) and Force (N)
displacement = data{:,3};
force = data{:,4};
% Shift displacement values to start at zero
displacement = displacement - displacement(1);

% Cut off data after 3 mm displacement
maxDisplacement = 3;  % Maximum displacement in mm
cutoffIdx = find(displacement > maxDisplacement, 1);
if ~isempty(cutoffIdx)
    displacement = displacement(1:cutoffIdx);
    force = force(1:cutoffIdx);
end

% Plot the force-displacement graph
figure;

plot(displacement, force, 'b-', 'LineWidth', 1.5);
xlabel('Displacement (mm)');
ylabel('Force (N)');
title('Force-Displacement Graph');
grid on;

% Calculate stress and strain
strain = displacement / initialHeight;  % Strain (dimensionless)
stress = (force*1000) / crossSectionalArea;    % Stress (N/mm^2 or MPa if area is in mm^2)

% Plot the stress-strain graph
figure;
plot(strain, stress, 'r-', 'LineWidth', 1.5);
xlabel('Strain (dimensionless)');
ylabel('Stress (MPa)');
title(filename(1:end-6));
grid on;

% Calculate the derivative of the stress-strain curve
dStress_dStrain = diff(stress) ./ diff(strain);

% Find the end of the toe region
toeRegionEndIdx = find(dStress_dStrain > toeSlopeThreshold, 1);

% Find the end of the linear portion
linearRegionEndIdx = find(abs(dStress_dStrain(toeRegionEndIdx:end) - mean(dStress_dStrain(toeRegionEndIdx:toeRegionEndIdx+10))) > linearSlopeThreshold, 1);
linearRegionEndIdx = linearRegionEndIdx + toeRegionEndIdx - 1; % Adjust index

% Ensure indices are within bounds
toeRegionEndIdx = max(toeRegionEndIdx, 1);
linearRegionEndIdx = min(linearRegionEndIdx, length(strain));

% Fit the linear portion
p = polyfit(strain(toeRegionEndIdx:linearRegionEndIdx), stress(toeRegionEndIdx:linearRegionEndIdx), 1);
youngsModulus = p(1);  % Young's Modulus (slope of linear portion)
x_intercept = -p(2) / p(1);  % X-intercept of the fitted line
disp(['Young''s Modulus: ', num2str(youngsModulus), ' MPa']);

% Calculate the fitted values for the linear portion
fittedStress = polyval(p, strain(toeRegionEndIdx:linearRegionEndIdx));

% Plot the polyfit line on the stress-strain graph
hold on;
plot(strain(toeRegionEndIdx:linearRegionEndIdx), fittedStress, 'b--', 'LineWidth', 1.5);
hold off;


% Define the regions
strain_region = strain(toeRegionEndIdx:linearRegionEndIdx);
stress_region = stress(toeRegionEndIdx:linearRegionEndIdx);

% Fit a linear polynomial to the data in the specified region
p = polyfit(strain_region, stress_region, 1);  % 1 indicates a linear fit

% Get the fitted values
fittedStress = polyval(p, strain_region);

% Plot the polyfit line on the stress-strain graph
hold on;
plot(strain_region, fittedStress, 'b--', 'LineWidth', 1.5);
hold off;

% Calculate R^2
stress_mean = mean(stress_region);
SST = sum((stress_region - stress_mean).^2);  % Total Sum of Squares
SSR = sum((stress_region - fittedStress).^2); % Residual Sum of Squares
R2 = 1 - (SSR/SST);

% Display the R^2 value
disp(['R^2 = ', num2str(R2)]);

% Find and display the start and end points on the x-axis
start_x = strain(toeRegionEndIdx);
end_x = strain(linearRegionEndIdx);
disp(['Start point on x-axis: ', num2str(start_x)]);
disp(['End point on x-axis: ', num2str(end_x)]);


% Find the peak force and peak stress
[peakForce, peakIdx] = max(force);
peakDisplacement = displacement(peakIdx);
peakStress = stress(peakIdx);
disp(['Peak Force: ', num2str(peakForce), ' N']);
disp(['Peak Stress: ', num2str(peakStress), ' MPa']);

% Find the 2% offset yield stress
offsetStrain = 0.02;
offsetLine = youngsModulus * (strain - (x_intercept + offsetStrain));  % Offset line
% Intersect offset line with the stress-strain curve
diff = stress - offsetLine;
[~, idx] = min(abs(diff));
offsetYieldStress = stress(idx);
offsetYieldForce = offsetYieldStress * crossSectionalArea;  % Convert to force
disp(['2% Offset Yield Stress: ', num2str(offsetYieldStress), ' MPa']);
disp(['2% Offset Yield Force: ', num2str(offsetYieldForce), ' N']);

% Plot the stress-strain graph
hold on;
h1 = plot(strain, stress, 'r-', 'LineWidth', 1.5);  % Stress-Strain Curve (red line)
h2 = plot(strain, offsetLine, 'k--', 'LineWidth', 1.5);  % 2% offset line (black dotted)
h3 = plot(strain(idx), stress(idx), 'go', 'MarkerSize', 8, 'MarkerFaceColor', 'g');  % 2% offset yield point (green circle)
h4 = plot(strain(peakIdx), peakStress, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');  % Peak stress (red circle)
h5 = plot(strain_region, fittedStress, 'b--', 'LineWidth', 1.5);  % Fitted linear portion (blue dotted)

% Update legend with the correct order
legend([h1, h5, h2, h3, h4], 'Stress-Strain Curve', 'Fitted Linear Portion', '2% Offset Line', '2% Offset Yield Point', 'Peak Stress');
hold off;

