clc
clear
E=input('Elastic modulus=');
v=input('poison ratio=');
d=input('Diameter of cross section area=');
A=pi*(d^2)/4;
noelem=input('input number of element :');     %Define number of elements
nonode=input('input number of nodes :');       %Define number of nodes
sysdof=nonode*2;                               %Define DOF of truss
K=zeros(sysdof);                               %Define stiffnes matrix
frc=zeros(sysdof,1);                           %Define force matrix
disp=zeros(sysdof,1);                          %Define displacement matrix
for i=1:nonode                                 %Define coordinates of nodes
    coord(i,1)=input(sprintf('input X coordinate of node %d=  ',i));
    coord(i,2)=input(sprintf('input Y coordinate of node %d=  ',i));
end
for j=1:noelem
    node1=input(sprintf('input 1st. node of %d element: ',j));     %Define First node of each element
    node2=input(sprintf('input 2nd. node of %d element: ',j));     %Define Sedonde node of each element
    nconnct(j,1)=node1;
    nconnct(j,2)=node2;
    x1=coord(node1,1);
    y1=coord(node1,2);
    x2=coord(node2,1);
    y2=coord(node2,2);
    L=sqrt(((x2-x1)^2)+(y2-y1)^2);               %Calculate lengh of element
    Lengh(j)=L;
    deg=atan((y2-y1)/(x2-x1));                  %Calculate degree of element
    tetha(j)=deg;
    c=cos(deg)/L;
    s=sin(deg)/L;
    Ke=A*E;
    k11=A*E*[c^2 s*c;
        s*c s^2];
    k22=k11;
    k12=A*E*[-(c^2) -s*c;
        -s*c -(s^2)];
    k21=k12;
    if node1<node2
        n1=node1;
        n2=node2;
    else
        n1=node2;
        n2=node1;
    end
    kk11=K(((2*n1)-1):(2*n1),((2*n1)-1):(2*n1));    %input k11 of element matrix in system matrix
    KKe11=k11+kk11;
    K(((2*n1)-1):(2*n1),((2*n1)-1):(2*n1))=KKe11;
    kk12=K(((2*n1)-1):(2*n1),((2*n2)-1):(2*n2));    %input k12 of element matrix in system matrix
    KKe12=k12+kk12;
    K(((2*n1)-1):(2*n1),((2*n2)-1):(2*n2))=KKe12;
    kk21=K(((2*n2)-1):(2*n2),((2*n1)-1):(2*n1));    %input k21 of element matrix in system matrix
    KKe21=k21+kk21;
    K(((2*n2)-1):(2*n2),((2*n1)-1):(2*n1))=KKe21;
    kke22=K(((2*n2)-1):(2*n2),((2*n2)-1):(2*n2));   %input k22 of element matrix in system matrix
    KKe22=k22+kke22;
    K(((2*n2)-1):(2*n2),((2*n2)-1):(2*n2))=KKe22;
end
% Apply constraints
cons=input('input number of nodes that constraints:  ');
for c=1:cons
    nconsno=input('which node is costraint?  ');
    dir=input('how many direction of node was constraint?  ');
    if dir==2;
        disp((2*nconsno)-1,1)=0;
        disp((2*nconsno),1)=0;
    elseif dir==1
        x=1;
        y=2;
        cdir=input('which direction is constraint? x or y?  ');
        if cdir==1
            disp((2*nconsno)-1,1)=0;
        elseif cdir==2
            disp((2*nconsno),1)=0;
        end
    end
end
% Apply loads
lds=input('how many node is forced?  ');
for l=lds
    nfrc=input('which node is forced?  ');
    fdir=input('how many direction of node was forced?  ');
    if fdir==2
        frc((2*nfrc)-1,1)=input('input X component of force');
        frc((2*nfrc),1)=input('input Y component of force');
    elseif dir==1
        x=1;
        y=2;
        ffdir=input('which direction is forced? x or y?  ');
        if ffdir==1
            frc((2*nfrc)-1,1)=input('input X component of force');
        elseif ffdir==2
            frc((2*nfrc),1)=input('input Y component of force');
        end
    end
end