clc;
clear;
E=2.1*10.^6;
LowerChordSegment=input('Enter The Number of Segment for Lower Chord: ');
ElementLenghtofLowerChord=input('Enter the Element Lenght of Lower Chord: ');
TrussHight=input('Enter The Hight of Truss: ');
ElementArea=input('Enter The Area of Elements: ');

SinTeta=TrussHight/sqrt((ElementLenghtofLowerChord/2).^2+TrussHight.^2);
CosTeta=(ElementLenghtofLowerChord/2)/sqrt((ElementLenghtofLowerChord/2).^2+TrussHight.^2);

%% In this Section Cordinates of Each Node will be Created
NodeCordinate=zeros(LowerChordSegment+1+(LowerChordSegment-1)+1,3);
x=0;
for i=1:LowerChordSegment+1
    NodeCordinate(i,1)=i;
    NodeCordinate(i,2)=x;
    NodeCordinate(i,3)=0;
    x=x+ElementLenghtofLowerChord;
end
x=ElementLenghtofLowerChord/2;
for i=LowerChordSegment+2:(LowerChordSegment+1+(LowerChordSegment-1)+1)
    NodeCordinate(i,1)=i;
    NodeCordinate(i,2)=x;
    NodeCordinate(i,3)=TrussHight;
    x=x+ElementLenghtofLowerChord;
end
NodeCordinate;

%% In This Section Conectivity Between Each Node Will Be Created
%[m,n]=size(NodeCordinate);
NodeConectivity=zeros(LowerChordSegment+(LowerChordSegment-1)+(LowerChordSegment*2),10);
j=1;
k=1;
for i=1:LowerChordSegment+(LowerChordSegment-1)
    if i==LowerChordSegment+1
        j=j+1;
        NodeConectivity(k,1)=NodeCordinate(j,1);
        NodeConectivity(k,2)=NodeCordinate(j+1,1);
        NodeConectivity(k,3)=NodeCordinate(j,2);
        NodeConectivity(k,4)=NodeCordinate(j,3);
        NodeConectivity(k,5)=NodeCordinate(j+1,2);
        NodeConectivity(k,6)=NodeCordinate(j+1,3);
    else
        NodeConectivity(k,1)=NodeCordinate(j,1);
        NodeConectivity(k,2)=NodeCordinate(j+1,1);
        NodeConectivity(k,3)=NodeCordinate(j,2);
        NodeConectivity(k,4)=NodeCordinate(j,3);
        NodeConectivity(k,5)=NodeCordinate(j+1,2);
        NodeConectivity(k,6)=NodeCordinate(j+1,3);
    end
    j=j+1;
    k=k+1;
end

k=1;
j=LowerChordSegment+2;
for i=LowerChordSegment+(LowerChordSegment-1)+1:LowerChordSegment+(LowerChordSegment-1)+LowerChordSegment
     NodeConectivity(i,1)=NodeCordinate(k,1);
     NodeConectivity(i,2)=NodeCordinate(j,1);
     NodeConectivity(i,3)=NodeCordinate(k,2); 
     NodeConectivity(i,4)=NodeCordinate(k,3);
     NodeConectivity(i,5)=NodeCordinate(j,2);
     NodeConectivity(i,6)=NodeCordinate(j,3);
     k=k+1;
     j=j+1;
end


k=LowerChordSegment+2;
j=2;
for i=LowerChordSegment+(LowerChordSegment-1)+1+LowerChordSegment:LowerChordSegment+LowerChordSegment+(LowerChordSegment-1)+LowerChordSegment
     NodeConectivity(i,1)=NodeCordinate(k,1);
     NodeConectivity(i,2)=NodeCordinate(j,1);
     NodeConectivity(i,3)=NodeCordinate(k,2); 
     NodeConectivity(i,4)=NodeCordinate(k,3);
     NodeConectivity(i,5)=NodeCordinate(j,2);
     NodeConectivity(i,6)=NodeCordinate(j,3);
     k=k+1;
     j=j+1;
end

