# Implements plasticity algorithm

import scipy as sc

def radialReturn(deps, E, sigY, H = 0):
    """
    Simulate a traction experience.
    :param deps: A iteration vector for epsilon.
    :param E: Hook factor.
    :param sigY: Max elastic stress.
    :param H: Ecrouissage.
    :return: Three vectors (eps, sig, epsP)
    """
    # Announcement
    print("Computing results ...")

    # Prepare results
    n = len(deps) + 1
    eps = sc.zeros(n)
    sig = sc.zeros(n)
    epsP = sc.zeros(n)

    # Algorithm loop
    for k in range(0, n - 1):
        # Add deformation
        eps[k + 1] = eps[k] + deps[k]

        # Plastic prediction
        sig_pred = E*(eps[k + 1] - epsP[k])
        f_pred = abs(sig_pred - H*epsP[k]) - sigY

        # Plasticity test
        if f_pred < 0.0 :
            # Elastic domain
            epsP[k + 1] = epsP[k]
            sig[k + 1] = sig_pred
        else:
            # Radial return
            lda = f_pred / (E + H)
            epsP[k + 1] = epsP[k] + sc.sign(sig_pred - H*epsP[k])*lda
            sig[k + 1] = E*(eps[k + 1] - epsP[k + 1])

    # End of computation
    return (eps, sig, epsP)