function [elem, QGlobal_elem, KMaterialGlobal_elem, KGeomGlobal_elem, QLocal_elem]=elemStateDetermin(elemIdx,elem,numEquations, v)
dofs_elem = numEquations{elemIdx};

T = computeRotationMatrix(elem);
u = T * v(dofs_elem);

if strcmp(elem{1}, '2dElasticBeam')
    [qBar, KBar] = elemStateDetermin_2dElasticBeam(elem, u);

    l = norm(elem{4} - elem{3});
    [LTransfo, Kgeom] = computeTransfoBasicToLocal(elem{8}, l, u, qBar);

    QLocal_elem = LTransfo' * qBar;
    QGlobal_elem = T' * QLocal_elem;
    KMaterialGlobal_elem = T' * (LTransfo' * KBar * LTransfo) * T;
    KGeomGlobal_elem = T' * Kgeom * T;
    
elseif strcmp(elem{1}, 'spring')
    [elem, QLocal_elem, KLocal]=elemStateDetermin_spring(elem,u);
    QGlobal_elem = QLocal_elem;
    KMaterialGlobal_elem = KLocal;
    KGeomGlobal_elem = zeros(2,2);
    
elseif strcmp(elem{1}, '2dModElasticBeam')
    [qBar, KBar] = elemStateDetermin_2dModElasticBeam(elem, u);

    l = norm(elem{4} - elem{3});
    [LTransfo, Kgeom] = computeTransfoBasicToLocal(elem{12}, l, u, qBar);

    QLocal_elem = LTransfo' * qBar;
    QGlobal_elem = T' * QLocal_elem;
    KMaterialGlobal_elem = T' * (LTransfo' * KBar * LTransfo) * T;
    KGeomGlobal_elem = T' * Kgeom * T;  
    
elseif strcmp(elem{1}, '2dDispBeamColumn')
    [qBar, KBar] = elemStateDetermin_2dDispBeamColumn(elem, u);

    l = norm(elem{4} - elem{3});
    [LTransfo, Kgeom] = computeTransfoBasicToLocal(elem{10}, l, u, qBar);

    QLocal_elem = LTransfo' * qBar;
    QGlobal_elem = T' * QLocal_elem;
    KMaterialGlobal_elem = T' * (LTransfo' * KBar * LTransfo) * T;
    KGeomGlobal_elem = T' * Kgeom * T; 
 
end
end


function [qBar, KBar]=elemStateDetermin_2dElasticBeam(elem,u)
l = norm(elem{4} - elem{3});
E = elem{5};
A = elem{6};
I = elem{7};

KBar = [
    E * A / l,           0,                0;
    0,           4 * E * I / l,     2 * E * I / l;
    0,           2 * E * I / l,     4 * E * I / l
];

uBar = computeUbar(elem{8}, l, u);

qBar = KBar * uBar;
end


function [qBar, KBar] = elemStateDetermin_2dModElasticBeam(elem, u)
    E = elem{5};           
    A = elem{6};           
    Ie = elem{11};         
    l = norm(elem{4} - elem{3});  
    S22mod = elem{7};      
    S23mod = elem{8};      
    S32mod = elem{9};      
    S33mod = elem{10};      

    KBar = [
        E * A / l, 0, 0;
        0, S22mod * E * Ie / l, S23mod * E * Ie / l;
        0, S32mod * E * Ie / l, S33mod * E * Ie / l
    ];

    uBar = computeUbar(elem{12}, l, u);

    qBar = KBar * uBar;

end


function [elem, QLocal, KLocal]=elemStateDetermin_spring(elem,u)
% Compute spring rotation
theta = u(2) - u(1);

% Unpack the spring material properties from elem{3}
spring_MatProps_input = elem{4};

% Check for correct number of material properties
if numel(spring_MatProps_input) ~= 22
    error('Invalid number of spring material properties.');
end

