function wn=vtb6_4(n,E,G,I,rho,A,L,K)%VTB6_4 Comparison of natural frequencies for an Euler-Bernoulli%beam, a beam including rotary inertia, and a beam including %shear deformation.% w = VTB6_4(n,E,G,I,rho,A,L,K) will return the first n natural % frequencies for the three types of beams mentioned above. The input % parameters are:% % E = Young's Modulus% G = Shear Modulus% I = moment of inertia about bending axis% rho = density% A = cross-sectional area% L = length% K = shear coefficient%% The output w is a n x 3 matrix.  The first column is the% frequencies for an E-B beam, the second column corresponds to% a beam with rotary inertia, and the third column is the beam% with shear deformation.alpha=sqrt(E*I/(rho*A));r=sqrt(I/A);for i=1:n    wn(i,1)=sqrt((alpha^2*i^4*pi^4)/L^4);   %E_B beam Note (1)    wn(i,2)=(alpha^2*i^4*pi^4)/L^4;    wn(i,2)=sqrt(wn(i,2)*(1/(1+i^2*pi^2*r^2/L^2)));  %Rotary inertia    wn(i,3)=(alpha^2*i^4*pi^4)/L^4;  %Shear deformation    wn(i,3)=sqrt(wn(i,3)*(1/(1+(i^2*pi^2*r^2/L^2)*(E/(K*G)))));end