%Script to create stress-strain curves from a CSV file
%Author: Victor Kojenov, vkojenov@pdx.edu

%Define material geometry (mm) and data file
datafile = 'sargentd_2_t3.csv';
l_initial = 25.4;
w_initial = 6.41;
t_initial = 1.25;

%Loading CSV files with columns [Time (min), Force (N), Position (mm)] into
%individual row vectors
data = csvread(datafile, 1, 0);   %Row offset is necessary for column labels
time = data(:,1)'./60;  %Convert to seconds
force = data(:,2)';
position = data(:,3)';

%Normalizing the force and position data into stress and strain
stress = force ./ (t_initial * w_initial);
strain = position ./ (position + l_initial);

%Plotting force and position
plot(strain, stress, '.')  
grid('on')
title('Stress-strain')
xlabel('Strain')
ylabel('Stress (MPa)')

%Have user select two values of the linear portion of the plot
[X, ~] = ginput(2);

%Find indices of the linear range based on the inputted x values
if X(1) > X(2)
    lin_index = find((strain < X(1)) & (strain > X(2)));
else
    lin_index = find((strain > X(1)) & (strain < X(2)));
end

%Simple linear fit of linear range based on max and min x values (effectively calculating modulus of elasticity)
lin_min = [strain(min(lin_index)), stress(min(lin_index))];
lin_max = [strain(max(lin_index)), stress(max(lin_index))];
elastic_modulus = (lin_max(2)-lin_min(2))/(lin_max(1)-lin_min(1));

%Finding x-intercept assuming perfect linear behavior to determine the offset due to jaw slippage
jaw_offset = lin_min(1)-(lin_min(2)/elastic_modulus);

%Correcting data for jaw slippage
strain = strain - jaw_offset;
new_origin = find(strain >= 0, 1, 'first');
time = time(1,new_origin:size(time,2));
stress = stress(1,new_origin:size(stress,2));
strain = strain(1,new_origin:size(strain,2));

%Plot force and position again to correct for sharp dropoff at the end
plot(strain, stress, '.')  
grid('on')
title('Stress-strain')
xlabel('Strain')
ylabel('Stress (MPa)')

%Have user select just above the first point that drops off sharply, and thereby find fracture strength
[X, Y] = ginput(1);
if X > max(strain)
    fracture_index = find(stress > Y, 1, 'last');
    sigma_fracture = stress(fracture_index);
    fracture_strain = strain(fracture_index);
else
    fracture_index = find((stress < Y) & (strain > X), 1, 'first') - 1;
    sigma_fracture = stress(fracture_index);
    fracture_strain = strain(fracture_index);
end

%Correct for dropoff
time = time(1:fracture_index);
stress = stress(1:fracture_index);
strain = strain(1:fracture_index);

%Calculating 0.2% offset line by minimizing error between linear line and actual data
lin_data = (strain - .002) .* elastic_modulus;
[err, sigma_y_index] = min(abs(lin_data - stress)); %find index with least error
sigma_y = stress(sigma_y_index);

%Calculating ultimate tensile strength and reduction of area
sigma_uts = max(stress);
reduction_area_fracture = 1/(1+fracture_strain); %using conservation of volume

%Report all crucial values
fprintf('0.2 percent offset yield strength: %.2f MPa\n', sigma_y)
fprintf('Ultimate tensile strength: %.2f MPa\n', sigma_uts)
fprintf('Percent elongation at fracture: %.2f\n', fracture_strain*100)
fprintf('Percent reduction of area at fracture: %.2f\n', reduction_area_fracture*100)
fprintf('Breaking stress: %.2f MPa\n', sigma_fracture)

%Final plot with line at 0.2% yield
yield_line_range = strain(1:sigma_y_index)-.002;
yield_line = elastic_modulus.*(yield_line_range-.002);
plot(strain, stress, '.', yield_line_range, yield_line)
title('Stress-strain')
xlabel('Strain')
ylabel('Stress (MPa)')
grid('on')
xlim([0 inf])
ylim([0 inf])