% Extract input material properties
theta_previous_input = spring_MatProps_input(1);
M_previous_input = spring_MatProps_input(2);
theta_M0_currentPos_input = spring_MatProps_input(3);
theta_M0_currentNeg_input = spring_MatProps_input(4);
theta_M0_projected_input = spring_MatProps_input(5);
ke_s_input = spring_MatProps_input(6);
ks_s_input = spring_MatProps_input(7);
kpc_s_input = spring_MatProps_input(8);
thetaY_s_input = spring_MatProps_input(9);
thetaC_s_input = spring_MatProps_input(10);
thetaU_s_input = spring_MatProps_input(11);
My_input = spring_MatProps_input(12);
Mu_input = spring_MatProps_input(13);
MmaxPos_input = spring_MatProps_input(14);
MmaxNeg_input = spring_MatProps_input(15);
yieldFlag_Pos_input = spring_MatProps_input(16);
cappingFlag_Pos_input = spring_MatProps_input(17);
yieldFlag_Neg_input = spring_MatProps_input(18);
cappingFlag_Neg_input = spring_MatProps_input(19);
reversalFlag_input = spring_MatProps_input(20);
residualFlag_input = spring_MatProps_input(21);
Di_previous_input = spring_MatProps_input(22);

% Call the constitutive law for the spring (you'll need to implement this function)
[M, k, MmaxPos, MmaxNeg, yieldFlag_Pos, cappingFlag_Pos, ...
    yieldFlag_Neg, cappingFlag_Neg, reversalFlag, residualFlag, ...
    theta_M0_currentPos, theta_M0_currentNeg, theta_M0_projected, Di] = trilinearResponse( ...
    theta, theta_previous_input, M_previous_input, theta_M0_currentPos_input, ...
    theta_M0_currentNeg_input, theta_M0_projected_input, ke_s_input, ks_s_input, ...
    kpc_s_input, thetaY_s_input, thetaC_s_input, thetaU_s_input, My_input, ...
    Mu_input, MmaxPos_input, MmaxNeg_input, yieldFlag_Pos_input, cappingFlag_Pos_input, ...
    yieldFlag_Neg_input, cappingFlag_Neg_input, reversalFlag_input, residualFlag_input, Di_previous_input);

% Spring resisting force and tangent stiffness
QLocal = [-M; M];
KLocal = [k, -k; -k, k];

% Assign the updated material properties to elem{4}
spring_MatProps_Output = [ ...
    theta, M, theta_M0_currentPos, theta_M0_currentNeg, theta_M0_projected, ...
    ke_s_input, ks_s_input, kpc_s_input, thetaY_s_input, thetaC_s_input, ...
    thetaU_s_input, My_input, Mu_input, MmaxPos, MmaxNeg, yieldFlag_Pos, ...
    cappingFlag_Pos, yieldFlag_Neg, cappingFlag_Neg, reversalFlag, residualFlag, Di];

elem{5} = spring_MatProps_Output;

end


function [qBar, KBar] = elemStateDetermin_2dDispBeamColumn(elem, u)
    l = norm(elem{4} - elem{3});

    uBar = computeUbar(elem{10}, l, u);
    
    xi = elem{5};   % x/l
    wt = elem{6};   % weight
    wtL = wt * l;   
    numSection = length(xi);   

    % Initialize element state variables
    qBar = zeros(3, 1);
    KBar = zeros(3, 3);

    % Loop over the integration sections
    for iSec = 1:numSection
        Bbar = [1/l, 0, 0;
                0, (6*xi(iSec) - 4)/l, (6*xi(iSec) - 2)/l];
        
        ds_iSec = Bbar * uBar;

        % Section state determination 
        materialInput = elem{7};
        coordinates_fibers_columns = elem{8};
        area_fibers_columns = elem{9};
        [Qsr, ks] = sectionStateDetermin_uniaxialFiber(materialInput, coordinates_fibers_columns, area_fibers_columns, ds_iSec);

        % Integrate over the element length
        qBar = qBar + Bbar' * Qsr * wtL(iSec);
        KBar = KBar + Bbar' * ks * Bbar * wtL(iSec);
    end
end