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, M_previous, theta_M0_currentPos, theta_M0_currentNeg, ...
    theta_M0_projected, ke, ks, kpc, thetaY, thetaC, thetaU, My, Mu, ...
    MmaxPos, MmaxNeg, yieldFlag_Pos, cappingFlag_Pos, yieldFlag_Neg, ...
    cappingFlag_Neg, reversalFlag, residualFlag, Di_previous)

    % Determine the direction of the current increment
    if theta >= theta_previous
        Di = 1;
    else
        Di = -1;
    end
    
    % Trial elastic moment
    M = M_previous + ke * (theta - theta_previous);
    
    % Quick return if residual state is active
    if residualFlag == 1
        M = 0;
        k = 10;
        return;
    end
    
    % Inelastic behavior for positive rotation
    if Di == 1 && M > MmaxPos
        yieldFlag_Neg = 0;
        cappingFlag_Neg = 0;
        
        if ~(yieldFlag_Pos || cappingFlag_Pos)
            theta_M0_currentPos = theta_M0_projected;
        end
        
        if abs(theta - theta_M0_currentPos) >= thetaY && abs(theta - theta_M0_currentPos) < thetaC
            % Hardening path
            yieldFlag_Pos = 1;
            theta_M0_currentNeg = theta_M0_projected;
            M = Di * My + ks * (theta - theta_M0_currentPos - thetaY);
            MmaxNeg = -My;
        elseif abs(theta - theta_M0_currentPos) >= thetaC && abs(theta - theta_M0_currentPos) < thetaU
            % Softening path
            yieldFlag_Pos = 0;
            cappingFlag_Pos = 1;
            M = Di * Mu + kpc * (theta - theta_M0_currentPos - thetaC);
        else
            % Residual path
            residualFlag = 1;
            M = 0;
        end
        
        MmaxPos = M;
    end
    
    % Inelastic behavior for negative rotation
    if Di == -1 && M < MmaxNeg
        yieldFlag_Pos = 0;
        cappingFlag_Pos = 0;
        
        if ~(yieldFlag_Neg || cappingFlag_Neg)
            theta_M0_currentNeg = theta_M0_projected;
        end
        
        if abs(theta - theta_M0_currentNeg) >= thetaY && abs(theta - theta_M0_currentNeg) < thetaC
            % Hardening path
            yieldFlag_Neg = 1;
            theta_M0_currentPos = theta_M0_projected;
            M = Di * My + ks * (theta - theta_M0_currentNeg + thetaY);
            MmaxPos = My;
        elseif abs(theta - theta_M0_currentNeg) >= thetaC && abs(theta - theta_M0_currentNeg) < thetaU
            % Softening path
            yieldFlag_Neg = 0;
            cappingFlag_Neg = 1;
            M = Di * Mu + kpc * (theta - theta_M0_currentNeg + thetaC);
        else
            % Residual path
            residualFlag = 1;
            M = 0;
        end
        
        MmaxNeg = M;
    end
    
    % Check for unloading
    if Di_previous * Di < 0
        theta_M0_projected = theta_previous - M_previous / ke;
    end
    
    % Compute the tangent stiffness
    if theta == theta_previous
        k = 0;
    elseif M == M_previous || residualFlag == 1
        k = 1e4;
    else
        k = (M - M_previous) / (theta - theta_previous);
    end
end