Before plotting linear regression line to find shear center

final_lab.m

% distance that leading and trailing edge fell
% use arctan(LE TE) formula to find the angle and shear center

% strain vs applied load relationship

clc; clear all; close all;

data = readtable('LeadingEdge.xlsx');

zero = data(1,:);
data = data(2:8,:);
data.Voltage = data.Voltage-zero.Voltage;

y = data.Distance;
x = data.Voltage;

disp('Leading Edge')
[leading_edge_m, leading_edge_b] = linear_regression(x, y);

data = readtable('TrailingEdge.xlsx');

zero = data(1,:);
data = data(2:6,:);
data.Voltage = data.Voltage-zero.Voltage;

y = data.Distance;
x = data.Voltage;

disp('Trailing Edge')
[trailing_edge_m, trailing_edge_b] = linear_regression(x, y);

data = readtable('PotentiometerData.xlsx');
zero = data(1,:);
data = data(2:end,:);
data.TE_Mean = data.TE_Mean-zero.TE_Mean;
data.LE_Mean = data.LE_Mean-zero.LE_Mean;

LE_displacement = data.LE_Mean*leading_edge_m+leading_edge_b; %should I abs this?
TE_displacement = data.TE_Mean*trailing_edge_m+trailing_edge_b; %should I abs this?

hold on;
plot(str2double(data.TE), LE_displacement, '-o', 'LineWidth', 1.5)
plot(str2double(data.TE), TE_displacement, '-o', 'LineWidth', 1.5)
xlabel('Distance from LE that weight was applied at (in)')
ylabel('Displacement of potentiometer (in)')
title('Displacement vs Point of weight application')
legend('Leading Edge', 'Trailing Edge')
xticks([0.5 3 6 9 11.5])
hold off;

twist_angle = atan((LE_displacement-TE_displacement)/12)

function [m, b] = linear_regression(x, y)
    N = length(x);
    q = zeros(N,1);

    S1 = 1/(N-1)*(sum(x.^2)-1/N*(sum(x))^2);
    m = (N*sum(x.*y)-sum(x)*sum(y))/(N*(N-1)*S1);
    b = mean(y)-m*mean(x);

    gamma = 0.95;
    p = (1-gamma)/2;
    nu = N - 1;
    z = abs(tinv(p,nu));

    S1 = sqrt(1/(N-1)*(sum(x.^2)-1/N*(sum(x))^2));

    S2 = sqrt(1/(N-1)*(sum(y.^2)-1/N*(sum(y))^2));
    q0 = (N-1)*(S2^2-m^2*S1^2);

    SB = sqrt(q0/((N-2)*(N-1)*(S1^2)));
    k = z*SB;

    fprintf('95%% confidence interval for the slope is %f +- %f\n',m,k)

    SA = sqrt(SB^2/N*sum(x.^2));
    l = z*SA;

    fprintf('95%% confidence interval for the intercept is %f +- %f\n\n',b,l)
end