'Apply on a scalar field.
The action of a vector field on a scalar field is a directional
differentiation.
If the argument is not a scalar field the behaviour is undefined.'

def __call__(self, field):

    if (covariant_order(field) != 0):
        raise NotImplementedError()
    field = vectors_in_basis(field, self._coord_sys)
    wrt_vector = self._coord_sys.base_vector(self._index)
    wrt_scalar = self._coord_sys.coord_function(self._index)
    vectors = list(field.atoms(BaseVectorField))
    d_funcs = [Function(('_#_%s' % i))(wrt_scalar) for (i, b) in enumerate(vectors)]
    d_result = field.subs(list(zip(vectors, d_funcs)))
    d_result = wrt_vector(d_result)
    d_result = d_result.subs(list(zip(d_funcs, vectors)))
    derivs = []
    for v in vectors:
        d = Add(*[(self._christoffel[(k, wrt_vector._index, v._index)] * v._coord_sys.base_vector(k)) for k in range(v._coord_sys.dim)])
        derivs.append(d)
    to_subs = [wrt_vector(d) for d in d_funcs]
    result = d_result.subs(list(zip(to_subs, derivs)))
    return result