%% In This Section The Shape of truss Is Bieng Painted
[m,n]=size(NodeConectivity);
[p,q]=size(NodeCordinate);
for i=1:m
    plot([NodeConectivity(i,3) NodeConectivity(i,5)],[NodeConectivity(i,4) NodeConectivity(i,6)],'-o');
    hold on
end
axis([-200,(LowerChordSegment+1)*(ElementLenghtofLowerChord),-(TrussHight+100),(TrussHight+100)]);

for i=1:p
    a=int2str(NodeCordinate(i,1));
    %b=int2str(NodeConectivity(i,2));
    text(NodeCordinate(i,2),NodeCordinate(i,3)-20,a);
end

plot([NodeCordinate(1,2) NodeCordinate(1,2)],[-ElementLenghtofLowerChord/10 ElementLenghtofLowerChord/10],'r');
plot([NodeCordinate(LowerChordSegment+1,2) NodeCordinate(LowerChordSegment+1,2)],[-ElementLenghtofLowerChord/10 ElementLenghtofLowerChord/10],'r');

%% In This Section Global Stiffnes Matrix is Produced

for i=1:m
    NodeConectivity(i,7)=NodeConectivity(i,1)*2-1;
    NodeConectivity(i,8)=NodeConectivity(i,1)*2;
    
     NodeConectivity(i,9)=NodeConectivity(i,2)*2-1;
     NodeConectivity(i,10)=NodeConectivity(i,2)*2;     
end


[m,n]=size(NodeConectivity);
[p,q]=size(NodeCordinate);
K=zeros(p*2);
for i=1:m
    L=sqrt((NodeConectivity(i,6)-NodeConectivity(i,4)).^2+(NodeConectivity(i,5)-NodeConectivity(i,3)).^2); % L is Lenght of Each Element
    KLocal=(ElementArea*E)/(L)*[1 -1;-1 1]; % KLocal is The Local stiffnes Matrix for Each Element
    S=(NodeConectivity(i,6)-NodeConectivity(i,4))/L; % Sin(Teta)
    C=(NodeConectivity(i,5)-NodeConectivity(i,3))/L; % Cos(Teta)
    T=[C S 0 0;0 0 C S]; % T is Matrix of Transformayion From Local Cordinates to Global One
    KGlobal=T'*KLocal*T; % KGlobal is The Global Stiffnes Matrix for Each Element
    for j=1:4
        for k=1:4
            K(NodeConectivity(i,j+6),NodeConectivity(i,k+6))=KGlobal(j,k)+ K(NodeConectivity(i,j+6),NodeConectivity(i,k+6)); % Global Stiffnes Matrix Is Assembled
        end
    end
end

%% In This Section External Force is Produced

% -- Forces Will be implied to Lower Chord of Truss --%
% -- Every Force will also be Implied to Nodes -- %

F=zeros(p*2,1);
for i=4:2:(LowerChordSegment+1)*2-2
    F(i,1)=input('Enter The Force in Y Direction: ');
end
F;


%% In This Section The The EQ KX=F is Going to be Solved

% -- To Imply the Boundery Counditions, Rows and Colums of K Matrix
% -- Representing Restrained Position should be Replaced by Zero -- %

% -- Diagonal Elements of K Matrix Related to Boundery Condition Should --
% -- Replaced By a Logical Number to Avoid a Singular K Matrix --%
sum=0;
for i=1:p*2
    sum=K(i,i)+sum;
end
sum=sum/(p*2);


K(1,:)=0;
K(:,1)=0;
K(2,:)=0;
K(:,2)=0;
K((LowerChordSegment+1)*2-1,:)=0;
K(:,(LowerChordSegment+1)*2-1)=0;
K((LowerChordSegment+1)*2,:)=0;
K(:,(LowerChordSegment+1)*2)=0;    

