%% Determines the values of Nu for all columns in the frame with considering the tangent modulus reduction
% E0 = scalar, ambient modulus of elasticity
% Ic = [n_stoerys × n_cols] moment of inertia of columns
% Hk = [n_storeys × 1] storey heights
% RL = [n_storeys × n_bays] beam left end connection rotational stiffness
% RR = [n_storeys × n_bays] beam right end connection rotational stiffness
% Rl = [1 × n_cols] ground connection rotational stiffness
% vc_buckle = [n_storeys × n_cols × n_storeys] ratio of upper to lower connection rotations for each column (1st and 2nd index) - DURING ROTATIONAL BUCKLING FOR EACH STOREY (3rd index)
% su = [n_storeys × n_cols] upper end splice fixity of columns (1 for continuous splice)
% sl = [n_storeys × n_cols] lower end splice fixity of columns (1 for continuous splice)
% incl_cols = toggle 1 to include column contribution to end fixity, 0 to exclude
% bool_inel = toggle 1 to include inelasticity tangent modulus considerations
% Ac = [n_stoerys × n_cols] area of columns
% Fy = scalar, ambient yield stress
% bool_prop = toggle 1 to consider pseudo-proportional loads (also specify Pl, Pu, X, mask), 0 for individual loading
% Pl = [n_stoerys × n_cols] lower bound loads
% Pu = [n_storeys × n_cols] upper bound loads
% X = [n_vars] vector of pseudo-proportional loading factors
% mask = [n_storeys × n_cols] boolean matrix of which loads to apply proportional loading on
function [Nu] = Iterate_Nu_Multi(E0,Ic,Hk,RL,RR,Rl,vc_buckle,su,sl,incl_cols,bool_inel,Ac,Fy,...
                    bool_prop,Pl,Pu,X,mask,bool_decomp,wul,wlu,...
                     bool_shear,shear_assumption,kappa,v0)
    
    %% Problem setup
    n_storeys = size(Ic,1);
    n_cols = size(Ic,2);
    tol = 0.001;
    Nu = zeros(n_storeys,n_cols);
    if isempty(wul) || isempty(wlu)
        wul = zeros(n_storeys,n_cols); % No drift during rotational buckling
        wlu = zeros(n_storeys,n_cols); % No drift during rotational buckling
    end
    max_iters = 250;
    Py = Fy*Ac;
    
    if bool_prop == 1
        % mask = [] will let it search using all variables
        if isempty(mask)
            mask = ones(size(Pl,1),size(Pl,2));
        end
        indices = find(mask==1);
        X_matrix = zeros(n_storeys,n_cols);
        for index = 1:length(indices)
            X_matrix(indices(index)) = X(index);
        end
    end

    for i = 1:n_storeys
        vc = vc_buckle(:,:,i); % Buckling mode is different for each storey
        
        for j = 1:n_cols
            
            weight = 1; % Convergence setting - 1 for Nu to move entirely to the new value; 0.5 to take a weighted average
            retry = 1;
            while retry > 0

                % Iteration 0 - Start with no load to get initial end fixity factors
                if bool_prop == 0
                    P = zeros(n_storeys,n_cols);
                else
                    P = Pl;
                end
                converge = 0;
                num_iters = 0;
                force_elastic = 0;

                % Iterate the value of Nu
                while (num_iters < max_iters && converge == 0)

                    % Update the end fixities
                    [ru, rl, phi, omega] = Get_End_Fixities(E0,Ic,Hk,P,RL,RR,Rl,vc,wul,wlu,su,sl,incl_cols,bool_inel,Ac,Fy,bool_decomp,...
                         bool_shear,shear_assumption,kappa,v0,1);
                    
                    % Update EI, kAG, neta
                    N = 0;
                    for ss = i:1:n_storeys
                        N = N + P(ss,j);
                    end
                    if bool_inel == 1
                        E = Get_E_IE(N,Ac(i,j),E0,Fy);
                        bool_buckle = 1;
                        if bool_buckle == 1
                            e = 0.5*(E0/E-1);
                        else
                            e = 0;
                        end
                        G = E0/(2*(1+v0) + 3*e);
                    else
                        E = E0;
                        G = E0/(2*(1+v0));
                    end
                    EI = E*Ic(i,j);
                    kAG = G*kappa*Ac(i,j);
                    if bool_shear == 1
                        eta = EI/Hk(i)^2/kAG;
                    else %%bool_shear == 0
                        eta = 0;
                    end
                    
                    % Estimate the new value of Nu based on the updated end fixities
                    K_approx = Get_K(ru(i,j), rl(i,j));
                    phiu = Get_Phiu(ru(i,j), rl(i,j), K_approx, eta, shear_assumption);
                    min_phiu = abs(phiu);
                    if phiu == -100 || phiu > 8 % If phiu doesn't initially converge or converges to a higher mode try multistart / brute force
                        phi_try = 0.01:0.05:11;
                        for kk = 1:length(phi_try)
                            phiu = Get_Phiu(ru(i,j), rl(i,j), 3.141592654/phi_try(kk), eta,shear_assumption);
                            if phiu > 0.01 && phiu<min_phiu
                                min_phiu = phiu;
                            end
                        end
                    end
                    if min_phiu == 100
                        phiu = -100;
                    else
                        phiu = min_phiu;
                    end
                    % Note: The value of phiu obtained from Get_Phiu does not include omega!! (it is not phi'u)

                    if phiu == -100
                        Nu(i,j) = -142; % phiu didn't converge
                        num_iters = max_iters;
                        continue;
                    end

                    if bool_inel == 1 && force_elastic == 0 % Inelastic analysis
                        Fe0 = phiu^2*(E0*Ic(i,j)/(Hk(i)^2*Ac(i,j))); % Based on the derivation, Fe0 shall be calculated using phiu, not phi'u
                        new_Nu = 0.85*Py(i,j)*((10 ^ (-1 / 7.38)) ^ (Fy / Fe0));
                    else % Elastic analysis
                        new_Nu = phiu^2*E0*Ic(i,j)/Hk(i)^2;
                    end

                    if abs(new_Nu - Nu(i,j)) < tol 
                        converge = 1;
                        if bool_inel == 1
                            if new_Nu < Py(i,j)/3 && bool_inel == 1 % Converged in the elastic range
                                force_elastic = 1; % For an elastic analysis; restart
                                num_iters = 0;
                                converge = 0;
                                new_Nu = 0;
                            end
                        end
                    end
                    Nu(i,j) = Nu(i,j) + (new_Nu-Nu(i,j))*weight; % Move Nu closer to the new value

                    % Update the loading matrix such that Nu occurs
                    if bool_prop == 0
                        P = zeros(n_storeys,n_cols);
                        P(i,j) = Nu(i,j);
                    else
                        P = Pl;
                        sumPl_above = 0;
                        sumPuPlX_above = 0;
                        for ss = i:n_storeys
                            sumPl_above = sumPl_above + Pl(ss,j);
                            sumPuPlX_above = sumPuPlX_above + (Pu(ss,j)-Pl(ss,j))*X_matrix(ss,j);
                        end
                        lambda_curr = (Nu(i,j) - sumPl_above)/sumPuPlX_above;
                        for index=1:length(indices)
                            if X_matrix(i,j) == 0
                                P(indices(index)) = Pl(indices(index));
                            else
                                P(indices(index)) = Pl(indices(index)) + (Pu(indices(index))-Pl(indices(index)))*X(index)*lambda_curr;
                            end
                        end
                    end
                    num_iters = num_iters + 1;
                end

                % Error 140 - Didn't converge?
                if converge == 0
                    Nu(i,j) = -140; 
                    retry = retry + 1;
                    weight = weight * 0.9;
                    if retry >= 15
                        retry = 0; % Give up
                    end
                else % Converged
                    retry = 0;
                end

            end
            
        end
    end
end