"""Original reference implementation of the documented original-form softening law.
Homogeneous, incompressible, uniaxial neo-Hookean base; not a fitted material model.
Beta in the modified Abaqus law is fixed to zero, not silently imported into 1999.
"""
import math

def energy_density(stretch,shear_Pa):
    if not math.isfinite(stretch) or stretch<=0 or not math.isfinite(shear_Pa) or shear_Pa<=0:raise ValueError('Positive finite stretch and shear modulus required')
    return .5*shear_Pa*(stretch*stretch+2/stretch-3)

def softening(stretch,maximum_energy_Pa,shear_Pa=100000.,r=2.,scale_Pa=30000.):
    if not all(math.isfinite(x) for x in [maximum_energy_Pa,r,scale_Pa]) or maximum_energy_Pa<0 or r<=1 or scale_Pa<=0:raise ValueError('Require Wmax>=0, r>1 and positive energy-density scale')
    w=max(0.,energy_density(stretch,shear_Pa));wm=max(maximum_energy_Pa,w)
    eta=1-math.erf((wm-w)/scale_Pa)/r
    nominal_stress_Pa=eta*shear_Pa*(stretch-stretch**-2)
    return nominal_stress_Pa,wm,eta

def path(stretches,shear_Pa=100000.,r=2.,scale_Pa=30000.,initial_maximum_energy_Pa=0.):
    wm=initial_maximum_energy_Pa;out=[]
    for stretch in stretches:
        stress,wm,eta=softening(float(stretch),wm,shear_Pa,r,scale_Pa)
        out.append(dict(stretch=float(stretch),nominal_stress_Pa=stress,maximum_energy_density_Pa=wm,eta=eta))
    return out