%{ 

Simple harmonic motion of oscillating horizontal bar 
Written by Dave Goeke, July 2019

Program takes in, and plots, five arrays of oscillations 
and wheel base
 
The linear regression and R-Squared analysis prove 
the formula for oscillation is the square root of two
times the coefficient of kinetic friction times acceleration
due to gravity divided by wheel base 

    w (omega or oscillation frequency) = 
    sqrt (2 * uK * g / L)

%}

clear all; close all

LargeBar_WB   = [ 2.6  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21] % Wheel base in inches
LargeBar_CC   = [59   59 56 42 42 41 31 30 29 27 25 22 19 17 16 13 17 20 18 21] % Oscillations in 30 seconds

ShortBar_WB   = [ 2.6  3  4  5  6  7  8  9] % Wheel base in inches
ShortBar_CC   = [63   64 48 40 37 35 30 32] % Oscillations in 30 seconds

WoodStick_WB  = [ 2.6  3  4  5  6  7  8  9] % Wheel base in inches
WoodStick_CC  = [46   45 34 28 27 25 26 25] % Oscillations in 30 seconds

ThinBar_WB    = [ 2.6  3  4  5  6  7  8] % Wheel base in inches
ThinBar_CC    = [51   47 44 38 41 34 32] % Oscillations in 30 seconds

PlasticRod_WB = [ 2.6  3  4  5  6  7  8  9 10 11 12] % Wheel base in inches
PlasticRod_CC = [57   50 43 37 36 34 31 29 33 28 30] % Oscillations in 30 seconds



figure(1)

subplot(221)
% Plot all four bars against. Oscillations per minute vs length of wheel
% base
plot(LargeBar_WB,2*LargeBar_CC,'-o','LineWidth',4) % Double the oscillation count so it is in minutes
hold on
grid on
plot(ShortBar_WB,2*ShortBar_CC,'-*g','LineWidth',4)
plot(WoodStick_WB,2*WoodStick_CC,'-+','LineWidth',4)
plot(ThinBar_WB,2*ThinBar_CC,'-x','LineWidth',4)
plot(PlasticRod_WB,2*PlasticRod_CC,'-x','LineWidth',4)
xlabel('Wheel base in inches','FontSize',23)
xticks([2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22])
ylabel('Oscillations per minute','FontSize',23)
legend({'Large Bar','Short Bar','Wood Stick','Thin Bar','Plastic Rod'},'Location','northeast','FontSize',23)
hold off



subplot(222)
% Plot short and long bar, same material. Oscillations per minute
plot(LargeBar_WB,2*LargeBar_CC,'-o','LineWidth',4)
hold on
grid on
plot(ShortBar_WB,2*ShortBar_CC,'-dg','LineWidth',4)
xlabel('Wheel base in inches','FontSize',23)
xticks([2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22])
ylabel('Oscillations per minute','FontSize',23)
legend({'Large Bar','Short Bar'},'Location','northeast','FontSize',23)



subplot(223)
% Plot all four bars. Oscillations per minute vs square root of wheel
% base
plot(sqrt(LargeBar_WB),2*LargeBar_CC,'-o','LineWidth',4)
hold on
grid on
plot(sqrt(ShortBar_WB),2*ShortBar_CC,'-*g','LineWidth',4)
plot(sqrt(WoodStick_WB),2*WoodStick_CC,'-+','LineWidth',4)
plot(sqrt(ThinBar_WB),2*ThinBar_CC,'-x','LineWidth',4)
plot(sqrt(PlasticRod_WB),2*PlasticRod_CC,'-x','LineWidth',4)
xlabel('Square root of wheel base','FontSize',23)
ylabel('Oscillations per minute','FontSize',23)
legend({'Large Bar','Short Bar','Wood Stick','Thin Bar','Plastic Rod'},'Location','northeast','FontSize',23)
hold off



subplot(224)
% Plot short and long bar, same material
plot(sqrt(LargeBar_WB),2*LargeBar_CC,'-o','LineWidth',4)
hold on
grid on
plot(sqrt(ShortBar_WB),2*ShortBar_CC,'-dg','LineWidth',4)
xlabel('Square root of wheel base','FontSize',23)
ylabel('Oscillations per minute','FontSize',23)
legend({'Large Bar','Short Bar'},'Location','northeast','FontSize',23)
hold off



figure(2)

subplot(321)  % Linear regression. Large Bar
              % ----------------------------              
% Remove last four data points 
LargeBar_WB   = [ 2.6  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 ]
LargeBar_CC   = [59   59 56 42 42 41 31 30 29 27 25 22 19 17 16 13 ]

c=polyfit(sqrt(LargeBar_WB),2*LargeBar_CC,1)
y_est=polyval(c,sqrt(LargeBar_WB))

