Convert a possiblty time-dependent operator into the nested-list
    structure required by QuTiP

def _time_dependent_to_qutip(
        op, full_space=None, time_symbol=symbols("t", real=True),
        convert_as='pyfunc'):
    """Convert a possiblty time-dependent operator into the nested-list
    structure required by QuTiP"""
    if full_space is None:
        full_space = op.space
    if time_symbol in op.free_symbols:
        op = op.expand()
        if isinstance(op, OperatorPlus):
            result = []
            for o in op.operands:
                if time_symbol not in o.free_symbols:
                    if len(result) == 0:
                        result.append(convert_to_qutip(o,
                                                       full_space=full_space))
                    else:
                        result[0] += convert_to_qutip(o, full_space=full_space)
            for o in op.operands:
                if time_symbol in o.free_symbols:
                    result.append(_time_dependent_to_qutip(o, full_space,
                                  time_symbol, convert_as))
            return result
        elif (
                isinstance(op, ScalarTimesOperator) and
                isinstance(op.coeff, ScalarValue)):
            if convert_as == 'pyfunc':
                func_no_args = lambdify(time_symbol, op.coeff.val)
                if {time_symbol, } == op.coeff.free_symbols:
                    def func(t, args):
                        # args are ignored for increased efficiency, since we
                        # know there are no free symbols except t
                        return func_no_args(t)
                else:
                    def func(t, args):
                        return func_no_args(t).subs(args)
                coeff = func
            elif convert_as == 'str':
                # a bit of a hack to replace imaginary unit
                # TODO: we can probably use one of the sympy code generation
                # routines, or lambdify with 'numexpr' to implement this in a
                # more robust way
                coeff = re.sub("I", "(1.0j)", str(op.coeff.val))
            else:
                raise ValueError(("Invalid value '%s' for `convert_as`, must "
                                  "be one of 'str', 'pyfunc'") % convert_as)
            return [convert_to_qutip(op.term, full_space), coeff]
        else:
            raise ValueError("op cannot be expressed in qutip. It must have "
                             "the structure op = sum_i f_i(t) * op_i")
    else:
        return convert_to_qutip(op, full_space=full_space)