%% Project Information
    % 1.575 Computational Structural Design and Optimization
    % Final Project
    % Fall 2016
    % Nablul Haseeb
    % Version 1.0 (Beta)
    
%% Assumptions
    % Loads act along entire width of beam
    % No node carries more than four struts/ties
    % Beam is either fixed-fixed or pinned-pinned

%% CLearing exisitng outputs/figures
    clear all
    clc
    clf

%% Diagram Control Variables
    Pointload_scale = 1.00;
    Distload_scale = 2.00;
    Arrow_size = 2;
    BearingPL_depth = 1;
    offset_top = 1;
    offset_bot = 1;
    offset_hor = 6;
    increment = 12;
    
%% Geometric/Material Inputs
    %l = input('Please enter the clear length of the beam in ft: ');
    %h = input('Please enter the depth of the beam in inches: ');
    %w = input('Please enter the width of the beam in inches: ');
    %fc = input('Please enter the compressive strength of concrete in psi: ');
    %fy = input('Please enter the yielding strength of steel in ksi: ');
    %l_res = input('Please enter 1 if left end is fixed, else enter 0: ');
    %r_res = input('Please enter 1 if right end is fixed, else enter 0: ');
    %stmtype = input('Please enter STM degree of analysis (1 or 2): ');
    %disp('---------------------------------------------------------------')

%% Input Override (temp.)
    l = 20;
    h = 60;
    w = 24;
    fc = 4000;
    fy = 60;
    l_res = 1;
    r_res = 1;
    stmtype = 2;

%% Loading Inputs
    % User input of number of point/distributed loads. Initial matrix setup
    p_num = input('Please enter the number of point loads acting on the beam: ');
    w_num = input('Please enter the number of distributed loads acting on the beam: ');
    disp('---------------------------------------------------------------')
    temp = zeros(50,3);
    temp_w = zeros(w_num,5);
    temp_p = zeros(p_num,3);
    load_num = 0;
    
    % Populating point loads into relevant matrix
    for i = 1:p_num
        temp(i,1) = input(['Please enter x distance of point load #' num2str(i) ' in ft: ']);
        temp(i,2) = input(['Please enter magnitude of load #' num2str(i) ' in kips: ']);
        temp(i,3) = 5; %input('Please enter the bearing width of load in inches: ');
        temp_p(i,1) = temp(i,1);
        temp_p(i,2) = temp(i,2);
        temp_p(i,3) = temp(i,3);
        load_num = load_num + 1;
        disp('---------------------------------------------------------------')
    end
    
    % Populating distributed loads into relevant matrix
    for i = 1:w_num
        temp_w(i,1) = input(['Please enter starting x distance of distributed load #' num2str(i) ' in ft: ']);
        temp_w(i,2) = input(['Please enter ending x distance of distributed load #' num2str(i) ' in ft: ']);
        temp_w(i,3) = input(['Please enter number of segments to divide distributed load #' num2str(i) ':']);
        temp_w(i,4) = input(['Please enter magnitude of distributed load #' num2str(i) ' in kips/ft: ']);
        temp_w(i,5) = 5; %input('Please enter the bearing width of load in inches: ');
        temp_segment = 1;
        
        for i2 = (load_num + 1):(load_num + temp_w(i,3))
            temp(i2,1) = temp_w(i,1) + ((((temp_w(i,2)) - (temp_w(i,1)))/((temp_w(i,3)) + 1))*temp_segment);
            temp(i2,3) = temp_w(i,5);
            
            if (temp_segment == 1) || (temp_segment == temp_w(i,3))
                temp(i2,2) = 1.5*((temp_w(i,2) - temp_w(i,1))/(temp_w(i,3) + 1))*temp_w(i,4);
            else
                temp(i2,2) = 1.0*((temp_w(i,2) - temp_w(i,1))/(temp_w(i,3) + 1))*temp_w(i,4);
            end
            load_num = load_num + 1;
            temp_segment = temp_segment + 1;
        end
        disp('---------------------------------------------------------------')
    end

    % Populating final "loads" matrix
    load = zeros(load_num,3);
    for i = 1:load_num
        load(i,1) = temp(i,1);
        load(i,2) = temp(i,2);
        load(i,3) = temp(i,3);
    end
    
    % Calculating total number of loads
    load_num = size(load,1);
    
