%  Computations of beam bending vibrations 
%  using a general characteristic  equation
%
clear all;
%format long

% specified 5 case: free-free,fixed-free, s-s,  fix-fix, fix-s     
k=10^8; % kw and ktehta  are assumed to have excessively large numbers!
k_w1     =[ 0  k  k  k  k];   
k_theta1 =[ 0  k  0  k  k];  
k_w2     =[ 0  0  k  k  k];  
k_theta2 =[ 0  0  0  k  0];

% put them into an array
Kinput =[k_w1; k_theta1; k_w2; k_theta2];

% the bounday condition cases
FreeFree   ='(Free-Free beam)';
FixedFree  ='(Fixed-Free beam)';
SSbeam     ='(Simply supported-Simply supported beam)';
FixedFixed ='(Fixed-Fixed beam)';
FixedS     ='(Fixed-Simply supported beam)';

disp(' ')
disp('Please select the bounday conditions');
disp('free-free beam:   1');
disp('fixed-free beam:  2');
disp('fixed-simply supported beam: 3');
disp('fixed-fixed beam: 4');
disp('fixed-simply supported beam: 5');
bccase =input('boundary condition for this analysis  = ');
disp(' ');
n_modes =input('specify number of modes to be computed (limit to max. 10!) =  ');

%
% load the initial beta*L data
%
load vibdata betaL_initial

% pick the corresponding boundary condition  for title print out
if bccase==1 
	symbol = FreeFree;
elseif bccase==2
	symbol= FixedFree;
elseif  bccase==3
    symbol = SSbeam;
elseif bccase==4
    symbol= FixedFixed;
 elseif  bccase==5
	 symbol = FixedS;
end; 

%  select the  boundary conditions
	k_w1 =Kinput(1,bccase);
	k_theta1 =Kinput(2,bccase);
	k_w2 =Kinput(3,bccase);
	k_theta2 =Kinput(4,bccase);
	
% provide initial roots (betaL) for this case from 
% computed by  DeterminantSearchBeam.m

betaL  = betaL_initial(bccase,:);
%nmodes =  max(size(betaL));

%
% compute the mode and mode shape of the five modes
%
for mode =1: n_modes	

betaL_accurate = betaL(mode);  % the starting value  of beta*L

norm = 1; iter = 0;
%seek for the root for this starting value
while (norm>1.e-07 & iter < 100),
iter = iter +1;
dbetaL  = Newton_Increment ( k_w1, k_w2, k_theta1, k_theta2, betaL_accurate);

betaL_accurate = betaL_accurate + dbetaL;
end;

%
% compute mode shape  for this frequency
%
Cbeam = CmatrixBeamGeneral ( k_w1, k_w2, k_theta1, k_theta2, betaL_accurate);

% it is assumed that the minor Cbeam(3x3) matrix is non-signular
Cminor = Cbeam(1:3, 1:3);
b=-Cbeam(1:3,4);
% mode shape coefficients assuming c_4 =1
coef = Cminor\b;

x_coord=zeros(1,0);
y_coord =zeros(1,0);

for i=0:100,
span = i/100;
arg =betaL_accurate*span;
sx = sin(arg); cx = cos(arg);
shx =sinh(arg); chx = cosh(arg);
W =coef(1)*sx + coef(2)*cx + coef(3)*shx + chx; 
x_coord = [x_coord span];
y_coord =[y_coord W];
end;

%normalize y_coord
y_coord = y_coord/(max(abs(y_coord)));
figure(mode);
plot(x_coord, y_coord);
xlabel('Beam  span');
ylabel('Mode shape ');
%legend(['[k_w1  k_theta1  k_w2  k_theta2] =', num2str(Kinput(bccase,:)),  '  beta*L =  ', num2str(betaL_accurate,7)]);
legend(['Case of boundary conditions = ', num2str(bccase),  '  beta*L =  ', num2str(betaL_accurate,7)]);
title([symbol]);

end;  %end mode=1:nmodes loop