for i=1:p*2
    K(1,1)=sum;
    K(2,2)=sum;
    K((LowerChordSegment+1)*2-1,(LowerChordSegment+1)*2-1)=sum;
    K((LowerChordSegment+1)*2-1,(LowerChordSegment+1)*2-1)=sum;
    K((LowerChordSegment+1)*2,(LowerChordSegment+1)*2)=sum;
    K((LowerChordSegment+1)*2,(LowerChordSegment+1)*2)=sum;
end

x=inv(K)*F;

% -- x Matrix is being refined to XNew Matrix to determine --
% -- Every Node Deformation more Easily -- %

XNew=zeros(LowerChordSegment+1+(LowerChordSegment-1)+1,2);
j=1;
for i=1:2:p*2
    XNew(j,1)=x(i,1);
    j=j+1;
end
j=1;
for i=2:2:p*2
    XNew(j,2)=x(i,1);
    j=j+1;
end

%% In This Section The deformed Shape Of Truss Is bieng Painted

NodeCordinateFinal=zeros(LowerChordSegment+1+(LowerChordSegment-1)+1,3);
[m,n]=size(NodeCordinate);
for i=1:m
    NodeCordinate(i,2)=NodeCordinate(i,2)+XNew(i,1)*50;% Deformation is Scaled by 50 Times
    NodeCordinate(i,3)=NodeCordinate(i,3)+XNew(i,2)*50;% Deformation is Scaled by 50 Times
end
NodeCordinate;
NodeConectivity;
[m,n]=size(NodeConectivity);
for i=1:m
    k=1;
    for j=1:2:3
        NodeConectivity(i,j+2)=NodeCordinate(NodeConectivity(i,k),2);
        NodeConectivity(i,j+3)=NodeCordinate(NodeConectivity(i,k),3);
        k=k+1;
    end
end

figure(2);
[m,n]=size(NodeConectivity);
for i=1:m
    plot([NodeConectivity(i,3) NodeConectivity(i,5)],[NodeConectivity(i,4) NodeConectivity(i,6)],'-o');
    hold on
end

plot([NodeCordinate(1,2) NodeCordinate(1,2)],[-ElementLenghtofLowerChord/10 ElementLenghtofLowerChord/10],'r');
plot([NodeCordinate(LowerChordSegment+1,2) NodeCordinate(LowerChordSegment+1,2)],[-ElementLenghtofLowerChord/10 ElementLenghtofLowerChord/10],'r');
%plot([-200 (LowerChordSegment+1)*(ElementLenghtofLowerChord)],[0 0]);
axis([-200,(LowerChordSegment+1)*(ElementLenghtofLowerChord),-(TrussHight+100),(TrussHight+100)]);

%% In This Section Global Deformation for Each Element is derived  and Saved to a New Matrix Named XGlobalForElement

 [m,n]=size(NodeConectivity);
[p,q]=size(NodeCordinate);
XGlobalForElement=zeros(m,6);
for i=1:m
    XGlobalForElement(i,1)=NodeConectivity(i,1);
    XGlobalForElement(i,2)=NodeConectivity(i,2);
    XGlobalForElement(i,3)=XNew(NodeConectivity(i,1),1);
    XGlobalForElement(i,4)=XNew(NodeConectivity(i,1),2);
    XGlobalForElement(i,5)=XNew(NodeConectivity(i,2),1);
    XGlobalForElement(i,6)=XNew(NodeConectivity(i,2),2);
end
XGlobalForElement;

%% In This Section Local Deformation is Derived and Subsequently Internal Forces And Stresses is bieng derived