%% Boundary-Reaction Calculations
    FEMi = 0;
    FEMf = 0;
    Ri = 0;
    Rf = 0;
    
    % Calculating end reactions/moments for fixed-fixed condition
    if (l_res == 1) && (r_res == 1)
        for i = 1:p_num
            a = temp_p(i,1);
            b = (l - a);
            P = temp_p(i,2);
            FEMi = FEMi + (P*(b^2)*a)/(l^2);
            FEMf = FEMf + (P*(a^2)*b)/(l^2);
            Ri = Ri + ((P*(b^2)*((3*a) + b))/(l^3));
            Rf = Rf + ((P*(a^2)*((3*b) + a))/(l^3));
        end
        
        for i = 1:w_num
            q = temp_w(i,4);
            d = (temp_w(i,2)) - (temp_w(i,1));
            a = (temp_w(i,1)) + (d/2);
            b = l - a;
            FEMi = FEMi + (((q*d)/(l^2))*((a*(b^2)) + (((a - (2*b))*(d^2))/12)));
            FEMf = FEMf + (((q*d)/(l^2))*((b*(a^2)) + (((b - (2*a))*(d^2))/12)));
            Ri = Ri + (((q*d)/(l^3))*((((2*a) + l)*(b^2)) + (((a - b)/4)*(d^2))));
            Rf = Rf + (((q*d)/(l^3))*((((2*b) + l)*(a^2)) - (((a - b)/4)*(d^2))));
        end
    end
    
    % Calculating end reactions/moments for pinned-pinned condition
    if (l_res == 0) && (r_res == 0)
        for i = 1:p_num
            a = temp_p(i,1);
            b = (l - a);
            P = temp_p(i,2);
            Ri = Ri + ((P*b)/l);
            Rf = Rf + ((P*a)/l);
        end
        
        for i = 1:w_num
            d = (temp_w(i,2)) - (temp_w(i,1));
            a = (temp_w(i,1)) + (d/2);
            b = l - a;
            P = d*temp_w(i,4);
            Ri = Ri + ((P*b)/l);
            Rf = Rf + ((P*a)/l);
        end
    end  
    
%% STM Formulation
    % Calculating number of joints and members. Initial matrix setup
    if stmtype == 1
        j = (2*load_num + 1);
        m = (4*load_num - 1);
    end
    if stmtype == 2
        j = (4*load_num + 3);
        m = (8*load_num + 3);
    end

    jointdata = zeros(j,7);
    memberdata = zeros(m,6);
    
    % Assigning joint/member ID
    for i = 1:j
        jointdata(i,1) = i;
    end
    for i = 1:m
        memberdata(i,1) = i;
        memberdata(i,4) = 1;
        memberdata(i,5) = 1;
    end
    
    maincounter = 1;
    
    if stmtype == 1
        big_temp = zeros(1000000,7);
    end
    if stmtype == 2
        big_temp = zeros(1000000,17);
    end

