clc;
clear;
E=2.1*10.^6;
v=0.3;
NumOfSpan=input('Enter The Number of Span: ');
NumOfStorey=input('Enter The Numbert of Storey: ');
HieghtOfEachStorey=input('Enter the Hieght Of Each Storey: ');
A=1;
I=10;
ElementLenght=HieghtOfEachStorey/4;
pause(0.5);
disp(' ');
disp(' **    **');
disp(' ');
disp('Area and Inertia of Each Element is the Same');
disp('Area and Inertia Can only Be Changed in the Code');
disp(' ');
disp(' **    **');
disp(' ');
pause(1);
%% In This Section Lenght Of Each Span Will Be Apllied
SpanLenght=zeros(1,NumOfSpan);
for i=1:NumOfSpan
    SpanLenght(1,i)=input('Enter The Lenght Of Each Span From Left To Right: ');
end
disp(' ');
%% In This Section Coordinates Of Each Element Will Be Applied
NodesCoordinate=zeros(((NumOfStorey)*4)*(NumOfSpan+1)+(NumOfSpan+1)+(NumOfSpan)*3*NumOfStorey,3);
size(NodesCoordinate);
NodesX1=zeros(1,NumOfSpan+1);
NodesX2=zeros(1,4*NumOfSpan);
for i=2:NumOfSpan+1
    if i==2
        NodesX1(1,i)=SpanLenght(1,i-1);
    else
        
        NodesX1(1,i)=SpanLenght(1,i-2)+SpanLenght(1,i-1)+NodesX1(1,i-2);
    end
end 
x=1;
    for i=1:NumOfSpan
        k=1;
        for j=x:x+3
            if i==1
                if j==1
                   NodesX2(1,j)=SpanLenght(1,i)/5;
                else
                NodesX2(1,j)=SpanLenght(1,i)/5*j;
                k=k+1;
                end
                else
            NodesX2(1,j)=NodesX1(1,i)+SpanLenght(1,i)/5*k;
            k=k+1;
            end
        end
        x=x+4; 
    end
x=1;
for i=0:(NumOfStorey)*4
    k=1;
    for j=x:x+NumOfSpan
        NodesCoordinate(j,1)=j;
        NodesCoordinate(j,2)=NodesX1(1,k);
        NodesCoordinate(j,3)=i*ElementLenght;
        k=k+1;
    end
    x=x+NumOfSpan+1;
end

[m,q]=size(NodesX2);
for i=1:NumOfStorey
    k=1;
    for j=x:x+q-1
        NodesCoordinate(j,1)=j;
        NodesCoordinate(j,2)=NodesX2(1,k);
        NodesCoordinate(j,3)=i*HieghtOfEachStorey;
        k=k+1;
    end
    x=x+q;
end
disp('Coordinate of Each Element Has Been Created');
disp(' ');
pause(0.5);
%% In This Section Conectivity Between Each Node Will Be Created
NodesConectivity=zeros(((NumOfSpan+1)*4)*NumOfStorey+((NumOfSpan*5)*NumOfStorey),6);
for i=1:((NumOfSpan+1)*4*NumOfStorey)
    NodesConectivity(i,1)=i;
    NodesConectivity(i,2)=i+(NumOfSpan+1);
    NodesConectivity(i,3)=NodesCoordinate(i,2);
    NodesConectivity(i,4)=NodesCoordinate(i,3);
    NodesConectivity(i,5)=NodesCoordinate(i+(NumOfSpan+1),2);
    NodesConectivity(i,6)=NodesCoordinate(i+(NumOfSpan+1),3);
end
j=1;
k=((NumOfSpan+1)*4)*NumOfStorey+1;
for i=(((NumOfSpan+1)*4*NumOfStorey)+(NumOfSpan+1))+1:((NumOfSpan+1)*4)*NumOfStorey+((NumOfSpan*4)*NumOfStorey)+(NumOfSpan+1)
    if j==4
        j=1;
    else
    NodesConectivity(k,1)=i;
    NodesConectivity(k,2)=i+1;
    NodesConectivity(k,3)=NodesCoordinate(i,2);
    NodesConectivity(k,4)=NodesCoordinate(i,3);
    NodesConectivity(k,5)=NodesCoordinate(i+1,2);
    NodesConectivity(k,6)=NodesCoordinate(i+1,3);
    j=j+1;
    k=k+1;
    end