XLocal=zeros(m,3);
FLocal=zeros(2,m);
Stress=zeros(1,m);
for i=1:m
    L=sqrt((NodeConectivity(i,6)-NodeConectivity(i,4)).^2+(NodeConectivity(i,5)-NodeConectivity(i,3)).^2); % L is Lenght of Each Element
    KLocal=(ElementArea*E)/(L)*[1 -1;-1 1]; % KLocal is The Local stiffnes Matrix for Each Element
    KLocalZigma=(E)/(L)*[-1 1]; % KLocalZigma is for Detemining Stress in each Element
    S=(NodeConectivity(i,6)-NodeConectivity(i,4))/L; % Sin(Teta)
    C=(NodeConectivity(i,5)-NodeConectivity(i,3))/L; % Cos(Teta)
    T=[C S 0 0;0 0 C S]; % T is Matrix of Transformayion From Local Cordinates to Global One
    XGlobal=[XGlobalForElement(i,3) XGlobalForElement(i,4) XGlobalForElement(i,5) XGlobalForElement(i,6)];
    XL=T*XGlobal';
    XLocal(i,1)=XGlobalForElement(i,1);
    XLocal(i,2)=XGlobalForElement(i,2);
     XLocal(i,3)=XL(1,1);
     XLocal(i,4)=XL(2,1);
     FL=KLocal*XL;
     FLocal(1,i)=FL(1,1);
     FLocal(2,i)=FL(2,1);
     ZL= KLocalZigma*XL;
     Stress(1,i)=ZL(1,1);
end

%% In This Section a Specified Colour is being Matched to Every Element

sum=0;
for i=1:m
    sum=sum+abs(Stress(1,i));
end
sum=sum/m;

figure(3);
[m,n]=size(NodeConectivity);
for i=1:m
    if abs(Stress(1,i))<=sum
         if abs(Stress(1,i))<=(sum/2)
            plot([NodeConectivity(i,3) NodeConectivity(i,5)],[NodeConectivity(i,4) NodeConectivity(i,6)],'g'); % "Green" is Applied to those Elements stressed less than sum and higher than sum/2
             hold on
             axis([-200,(LowerChordSegment+1)*(ElementLenghtofLowerChord),-(TrussHight+100),(TrussHight+100)]);
             %pause(1);
         else
           plot([NodeConectivity(i,3) NodeConectivity(i,5)],[NodeConectivity(i,4) NodeConectivity(i,6)],'y');  % "Yellow" is Applied to those Elements stressed less than sum/2
             hold on
             axis([-200,(LowerChordSegment+1)*(ElementLenghtofLowerChord),-(TrussHight+100),(TrussHight+100)]);
            % pause(1);
         end
    end
     if abs(Stress(1,i))>sum
         if abs(Stress(1,i))>((sum*1.5))
             plot([NodeConectivity(i,3) NodeConectivity(i,5)],[NodeConectivity(i,4) NodeConectivity(i,6)],'r'); % "Red" is Applied to those Elements stressed higher than 3sum/4
             hold on
             axis([-200,(LowerChordSegment+1)*(ElementLenghtofLowerChord),-(TrussHight+100),(TrussHight+100)]);
             %pause(1);
         else
            plot([NodeConectivity(i,3) NodeConectivity(i,5)],[NodeConectivity(i,4) NodeConectivity(i,6)],'m'); % "M" is Applied to those Elements stressed less than 3sum/4 and higher than sum
             hold on
             axis([-200,(LowerChordSegment+1)*(ElementLenghtofLowerChord),-(TrussHight+100),(TrussHight+100)]);
             % pause(1);
         end
     end
end
plot([NodeCordinate(1,2) NodeCordinate(1,2)],[-ElementLenghtofLowerChord/10 ElementLenghtofLowerChord/10],'k');
plot([NodeCordinate(LowerChordSegment+1,2) NodeCordinate(LowerChordSegment+1,2)],[-ElementLenghtofLowerChord/10 ElementLenghtofLowerChord/10],'k');
axis([-200,(LowerChordSegment+1)*(ElementLenghtofLowerChord),-(TrussHight+100),(TrussHight+100)]);

disp(' ');
disp('~ Press x to see The Defornation matrix: ');
disp('~ Press XNew to see The Deformation for Each Node:  ');
disp(' ');
disp('  Notice -- The Number of XNew Colum Shows The Number associated with Each node -- Notice');
disp(' ');
disp('~ Write FLocal to see The Internal forces: ');
disp(' ');
disp('  Notice -- The Number of FLocal Colum Shows The Number associated with Each node -- Notice');
disp(' ');
disp('~ Write Stress to See The internal Stress: ');