%     for i1 = 0:increment:(load(2,1) - load(1,1))*12
%     for i2 = 0:increment:(load(3,1) - load(2,1))*12
%     for i3 = 0:increment:(load(4,1) - load(3,1))*12
%     nodestep(1,1) = i1;
%     nodestep(2,1) = i2;
%     nodestep(3,1) = i3;
    
    for it1 = 0:increment:((load(1,1))*12) - offset_hor
    for it2 = 0:increment:(load(2,1) - load(1,1))*12
    for it3 = 0:increment:(load(3,1) - load(2,1))*12
    for it4 = 0:increment:(load(4,1) - load(3,1))*12
    for it5 = 0:increment:((l*12) - offset_hor) - ((load(4,1))*12)
    for i1 = 0:increment:((load(1,1))*12) - offset_hor
    for i2 = 0:increment:(load(2,1) - load(1,1))*12
    for i3 = 0:increment:(load(2,1) - load(1,1))*12
    for i4 = 0:increment:(load(3,1) - load(2,1))*12
    for i5 = 0:increment:(load(3,1) - load(2,1))*12
    for i6 = 0:increment:(load(4,1) - load(3,1))*12
    for i7 = 0:increment:(load(4,1) - load(3,1))*12
    for i8 = 0:increment:((l*12) - offset_hor) - ((load(4,1))*12)
    nodestep(1,1) = i1;
    nodestep(2,1) = i2;
    nodestep(3,1) = i3;
    nodestep(4,1) = i4;
    nodestep(5,1) = i5;
    nodestep(6,1) = i6;
    nodestep(7,1) = i7;
    nodestep(8,1) = i8;
    nodesteptop(1,1) = it1;
    nodesteptop(2,1) = it2;
    nodesteptop(3,1) = it3;
    nodesteptop(4,1) = it4;
    nodesteptop(5,1) = it5;
    
    % Assigning joint properties for top-chord joints
    if stmtype == 1
        for i = 1:load_num
            jointdata(i,2) = load(i,1);
            jointdata(i,3) = (h - offset_top)/12;
            if (l_res == 1) && (r_res == 1)
                if i == 1
                    jointdata(i,4) = -1*(FEMi/((h - offset_top - offset_bot)/12));
                end
                if i == load_num
                    jointdata(i,4) = (FEMf/((h - offset_top - offset_bot)/12));
                end
            else
                jointdata(i,4) = 0;
            end
            jointdata(i,5) = -1*(load(i,2));
        end
    end
        
    if stmtype == 2
         for i = 1:load_num
             jointdata((i*2),2) = load(i,1);
             jointdata((i*2),3) = (h - offset_top)/12;
             jointdata((i*2),5) = -1*(load(i,2));
         end
         for i = 1:load_num
             jointdata((i*2) + 1,2) = jointdata((i*2),2) + ((nodesteptop(i,1))/12);
             jointdata((i*2) + 1,3) = (h - offset_top)/12;
             
             if i == 1
                 jointdata(i,2) = (offset_hor/12) + ((nodesteptop(1,1))/12);
                 jointdata(i,3) = (h - offset_top)/12;
                 jointdata(i,4) = -1*(FEMi/((h - offset_top - offset_bot)/12));
             end
             if i == (2*load_num)
                 jointdata(i + 1,2) = (jointdata((i*2),2)) + ((nodesteptop(5,1))/12);
                 jointdata(i + 1,3) = (h - offset_top)/12;
                 jointdata(i + 1,4) = (FEMf/((h - offset_top - offset_bot)/12));   
             end
         end
    end

    % Assigning joint properties for bottom-chord joints
    if stmtype == 1
        temp_count_1 = 1;
        for i = (load_num + 2):(j - 1)
            jointdata(i,2) = load(temp_count_1,1) + ((nodestep(temp_count_1,1))/12);
            jointdata(i,3) = (offset_bot/12);
            temp_count_1 = temp_count_1 + 1;
        end
    end  
    if stmtype == 2
        temp_count_1 = 1;
        for i = (2*load_num + 3):(j - 1)
            jointdata(i,2) = jointdata(temp_count_1,2) + ((nodestep(temp_count_1,1))/12);
            jointdata(i,3) = (offset_bot/12);
            temp_count_1 = temp_count_1 + 1;
        end
    end

    % Assigning joint properties for left end boundary
    if stmtype == 1
        jointdata((load_num + 1),2) = (offset_hor/12);
        jointdata((load_num + 1),3) = (offset_bot/12);
        jointdata((load_num + 1),6) = 1;
        jointdata((load_num + 1),7) = 1;
        if (l_res == 1) && (r_res == 1)
            jointdata((load_num + 1),4) = (FEMi/((h - offset_top - offset_bot)/12));
        end
    end
    if stmtype == 2
        jointdata((2*load_num + 2),2) = (offset_hor/12);
        jointdata((2*load_num + 2),3) = (offset_bot/12);
        jointdata((2*load_num + 2),6) = 1;
        jointdata((2*load_num + 2),7) = 1;
        if (l_res == 1) && (r_res == 1)
            jointdata((2*load_num + 2),4) = (FEMi/((h - offset_top - offset_bot)/12));
        end
    end


    % Assigning joint properties for right end boundary
    jointdata(j,2) = l - (offset_hor/12);
    jointdata(j,3) = (offset_bot/12);
    jointdata(j,6) = 1;
    jointdata(j,7) = 1;
    if (l_res == 1) && (r_res == 1)
        jointdata(j,4) = -1*(FEMf/((h - offset_top - offset_bot)/12));
    end 
    
    % Assigning member properties for postively sloped diagonal members
    if stmtype == 1
        temp_count_2 = 1;
        for i = 1:(load_num)
            memberdata(i,2) = jointdata((load_num + temp_count_2),1);
            memberdata(i,3) = jointdata(temp_count_2,1);
            temp_count_2 = temp_count_2 + 1;
        end
    end
    if stmtype == 2
        temp_count_2 = 1;
        for i = 1:((2*load_num) + 1)
            memberdata(i,2) = jointdata(((2*load_num + 1) + temp_count_2),1);
            memberdata(i,3) = jointdata(temp_count_2,1);
            temp_count_2 = temp_count_2 + 1;
        end
    end

    % Assigning member properties for negatively sloped diagonal members
    if stmtype == 1
        temp_count_3 = 1;
        for i = (load_num + 1):(2*load_num)
            memberdata(i,2) = jointdata((load_num + temp_count_3 + 1),1);
            memberdata(i,3) = jointdata(temp_count_3,1);
            temp_count_3 = temp_count_3 + 1;
        end
    end
    if stmtype == 2
        temp_count_3 = 1;
        for i = ((2*load_num) + 2):(2*((2*load_num) + 1))
            memberdata(i,2) = jointdata(((2*load_num + 1) + temp_count_3 + 1),1);
            memberdata(i,3) = jointdata(temp_count_3,1);
            temp_count_3 = temp_count_3 + 1;
        end
    end
    
    % Assigning member properties for top-chord members
    if stmtype == 1
    temp_count_4 = 1;
        for i = ((2*load_num) + 1):((3*load_num) - 1)
            memberdata(i,2) = jointdata(temp_count_4,1);
            memberdata(i,3) = jointdata((temp_count_4 + 1),1);
            temp_count_4 = temp_count_4 + 1;
        end
    end
    if stmtype == 2
    temp_count_4 = 1;
        for i = ((2*((2*load_num) + 1)) + 1):(((2*((2*load_num) + 1)) + 1) + (2*load_num))
            memberdata(i,2) = jointdata(temp_count_4,1);
            memberdata(i,3) = jointdata((temp_count_4 + 1),1);
            temp_count_4 = temp_count_4 + 1;
        end
    end

    % Assigning member properties for bottom-chord members
    if stmtype == 1
        temp_count_5 = 1;
        for i = (3*load_num):m
            memberdata(i,2) = jointdata(temp_count_5 + load_num,1);
            memberdata(i,3) = jointdata((temp_count_5 + load_num + 1),1);
            temp_count_5 = temp_count_5 + 1;
        end
    end
    if stmtype == 2
        temp_count_5 = 1;
        for i = (((2*((2*load_num) + 1)) + 1) + (2*load_num) + 1):m
           memberdata(i,2) = jointdata(((2*load_num + 1) + temp_count_5),1);
           memberdata(i,3) = jointdata(((2*load_num + 1) + temp_count_5 + 1),1);
           temp_count_5 = temp_count_5 + 1;
        end 
    end
     
    % Calculating member lengths
    for i = 1:m
        memberdata(i,6) = (((((jointdata(memberdata(i,3),2)) - (jointdata(memberdata(i,2),2)))^2) + ...
                            (((jointdata(memberdata(i,3),3)) - (jointdata(memberdata(i,2),3)))^2))^0.5);
    end
    