y_residual = 2*LargeBar_CC-y_est;
SSresid = sum(y_residual.^2)
SStotal = (length(2*LargeBar_CC)-1) * var(2*LargeBar_CC)
rsq=1-SSresid/SStotal

plot(sqrt(LargeBar_WB),2*LargeBar_CC,'-o','LineWidth',4)
title(['Linear Regression, Large Bar. R-Squared = ' num2str(round(rsq,3))],'FontSize',18)
xlabel('Square root of wheel base','FontSize',14)
ylabel('Oscillations per minute','FontSize',14)

hold on
grid on
plot(sqrt(LargeBar_WB),y_est,'r--')
legend({'Oscillation Rate','Linear Best Fit'},'Location','northeast','FontSize',18)
hold off



subplot(322)  % Linear regression. Short Bar
              % ----------------------------
c=polyfit(sqrt(ShortBar_WB),2*ShortBar_CC,1)
y_est=polyval(c,sqrt(ShortBar_WB))

y_residual = 2*ShortBar_CC-y_est;
SSresid = sum(y_residual.^2)
SStotal = (length(2*ShortBar_CC)-1) * var(2*ShortBar_CC)
rsq=1-SSresid/SStotal

plot(sqrt(ShortBar_WB),2*ShortBar_CC,'-o','LineWidth',4)
title(['Linear Regression, Short Bar. R-Squared = ' num2str(round(rsq,3))],'FontSize',18)
xlabel('Square root of wheel base','FontSize',14)
ylabel('Oscillations per minute','FontSize',14)
hold on
grid on
plot(sqrt(ShortBar_WB),y_est,'r--')
legend({'Oscillation Rate','Linear Best Fit'},'Location','northeast','FontSize',18)
hold off



subplot(323)  % Linear regression. Wood Stick
              % -----------------------------
c=polyfit(sqrt(WoodStick_WB),2*WoodStick_CC,1)
y_est=polyval(c,sqrt(WoodStick_WB))

y_residual = 2*WoodStick_CC-y_est;
SSresid = sum(y_residual.^2)
SStotal = (length(2*WoodStick_CC)-1) * var(2*WoodStick_CC)
rsq=1-SSresid/SStotal

plot(sqrt(WoodStick_WB),2*WoodStick_CC,'-o','LineWidth',4)
title(['Linear regression, Wood Stick. R-Squared = ' num2str(round(rsq,3))],'FontSize',18)
xlabel('Square root of wheel base','FontSize',14)
ylabel('Oscillations per minute','FontSize',14)
hold on
grid on
plot(sqrt(WoodStick_WB),y_est,'r--')
legend({'Oscillation Rate','Linear Best Fit'},'Location','northeast','FontSize',18)
hold off



subplot(324)  % Linear regression. Thin Bar
              % ---------------------------
c=polyfit(sqrt(ThinBar_WB),2*ThinBar_CC,1)
y_est=polyval(c,sqrt(ThinBar_WB))

y_residual = 2*ThinBar_CC-y_est;
SSresid = sum(y_residual.^2)
SStotal = (length(2*ThinBar_CC)-1) * var(2*ThinBar_CC)
rsq=1-SSresid/SStotal

plot(sqrt(ThinBar_WB),2*ThinBar_CC,'-o','LineWidth',4)
title(['Linear Regression, Thin Bar. R-Squared = ' num2str(round(rsq,3))],'FontSize',18)
xlabel('Square root of wheel base','FontSize',14)
ylabel('Oscillations per minute','FontSize',14)
hold on
grid on
plot(sqrt(ThinBar_WB),y_est,'r--')
legend({'Oscillation Rate','Linear Best Fit'},'Location','northeast','FontSize',18)
hold off



subplot(325)  % Linear regression. Plastic Rod
              % ------------------------------
% Remove last three data points
PlasticRod_WB = [ 2.6  3  4  5  6  7  8  9]
PlasticRod_CC = [57   50 43 37 36 34 31 29]

c=polyfit(sqrt(PlasticRod_WB),2*PlasticRod_CC,1)
y_est=polyval(c,sqrt(PlasticRod_WB))

y_residual = 2*PlasticRod_CC-y_est;
SSresid = sum(y_residual.^2)
SStotal = (length(2*PlasticRod_CC)-1) * var(2*PlasticRod_CC)
rsq=1-SSresid/SStotal

plot(sqrt(PlasticRod_WB),2*PlasticRod_CC,'-o','LineWidth',4)
title(['Linear Regression, Plastic Rod. R-Squared = ' num2str(round(rsq,3))],'FontSize',18)
xlabel('Square root of wheel base','FontSize',14)
ylabel('Oscillations per minute','FontSize',14)
hold on
grid on
plot(sqrt(PlasticRod_WB),y_est,'r--')
legend({'Oscillation Rate','Linear Best Fit'},'Location','northeast','FontSize',18)
hold off