def hermite(m0,m1,t0,t1,u,domain=None):
    """Evaluate a cubic hermite curve at interpolation parameter u in [0,1].
    
    Endpoints are m0 and m1, with derivatives t0 and t1.  These are assumed to be Numpy arrays.
    
    Alternatively, if `domain` != None, then this will use domain=(a,b)
    as the interpolation domain
    """
    if domain is not None:
        assert isinstance(domain,(list,tuple)) and len(domain) == 2,"Need to provide a pair as a domain"
        scale = (domain[1]-domain[0])
        t = (u - domain[0])/scale
        return hermite(m0,m1,t0*scale,t1*scale,t)
    u2 = u**2
    u3 = u**3
    cm0 = 2*u3 - 3*u2 + 1
    cm1 = -2*u3 + 3*u2
    ct0 = u3 - 2*u2 + u
    ct1 = u3 - u2
    return cm0*m0 + cm1*m1 + ct0*t0 + ct1*t1