% beams.m
% brody dylan johnson
% 06.26.21

  clear; close all;
  
% data (cm)
  x = [0:5:85];
  un = [75.6 75.5 75.4 75.3 75.1 75.0 74.8 74.6 74.5 74.3 74.1 ...
        73.9 73.8 73.5 73.3 73.0 72.8 72.7];
  vn = [75.6 75.3 75.0 74.5 73.9 73.3 72.4 71.5 70.5 69.5 68.4 ...
        67.3 66.0 64.8 63.4 62.1 61.1 59.6];

% Old Data
%  x = [0 10 20 30 40 50 60 70 80 85];
%  un = [75.2475 75.2475 74.93 74.6125 74.295 73.9775 73.66 ...
%        73.18375 72.7075 72.39];
%  vn = [75.2475 74.9300 73.6600 73.3425 70.1675 67.9450 ...
%        65.4050 62.7063 59.6900 58.4200];
    
% observed deflections
  w1ob = un(1)-un;
  w2ob = un-vn;
    
% physical constants
  g = 9.81;       % N/kg
  rho = 7.8e-3;   % kg/cm^3
  E = 2e7;        % N/cm^2
  b = 0.5*2.54;   % cm
  h = 0.125*2.54; % cm
  I = b*h^3/12;   % cm^4
  ell = 85;       % cm
  W0 = rho*b*h*g; % N/cm
  W1 = 0.5*g;     % N
  a = 80;         % cm
  
% theoretical deflections
  w13dot = -W0*ell/(E*I);
  w12dot = W0*ell^2/(2*E*I);
  w1 = w13dot*x.^3/6 + w12dot*x.^2/2 ...
     + W0*x.^4/(24*E*I);
  w23dot = -W1/(E*I);
  w22dot = W1*a/(E*I);
  w2 = w23dot*x.^3/6 + w22dot*x.^2/2 ...
     + W1*(x-a).^3.*(x>a)/(6*E*I);
 
% graphs
  figure;
  plot(x,w1ob,'sk',x,w1,'-b');
  legend('Observed','Predicted','Location','southeast');
  ylabel('Deflection (cm)');
  xlabel('x (cm)');
  axis([0 85 0 3]);
  set(gca,'FontSize',16)

  figure;
  plot(x,w2ob,'sk',x,w2,'-b');
  legend('Observed','Predicted','Location','southeast');
  ylabel('Deflection (cm)');
  xlabel('x (cm)');
  axis([0 85 0 15]);
  set(gca,'FontSize',16)