end
x=(NumOfSpan+1)*4*NumOfStorey+(NumOfSpan+1)+1;
y=((NumOfSpan+1)*4)+1;
for i=1:NumOfStorey
    for j=1:NumOfSpan
        NodesConectivity(k,1)=y;
        NodesConectivity(k,2)=x;
        NodesConectivity(k,3)=NodesCoordinate(y,2);
        NodesConectivity(k,4)=NodesCoordinate(y,3);
        NodesConectivity(k,5)=NodesCoordinate(x,2);
        NodesConectivity(k,6)=NodesCoordinate(x,3);
        y=y+1;
        x=x+4;
        k=k+1;
    end
    y=y+(NumOfSpan+1)*4-NumOfSpan;
end
x=(NumOfSpan+1)*4+2;
y=(NumOfSpan+1)*4*NumOfStorey+(NumOfSpan+1)+4;
for i=1:NumOfStorey
    for j=1:NumOfSpan
        NodesConectivity(k,1)=y;
        NodesConectivity(k,2)=x;
        NodesConectivity(k,3)=NodesCoordinate(y,2);
        NodesConectivity(k,4)=NodesCoordinate(y,3);
        NodesConectivity(k,5)=NodesCoordinate(x,2);
        NodesConectivity(k,6)=NodesCoordinate(x,3);
        y=y+4;
        x=x+1;
        k=k+1;
    end
    x=x+(NumOfSpan+1)*4-NumOfSpan;
    %y=y+3;
end
disp('NodesConectivity of Each Element Has Been Created');
disp(' ');
pause(0.5);
%% In This Section The Shape of truss Is Bieng Painted
[m,n]=size(NodesConectivity);
[p,q]=size(NodesCoordinate);
for i=1:m
    plot([NodesConectivity(i,3) NodesConectivity(i,5)],[NodesConectivity(i,4) NodesConectivity(i,6)],'-o');
    hold on
end
axis([-200,sum(SpanLenght)+200,-(HieghtOfEachStorey+50),(NumOfStorey)*(HieghtOfEachStorey)+100]);
for i=1:p
    a=int2str(NodesCoordinate(i,1));
    %b=int2str(NodeConectivity(i,2));
    text(NodesCoordinate(i,2),NodesCoordinate(i,3)-20,a);
end
disp('Shape of Frame Has Been Created');
disp(' ');
pause(0.5);
%% In This Section Global Stiffnes Matrix is Produced
for i=1:m
    NodesConectivity(i,7)=NodesConectivity(i,1)*3-2;
    NodesConectivity(i,8)=NodesConectivity(i,1)*3-1;
    NodesConectivity(i,9)=NodesConectivity(i,1)*3;
    
    NodesConectivity(i,10)=NodesConectivity(i,2)*3-2;
    NodesConectivity(i,11)=NodesConectivity(i,2)*3-1; 
    NodesConectivity(i,12)=NodesConectivity(i,2)*3; 
end
[m,n]=size(NodesConectivity);
[p,q]=size(NodesCoordinate);
K=zeros(p*3);
for i=1:m
    L=sqrt((NodesConectivity(i,6)-NodesConectivity(i,4)).^2+(NodesConectivity(i,5)-NodesConectivity(i,3)).^2);
    KLocal=(E)/L*[A 0 0 -A 0 0;0 12*(I/(L)^2) 6*(I/L) 0 -12*(I/(L)^2) 6*(I/L);0 6*(I/L) 4*I 0 -6*(I/L) 2*I;-A 0 0 A 0 0;0 -12*(I/(L)^2) -6*(I/L) 0 12*(I/(L)^2) -6*(I/L);0 6*(I/L) 2*I 0 -6*(I/L) 4*I];
    S=(NodesConectivity(i,6)-NodesConectivity(i,4))/L; % Sin(Teta)
    C=(NodesConectivity(i,5)-NodesConectivity(i,3))/L; % Cos(Teta)
    T=[C S 0 0 0 0;-S C 0 0 0 0;0 0 1 0 0 0;0 0 0 C S 0;0 0 0 -S C 0;0 0 0 0 0 1];
    KGlobal=T'*KLocal*T; % KGlobal is The Global Stiffnes Matrix for Each Element
    for j=1:6
        for k=1:6
            K(NodesConectivity(i,j+6),NodesConectivity(i,k+6))=KGlobal(j,k)+ K(NodesConectivity(i,j+6),NodesConectivity(i,k+6)); % Global Stiffnes Matrix Is Assembled
        end
    end
end
disp('Global Stifness Matrix Has Been Assembled');
disp(' ');
pause(0.5);
disp(' **     ** ')
pause(0.5);
disp(' ')
%% In This Section External Force is Produced
F=zeros(p*3,1);
y=(NumOfSpan+1)*4+1;
x=(y)*3-2;
for i=1:NumOfStorey
    F(x,1)=input('Enter The Force in Y Direction: ');
    y=y+(NumOfSpan+1)*4;
    x=y*3-2;
