% compute the first three natural frequencies of a cable with
% flexible supports and end masses
% INPUT DATA YOU NEED TO PROVIDE
% delclare the beam dimensions and properties
L = 100;            % cable length in meters
A = 0.001;          % cable cross section area m^2
rho = 3*2700*A;     % cable mass per unit length in kg/m
mass_cable = rho*L; % total cable mass
T = 5000;           % cable tension in Newton
%
% Account for the tip mass;
%
mu_L = 100;  % the ratio of the concentrated tip mass at x=0 to mass_cable
mu_0 = 200;  % the ratio of the concentrated tip mass at x=L to mass_cable
k_0 = 10000; % left-end stiffness relative to tension T;
k_L = 10000; % right-end stiffness relative to tension T;
%
% ENDS YOUR INPUT DATA. THE REST DOES IT FOR YOU!
%
% beta_try approximate roots obtained from determinant search

for n=1: 5,
%  k = 10^n; % the support spring constants relative to the tension T, k_L = k_0 = k;
%
% compute the first natural frequency
%
beta_1 =beta_try(n);  % the starting value from completely free-free case;
x = beta_1;
norm = 1;
while (norm>1.e-06)
sx = sin(x); cx = cos(x); tx = tan(x);

det = ((k_0)  - mu_0*x^2)*((k_L)  - mu_L*x^2)*sin(x);
det = det + (((k_0)  - mu_0*x^2) +((k_L)  - mu_L*x^2))*x*cos(x);
det = det-x^2*sin(x);
norm = abs(det);

slope =        - (2*mu_0*x)*((k_L)  - mu_L*x^2)*sin(x);
slope = slope  - (2*mu_L*x)*((k_0)  - mu_0*x^2)*sin(x);
slope = slope  + ((k_0)  - mu_0*x^2)*((k_L)  - mu_L*x^2)*cos(x);
slope = slope  + (((k_0)  - mu_0*x^2) +((k_L)  - mu_L*x^2))*cos(x);
slope = slope  - (((k_0)  - mu_0*x^2) +((k_L)  - mu_L*x^2))*x*sin(x);
slope = slope  - 2*(   mu_0  +  mu_L )*x^2*cos(x);
slope = slope  - 2*x*sin(x) - x^2*cos(x);

if (abs(slope)>1.e-08)
xnew = x - det/slope;
else
xnew = x;
end;
x= xnew;
end;
% compute the first natural frequency in Hertz
beta_1 = x;
disp( ' ' );
%disp( [' Frequency factor, beta_1 =' , num2str(beta_1)]);
omega = (beta_1/L)* sqrt(T/rho);
freq = omega/(2*pi); %in Hertz
%disp( ' ' );
%disp([' Fundamental natural frequency in Hertz = ' , num2str(freq)]);
%
% compute the mode shape
%
%beta = beta_1/L;
x_coord=zeros(1,0);
y_coord =zeros(1,0);

nom = (k_L - mu_L*beta_1^2)*cos(beta_1) - beta_1*sin(beta_1)
denom = (k_L - mu_L*beta_1^2)*sin(beta_1) + beta_1*cos(beta_1)
coef = nom/denom;

for i=0:100,
span = i/100;
arg =beta_1*span;
sx = sin(arg); cx = cos(arg);
W = -coef*sx + cx;
x_coord = [x_coord span];
y_coord =[y_coord W];
end;

%normalize y_coord
y_coord = y_coord/(max(abs(y_coord)));
figure(n);
plot(x_coord, y_coord);
xlabel('cable span');
ylabel('mode shape ');
legend(['frequency in Hertz = ', num2str(freq), ' beta*L = ', num2str(beta_1)]);
title('Mode shape of the first natural frequency of cable with flexible supports');

end;  %end the outer loop