%% Optimization

            % Counter definitions for DOF
            counter5 = 1;
            counter6 = 0;
            counter7 = 1;
            counter8 = 1;     % Counts rows in [dof]
            counter9 = 1;     % Counts DOF in [dof]
            counter10 = 1;    % Counts rows in [R]

            % Counting number of restraints in truss system
            for counter5 = 1:j
                counter6 = jointdata(counter5,6) + jointdata(counter5,7) + counter6;
                counter5 = counter5 + 1;
            end

            nr = counter6;          % Defining number of restraints
            NDOF = (2*j) - nr;      % Defining number of degrees of freedom
            counter12 = NDOF + 1;   % Counts number of restraints in the [R]
            dof = zeros(NDOF, 3);   % Defining empty [dof]
            R = zeros(nr, 3);       % Defining empty [R] 

            for counter7 = 1:j

            % Checking restain in the x direction

            % Populating [dof] accordingly
            if jointdata(counter7,6) == 0
                       dof(counter8,1) = counter7;
                       dof(counter8,2) = 1;
                       dof(counter8,3) = counter9;
                       counter8 = counter8 + 1;
                       counter9 = counter9 + 1;

            % Populating [R] accordingly
            else
                        R(counter10,1) = counter7;
                        R(counter10,2) = 1;
                        R(counter10,3) = counter12;
                        counter10 = counter10 + 1;
                        counter12 = counter12 + 1;
            end

            % Checking restain in the y direction

            % Populating [dof] accordingly
            if jointdata(counter7,7) == 0
                       dof(counter8,1) = counter7;
                       dof(counter8,2) = 2;
                       dof(counter8,3) = counter9;
                       counter8 = counter8 + 1;
                       counter9 = counter9 + 1;

            % Populating [R] accordingly  
            else
                        R(counter10,1) = counter7;
                        R(counter10,2) = 2;
                        R(counter10,3) = counter12;
                        counter10 = counter10 + 1;
                        counter12 = counter12 + 1;
            end
            end
            
            % Counter definitions for Stiffness and Transformation Matrices
            counter13 = 1;

            % For each member, calculates cos, sin, [T], [k] and [K], respectively.
            for counter13 = 1:m
                c = ((jointdata(memberdata(counter13,3),2)) - (jointdata(memberdata(counter13,2),2)))/(memberdata(counter13,6));
                s = ((jointdata(memberdata(counter13,3),3)) - (jointdata(memberdata(counter13,2),3)))/(memberdata(counter13,6));
                T{counter13} = [c, s, 0, 0; -s, c, 0, 0; 0, 0, c, s; 0, 0, -s, c];
                localk{counter13} = (((memberdata(counter13,4))*(memberdata(counter13,5)))/(memberdata(counter13,6))).*[1, 0, -1, 0; 0, 0, 0, 0;...
                            -1, 0, 1, 0; 0, 0, 0, 0];
                globalk{counter13} = (transpose(T{counter13}))*(localk{counter13})*(T{counter13});
            end

            % Counter definitions for Member ID Vectors
            counter14 = 1;
            counter15 = 1;
            counter16 = 1;

            for counter14 = 1:m

            % Determining beginning and ending joints, respectively.
                jb = memberdata(counter14,2);
                je = memberdata(counter14,3);

            % Populating [dof]
            for counter15 = 1:NDOF

            % Checking if DOF corresponds to beginning joint
            if dof(counter15,1) == jb
            if dof(counter15,2) == 1
                            ID1 = dof(counter15,3);
            else
                            ID2 = dof(counter15,3);
            end
            end

            % Checking if DOF corresponds to ending joint
            if dof(counter15,1) == je
            if dof(counter15,2) == 1
                            ID3 = dof(counter15,3);
            else
                            ID4 = dof(counter15,3);
            end
            end
            end

            % Populating [R]
            for counter16 = 1:nr

                % Checking if DOF corresponds to beginning joint
                if R(counter16,1) == jb
                    if R(counter16,2) == 1
                            ID1 = R(counter16,3);
                        else
                            ID2 = R(counter16,3);
                        end
                    end

                    % Checking if DOF corresponds to ending joint
                    if R(counter16,1) == je
                        if R(counter16,2) == 1
                            ID3 = R(counter16,3);
                        else
                            ID4 = R(counter16,3);
                    end
                end
            end
                V{counter14} = [ID1; ID2; ID3; ID4]; % Defining Member ID Vector
            end

            S = zeros(NDOF,NDOF);   % Defining empty [S]

            % Counter definitions for Structural Stiffness Matrix
            counter17 = 1;
            counter18 = 1;
            counter19 = 1;

            % Using Member ID Vectors to populate [S]
            for counter17 = 1:m
                for counter18 = 1:4
                    for counter19 = 1:4
                        if (V{counter17}(counter18) <= NDOF) && (V{counter17}(counter19) <= NDOF)
                            column = V{counter17}(counter18);
                            row = V{counter17}(counter19);
                            S(row,column) = (S(row,column)) + (globalk{counter17}(counter19,counter18));
                        else
                            S = S;
                        end
                    end
                end
            end

            % Counter definition for P Matrix
            counter20 = 1;

            P = zeros(NDOF,1);  % Defining empty [P]

            % Populating [P] based on [dof] and inputted joint data
            for counter20 = 1:NDOF
            if dof(counter20,2) == 1
                    P(counter20,1) = jointdata((dof(counter20,1)),4);
            else
                    P(counter20,1) = jointdata((dof(counter20,1)),5);
            end
            end

            d = zeros(NDOF,1); % Defining empty [d]
            d = S\P;           % Calculating [d]

            % Counter definitions for Displacement Matrices
            counter21 = 1;
            counter22 = 1;
            counter27 = 1;

            for counter21 = 1:m

            % For each member, defining empty [v] 
               globalu{counter21} = [0;0;0;0];

            % Populating [v] based Member ID Vector
            for counter22 = 1:4
            if V{counter21}(counter22) <= NDOF
                       globalu{counter21}(counter22) = d(V{counter21}(counter22));
            end
            end
            end

            ux = zeros(j,1);   % Defining empty [ux]
            uy = zeros(j,1);   % Defining empty [uy]

            % Populating [ux] and [uy] based on beginning and end joints and [globalu]
            for counter27 = 1:m
                jointbe = memberdata(counter27,2);
                jointee = memberdata(counter27,3);
                ux(jointbe,1) = globalu{counter27}(1);
                ux(jointee,1) = globalu{counter27}(3);
                uy(jointbe,1) = globalu{counter27}(2);
                uy(jointee,1) = globalu{counter27}(4);
            end

            % Defining counters for Force Matrices
            counter23 = 1;

            % For each member, calculating [u], [Q] and [F]
            for counter23 = 1:m
                u{counter23} = (T{counter23})*(globalu{counter23});
                Q{counter23} = (localk{counter23})*(u{counter23});
                F{counter23} = (transpose(T{counter23}))*(Q{counter23});
            end

            Rx = zeros(j,1);     % Defining empty [Rx]
            Ry = zeros(j,1);     % Defining empty [Ry]

            % Defining counters for Reactions
            counter25 = 1;
            counter26 = 1;

            % Populating [Rx] and [Ry] based on beginning and end joints, and [F]
            for counter25 = 1:m
                jointb = memberdata(counter25,2);
                jointe = memberdata(counter25,3);
                Rx(jointb,1) = Rx(jointb,1) + F{counter25}(1);
                Rx(jointe,1) = Rx(jointe,1) + F{counter25}(3);
                Ry(jointb,1) = Ry(jointb,1) + F{counter25}(2);
                Ry(jointe,1) = Ry(jointe,1) + F{counter25}(4);
            end

            % Calibrating [Rx] and [Ry] with respect to [P]
            for counter26 = 1:j
                Rx(counter26,1) = Rx(counter26,1) - jointdata(counter26,4);
                Ry(counter26,1) = Ry(counter26,1) - jointdata(counter26,5);
            end

            % Defining counter for Joint Output
            counter24 = 1;

            joint = zeros(j,1);  % Defining empty [joint]

            % Populating [joint]
            for counter24 = 1:j
                joint(counter24,1) = round(counter24);
            end
            % Defining conter for Member Output
            counter25 = 1;

            member = zeros(m,1);       % Defining empty [member]
            Axial_Force = zeros(m,1);  % Defining empty [Axial_Force]
            Stress = zeros(m,1);       % Defining empty [Stress]

            % Populating [member], [Axial_Force] and [Stress]
            for counter25 = 1:m
                    member(counter25,1) = round(counter25);
                    Axial_Force(counter25,1) = (((Q{counter25}(1))^2)^0.5);
                    Stress(counter25,1) = (Axial_Force(counter25,1))/(memberdata(counter25,5));
                    if (Q{counter25}(1)) > 0
                            Axial_Force(counter25,1) = -1*(Axial_Force(counter25,1));
                    end
            end
            
            Q_total = 0;
            L_total = 0;
            QL_total = 0;
            A_total = 0;

            for i = 1:m
                Q_total = Q_total + abs(Axial_Force(i,1));
                L_total = L_total + memberdata(i,6);
                QL_total = QL_total + abs((memberdata(i,6))*(Axial_Force(i,1)));
                A_total = A_total + ((abs(Axial_Force(i,1)))/fy)*memberdata(i,6)*12;
            end
            
            if stmtype == 1
                big_temp(maincounter,1) = Q_total;
                big_temp(maincounter,2) = L_total;
                big_temp(maincounter,3) = QL_total;
                big_temp(maincounter,4) = nodestep(1,1);
                big_temp(maincounter,5) = nodestep(2,1);
                big_temp(maincounter,6) = nodestep(3,1);
                big_temp(maincounter,7) = A_total;
            end
           
            if stmtype == 2
                big_temp(maincounter,1) = Q_total;
                big_temp(maincounter,2) = L_total;
                big_temp(maincounter,3) = QL_total;
                big_temp(maincounter,4) = nodestep(1,1);
                big_temp(maincounter,5) = nodestep(2,1);
                big_temp(maincounter,6) = nodestep(3,1);
                big_temp(maincounter,7) = nodestep(4,1);
                big_temp(maincounter,8) = nodestep(5,1);
                big_temp(maincounter,9) = nodestep(6,1);
                big_temp(maincounter,10) = nodestep(7,1);
                big_temp(maincounter,11) = nodestep(8,1);
                big_temp(maincounter,12) = nodesteptop(1,1);
                big_temp(maincounter,13) = nodesteptop(2,1);
                big_temp(maincounter,14) = nodesteptop(3,1);
                big_temp(maincounter,15) = nodesteptop(4,1);
                big_temp(maincounter,16) = nodesteptop(5,1);
                big_temp(maincounter,17) = A_total;
            end
            maincounter = maincounter + 1;
    end
    end
    end
    end
    end
    end
    end
    end
    end
    end
    end
    end
    end

    if stmtype == 1
        mainmain = zeros(maincounter,6);
        for i = 1:maincounter
            mainmain(i,1) = big_temp(i,1);
            mainmain(i,2) = big_temp(i,2);
            mainmain(i,3) = big_temp(i,3);
            mainmain(i,4) = big_temp(i,4);
            mainmain(i,5) = big_temp(i,5);
            mainmain(i,6) = big_temp(i,6);
            mainmain(i,7) = big_temp(i,7);
        end
    end
    
    if stmtype == 2
        mainmain = zeros(maincounter,16);
        for i = 1:maincounter
            mainmain(i,1) = big_temp(i,1);
            mainmain(i,2) = big_temp(i,2);
            mainmain(i,3) = big_temp(i,3);
            mainmain(i,4) = big_temp(i,4);
            mainmain(i,5) = big_temp(i,5);
            mainmain(i,6) = big_temp(i,6);
            mainmain(i,7) = big_temp(i,7);
            mainmain(i,8) = big_temp(i,8);
            mainmain(i,9) = big_temp(i,9);
            mainmain(i,10) = big_temp(i,10);
            mainmain(i,11) = big_temp(i,11);
            mainmain(i,12) = big_temp(i,12);
            mainmain(i,13) = big_temp(i,13);
            mainmain(i,14) = big_temp(i,14);
            mainmain(i,15) = big_temp(i,15);
            mainmain(i,16) = big_temp(i,16);
            mainmain(i,17) = big_temp(i,17);
        end
    end
    
    Mmax = max(mainmain);
    Mmin = min(mainmain);
    Qmax = Mmax(1);
    Lmax = Mmax(2);
    QLmin = Mmin(3);
    
    if stmtype == 1
        Amin = Mmin(7);
    end
    if stmtype == 2
        Amin = Mmin(17);
    end