end
F;

%% In This Section The The EQ KX=F is Going to be Solved
sum=0;
for i=1:p*2
    sum=K(i,i)+sum;
end
sum=sum/(p*2);
for i=1:(NumOfSpan+1)*3
    K(i,:)=0;
    K(:,i)=0;
end
for i=1:(NumOfSpan+1)*3
    K(i,i)=sum;
end
K;
x=inv(K)*F;
disp(' ');
disp('F=KX is Being solved');
disp(' . ');
pause(0.75);
disp(' . ');
pause(0.75);
disp(' . ');
pause(0.75);
disp(' . ');
pause(0.5);
[m,n]=size(NodesCoordinate);
Xnew=zeros(m,3);
[p,q]=size(x);
p=p/3;
for i=1:p
    Xnew(i,1)=i;
    Xnew(i,2)=x(i*3-2,1);
    Xnew(i,3)=x(i*3-1,1);
end
Xnew;
for i=1:m
    NodesCoordinate(i,2)=NodesCoordinate(i,2)+Xnew(i,2)*5;% Deformation is Scaled by 5 Times
    NodesCoordinate(i,3)=NodesCoordinate(i,3)+Xnew(i,3)*5;% Deformation is Scaled by 5 Times
end
[m,n]=size(NodesConectivity);
for i=1:m
    k=1;
    for j=1:2:3
        NodesConectivity(i,j+2)=NodesCoordinate(NodesConectivity(i,k),2);
        NodesConectivity(i,j+3)=NodesCoordinate(NodesConectivity(i,k),3);
        k=k+1;
    end
end
figure(2);
[m,n]=size(NodesConectivity);
for i=1:m
    plot([NodesConectivity(i,3) NodesConectivity(i,5)],[NodesConectivity(i,4) NodesConectivity(i,6)],'-o');
    hold on
end
%axis([-200,sum(SpanLenght)+200,-(HieghtOfEachStorey+50),(NumOfStorey)*(HieghtOfEachStorey)+100]);
disp('Shape of Deformed Frame has Been Ploted');
disp(' ');
pause(0.5);
%% Determine Internal Force of Each Element

[m,n]=size(NodesConectivity);
[p,q]=size(NodesCoordinate);
F=zeros(6,m);
XLocal=zeros(1,6);
for i=1:m
    XLocal(1,1)=x(NodesConectivity(i,1)*3-2,1);
    XLocal(1,2)=x(NodesConectivity(i,1)*3-1,1);
    XLocal(1,3)=x(NodesConectivity(i,1)*3,1);
    XLocal(1,4)=x(NodesConectivity(i,2)*3-2,1);
    XLocal(1,5)=x(NodesConectivity(i,2)*3-1,1);
    XLocal(1,6)=x(NodesConectivity(i,2)*3,1);
   
    L=sqrt((NodesConectivity(i,6)-NodesConectivity(i,4)).^2+(NodesConectivity(i,5)-NodesConectivity(i,3)).^2);
    KLocal=(E)/L*[A 0 0 -A 0 0;0 12*(I/(L)^2) 6*(I/L) 0 -12*(I/(L)^2) 6*(I/L);0 6*(I/L) 4*I 0 -6*(I/L) 2*I;-A 0 0 A 0 0;0 -12*(I/(L)^2) -6*(I/L) 0 12*(I/(L)^2) -6*(I/L);0 6*(I/L) 2*I 0 -6*(I/L) 4*I];
    S=(NodesConectivity(i,6)-NodesConectivity(i,4))/L; % Sin(Teta)
    C=(NodesConectivity(i,5)-NodesConectivity(i,3))/L; % Cos(Teta)
    T=[C S 0 0 0 0;-S C 0 0 0 0;0 0 1 0 0 0;0 0 0 C S 0;0 0 0 -S C 0;0 0 0 0 0 1];
    FElement=KLocal*T*XLocal';
    F(1,NodesConectivity(i,1))=FElement(1,1);
    F(2,NodesConectivity(i,1))=FElement(2,1);
    F(3,NodesConectivity(i,1))=FElement(3,1);
    F(4,NodesConectivity(i,1))=FElement(4,1);
    F(5,NodesConectivity(i,1))=FElement(5,1);
    F(6,NodesConectivity(i,1))=FElement(6,1);
end
disp('Press F to See The Internal Forces');
disp('Press x to See The Each Nodes Dispplacement');
disp(' ');