clear all;clc
% Lab7:- Eigenmodes of an Euler-Bernoulli Beam

% set up physical parameters
E = 210E9; % elastic modulus(N/m^2)
Izz = 5E-05; % moment of inertia (m^4)
p = 8.3333; % loading, force per length
L = 1.0; % length of the beam (m)
mass = 200; %mass of the beam (kg)
k = 2E06; %spring constant (N/m) 
for step = 5:1:5
nnodes = 2^step + 1

%rho*A*L = mass, where L = 1
rhoAL = mass*ones(nnodes-1,1);

% set up mesh
X1 = linspace(0,L,nnodes);

load femHlib; % load cubic Hermite element library matrices

% form the global load column matrix
MASS = ass1dH(1,rhoAL,[],1,A200H,X1);

%remove the first 2 dirichlet rows
MASSdiri = MASS(3:end,3:end);

% form the global stiffness matrix
K = ass1dH(E*Izz,[],[],-3,A222H,X1);

%Spring matrix
SPRING = zeros(2*nnodes);
SPRING(2*nnodes-1, 2*nnodes-2) = k;
K= K+SPRING;

%remove the first 2 dirichlet rows
Kdiri = K(3:end,3:end);

%solve for the eigenvalues and eigenvectors
[evectors evalues flag] = eigs(inv(MASSdiri)*Kdiri,4,'SM');

OMG_SQRD = diag(evalues);
format long
omg = sqrt(OMG_SQRD);
format short


%look at the first four eigenvalues
for jj=1:1:4

Q = evectors(:,jj);

%separate the slopes from the mode shapes
y = zeros(nnodes,1);
slope = zeros(nnodes,1);

count = 1;
for ii=1:1:nnodes
if ii==1 %re-apply BC at left
y(ii)=0;
slope(ii) = 0;
else
y(ii) = Q(count);
count = count + 1;
slope(ii) = Q(count);
count = count + 1;
end
end

y
slope
%fprintf('\n %2.6f %2.6f ',y,slope)


%plot modes
figure(1) %plot each mode
subplot(4,1,jj)
plot(X1,real(y),'bo')
if jj==1
title({'E-B Beam cantilever beam displacement,nnodes=' num2str(nnodes) ' omega=' num2str(omg(jj))})
else
title({'omega=' num2str(omg(jj))})
end

figure(2) %plot each mode
subplot(4,1,jj)
plot(X1,real(slope),'r^')
if jj==1
title({'E-B Beam cantilever beam slope, nnodes=' num2str(nnodes) 'omega=' num2str(omg(jj))})
else
title({'omega=' num2str(omg(jj))})
end

end %end mode loop

end %end mesh loop