%% Diagram
    subplot(1,2,1)
    hold on
    axis([-36 (((l*12) + 36)) -12 (h + 36)])
    axis('square')
    
    % Plotting outline of beam
    plot([0 (l*12)], [0 0], 'Linewidth', 4', 'Color','k')
    plot([0 0], [0 h], 'Linewidth', 4, 'Color','k')
    plot([(l*12) (l*12)], [0 h], 'Linewidth', 4, 'Color','k')
    plot([0 (l*12)], [h h], 'Linewidth', 4, 'Color','k')
    title('Initial STM Layout')
    
    % Plotting point loads and bearing plates
    for i = 1:load_num
        x1 = load(i,1)*12;
        x2 = x1;
        y1 = h + BearingPL_depth;
        y2 = (h + BearingPL_depth) + ((load(i,2))*Pointload_scale);
        x1a1 = x1;
        x2a1 = x1 - Arrow_size;
        y1a1 = h + BearingPL_depth;
        y2a1 = h + BearingPL_depth + Arrow_size; 
        x1a2 = x1;
        x2a2 = x1 + Arrow_size;
        y1a2 = h + BearingPL_depth;
        y2a2 = h + BearingPL_depth + Arrow_size;
        
        plot([x1 x2], [y1 y2], 'Linewidth', 1', 'Color','k')
        plot([x1a1 x2a1], [y1a1 y2a1], 'Linewidth', 1', 'Color','k')
        plot([x1a2 x2a2], [y1a2 y2a2], 'Linewidth', 1', 'Color','k')
        
        x1b = x1 - (load(i,3))*0.5;
        y1b = h;
        x2b = x1 - (load(i,3))*0.5;
        y2b = h + BearingPL_depth;
        x3b = x1 + (load(i,3))*0.5;
        y3b = h;
        x4b = x1 + (load(i,3))*0.5;
        y4b = h + BearingPL_depth;
        
        plot([x1b x2b], [y1b y2b], 'Linewidth', 1', 'Color','m')
        plot([x3b x4b], [y3b y4b], 'Linewidth', 1', 'Color','m')
        plot([x1b x3b], [y1b y3b], 'Linewidth', 1', 'Color','m')
        plot([x2b x4b], [y2b y4b], 'Linewidth', 1', 'Color','m')
    end
    
    % Plotting distributed loads
    for i = 1:w_num
        x1 = temp_w(i,1)*12;
        y1 = h;
        x2 = temp_w(i,1)*12;
        y2 = h + (temp_w(i,4)*Distload_scale);
        x3 = temp_w(i,2)*12;
        y3 = h;
        x4 = temp_w(i,2)*12;
        y4 = h + (temp_w(i,4)*Distload_scale);
        
        plot([x1 x2], [y1 y2], 'Linewidth', 1', 'Color','c')
        plot([x3 x4], [y3 y4], 'Linewidth', 1', 'Color','c')
        plot([x1 x3], [y1 y3], 'Linewidth', 1', 'Color','c')
        plot([x2 x4], [y2 y4], 'Linewidth', 1', 'Color','c')
    end
    
    % Plotting FEMi arrow
    if l_res == 1
        a = 12;
        b = (h/2);
        x0 = -12;
        y0 = (h/2);
        t = pi/2:0.01:3*pi/2;
        x = x0 + a*cos(t);
        y = y0 + b*sin(t);
        x1a3 = -12;
        x2a3 = x1a3 - Arrow_size;
        y1a3 = 0;
        y2a3 = Arrow_size; 
        x1a4 = x1a3;
        x2a4 = x1a3 - Arrow_size;
        y1a4 = 0;
        y2a4 = -Arrow_size;
        plot([x1a3 x2a3], [y1a3 y2a3], 'Linewidth', 1', 'Color','c')
        plot([x1a4 x2a4], [y1a4 y2a4], 'Linewidth', 1', 'Color','c')
        plot(x,y, 'Linewidth', 1', 'Color','c');   
    end
    
    % Plotting FEMf arrow
    if r_res == 1
        a = 12;
        b = (h/2);
        x0 = (l*12) + 12;
        y0 = (h/2);
        t = -pi/2:0.01:pi/2;
        x = x0 + a*cos(t);
        y = y0 + b*sin(t);
        x1a3 = (l*12) + 12;
        x2a3 = x1a3 + Arrow_size;
        y1a3 = 0;
        y2a3 = Arrow_size; 
        x1a4 = x1a3;
        x2a4 = x1a3 + Arrow_size;
        y1a4 = 0;
        y2a4 = -Arrow_size;
        plot([x1a3 x2a3], [y1a3 y2a3], 'Linewidth', 1', 'Color','c')
        plot([x1a4 x2a4], [y1a4 y2a4], 'Linewidth', 1', 'Color','c')
        plot(x,y, 'Linewidth', 1', 'Color','c');
        plot(x,y, 'Linewidth', 1', 'Color','c');
    end
    
    % Plotting struts and ties of STM
    for i = 1:m
        x1 = 12*jointdata((memberdata(i,2)),2);
        x2 = 12*jointdata((memberdata(i,3)),2);
        y1 = 12*jointdata((memberdata(i,2)),3);
        y2 = 12*jointdata((memberdata(i,3)),3);
        if Axial_Force(i,1) > 0
            plot([x1 x2], [y1 y2], 'Linewidth', 1', 'Color','b')
        end
        if Axial_Force(i,1) < 0
            plot([x1 x2], [y1 y2], 'Linewidth', 1', 'Color','r')
        end
    end
    
    hold on
    subplot(1,2,2)
        Q = mainmain(:,1);
        L = mainmain(:,2);
        scatter(L,Q)
        title('F vs. L of Node Increments')
        xlabel('Total Truss Length (ft)')
        ylabel('Total Truss Force (kips')