import numpy as np

def 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
):
    """
    Computes the moment and stiffness based on the trilinear material law.
    """
    # Determine the direction of the current increment
    Di = 1 if theta >= theta_previous else -1
    
    # 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 (
            M, k, MmaxPos, MmaxNeg, yieldFlag_Pos, cappingFlag_Pos, 
            yieldFlag_Neg, cappingFlag_Neg, reversalFlag, residualFlag, 
            theta_M0_currentPos, theta_M0_currentNeg, theta_M0_projected, Di
        )

    # Inelastic behavior for positive rotation
    if Di == 1 and M > MmaxPos:
        yieldFlag_Neg = 0
        cappingFlag_Neg = 0
        
        if not (yieldFlag_Pos or cappingFlag_Pos):
            theta_M0_currentPos = theta_M0_projected
        
        if 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
        elif 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
        
        MmaxPos = M

    # Inelastic behavior for negative rotation
    if Di == -1 and M < MmaxNeg:
        yieldFlag_Pos = 0
        cappingFlag_Pos = 0
        
        if not (yieldFlag_Neg or cappingFlag_Neg):
            theta_M0_currentNeg = theta_M0_projected
        
        if 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
        elif 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
        
        MmaxNeg = M

    # Check for unloading
    if Di_previous * Di < 0:
        theta_M0_projected = theta_previous - M_previous / ke

    # Compute the tangent stiffness
    if theta == theta_previous:
        k = 0
    elif M == M_previous or residualFlag == 1:
        k = 10
    else:
        k = (M - M_previous) / (theta - theta_previous)

    return (
        M, k, MmaxPos, MmaxNeg, yieldFlag_Pos, cappingFlag_Pos, 
        yieldFlag_Neg, cappingFlag_Neg, reversalFlag, residualFlag, 
        theta_M0_currentPos, theta_M0_currentNeg, theta_M0_projected, Di
    )


def bilinMonotonic(E,fy,ks,strain):
    # Trial stress
    sigma=E*strain
    tangentModulus=E

    # Plastic correction
    if abs(sigma) > fy:
        sigma = np.sign(sigma)*(fy+(abs(strain)-fy/E)*ks*E)
        tangentModulus = ks*E

    return sigma, tangentModulus