%
disp(' ');
disp(' plane_truss.m   ver 1.1  January 9, 2009 ');
disp(' by Tom Irvine  Email: tomirvine@aol.com ');
disp(' ');
disp(' This program calculates the natural frequencies of a plane truss. ');
disp(' Assume same mass density and elastic modulus for each element. ');
%
%  The input file must be an external txt file with the following format:
%
%   number of elements
%   number of nodes
%   element number     first node    second node   area    length    angle(deg)  (repeat for each element)    
%   number of constraints
%   constrained dof   
%
clear num_elem;
clear elem;
clear nc;
clear nnc;
%
disp(' Enter input filename ');
%
[filename, pathname] = uigetfile('*.*');
filename = fullfile(pathname, filename);
fid = fopen(filename,'r');
num_elem = fscanf(fid,'%d',[1 1]);
%
iflag=0;
if(num_elem >=1 )
    out1=sprintf('\n %d elements   \n',num_elem);
    disp(out1);    
else
    num_elem 
    disp(' data input error ');
    iflag=1;
end
%
if(iflag==0)
%
% num_nodes = fscanf(fid,'%d',[1 1]);
%
for(i=1:num_elem)
    elem=fscanf(fid,'%d %d %d %g %g %g ',[1 6])
    e_n1(i)     =elem(1,2);
    e_n2(i)     =elem(1,3);
    e_area(i)   =elem(1,4);
    e_length(i) =elem(1,5);
    e_theta(i)  =elem(1,6);
%
    nz(2*i-1) = e_n1(i);
    nz(2*i)   = e_n2(i);   
%
end
%
zzz=zeros(1000);
for(i=1:2*num_elem)
    zzz(nz(i))=zzz(nz(i))+1;
end
num_nodes=0;
for(i=1:1000)
    if(zzz(i)>=1)
        num_nodes=num_nodes+1;
    end
end
%
out1=sprintf('\n %d nodes  \n',num_nodes);
disp(out1);
%
nc = fscanf(fid,'%d',[1 1]);   % number of constraints
%
out1=sprintf('\n %d  constrained dof   \n',nc);
disp(out1);
%
for(i=1:nc)
    nnc(i)=0;
end    
for(i=1:nc)
    nnc(i)=fscanf(fid,'%d',[1 1]);
end
%
nnc
%
ndof=num_nodes*2;  % number of degrees of freedom, unconstrained
%
MG=zeros(ndof,ndof);
KG=zeros(ndof,ndof);
%
m=zeros(4,4);
m(1,1)=2;
m(2,2)=2;
m(3,3)=2;
m(4,4)=2;
m(1,3)=1;
m(2,4)=1;
m(3,1)=1;
m(4,2)=1;
%
k=zeros(4,4);
k(1,1)=1;
k(3,3)=1;
k(1,3)=-1;
k(3,1)=-1;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
for(i=1:num_elem)
%
    n1  = e_n1(i);
    n2  = e_n2(i);
    A   =e_area(i);
    L   =e_length(i);
    theta   =e_theta(i);
%
    theta=  theta*pi/180.;    
%
    n(1)=2*n1-1;  
    n(2)=2*n1;
    n(3)=2*n2-1;
    n(4)=2*n2;
%
    [MG,KG]=coordinate_transform(n,L,A,theta,m,k,MG,KG);
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
MG;
KG;
%
MG=MG/6;
%
mmm=ndof-nc;
%
out1=sprintf('\n %d  dof after constraints applied   \n',mmm);
disp(out1);
%
MGC=zeros(mmm,mmm);
KGC=zeros(mmm,mmm);
%
ir=0;
for(i=1:ndof)          % row
    ic=1;
    row_flag=0;
    for(nv=1:nc)
        if(i==nnc(nv))
            row_flag=1;
        end
    end
%
    if(row_flag==0)
        ir=ir+1;
%
        for(j=1:ndof)      % column
            column_flag=0;
            for(nv=1:nc)
                if(j==nnc(nv))
                    column_flag=1;
                end
            end                
%
            if(column_flag==0)
%                
                MGC(ir,ic)=MG(i,j);
                KGC(ir,ic)=KG(i,j);
                ic=ic+1;
%  
            end    
        end        
%
    end
% 
end
%
disp(' Constrained mass matrix ');
MGC
disp(' Constrained stiffness matrix ');
KGC
%
clear ModeShapes;
clear EigenValues;
[ModeShapes,EigenValues]=eig(KGC,MGC);
EigenValues;
%
disp(' ');
disp(' Natural Frequencies (rad/sec) ');
disp(' ');
%
for(i=1:mmm)
    omega(i)=sqrt(EigenValues(i,i)); 
    out1=sprintf('  %8.4g  (1/L)*sqrt(E/rho) ',omega(i));
    disp(out1);
end
%
MST=ModeShapes';
temp=zeros(mmm,mmm);
temp=MGC*ModeShapes;
QTMQ=MST*temp;
%   
for(i=1:mmm)
    scale(i)=1./sqrt(QTMQ(i,i));
    ModeShapes(:,i) = ModeShapes(:,i)*scale(i);    
end
%
MST=ModeShapes';
%
disp(' ');
disp('  Modes Shapes (column format)');
ModeShapes
%
end
fclose(fid);