def METHOD_NAME(L, wlist, rhoss, J_ops):
    J_ops = [op.data for op in J_ops]
    rhoss_vec = operator_to_vector(rhoss).data
    N_j_ops = len(J_ops)
    current = np.zeros(N_j_ops)
    noise = np.zeros((N_j_ops, N_j_ops, len(wlist)))
    tr_op = qeye(L.dims[0][0])
    tr_op_vec = operator_to_vector(tr_op)
    Pop = _data.kron(rhoss_vec, tr_op_vec.data.transpose())
    Iop = _data.identity(np.prod(L.dims[0][0])**2)
    Q = _data.sub(Iop, Pop)
    Q_ops = [_data.matmul(Q, _data.matmul(op, rhoss_vec)) for op in J_ops]
    for k, w in enumerate(wlist):
        if w != 0.0:
            L_temp = 1.0j * w * spre(tr_op) + L
        else:
            # At zero frequency some solvers fail for small systems.
            # Adding a small finite frequency of order 1e-15
            # helps prevent the solvers from throwing an exception.
            with CoreOptions(auto_tidyup=False):
                L_temp = 1e-15j * spre(tr_op) + L
        X_rho = [_solve(L_temp.data, op) for op in Q_ops]
        for i, j in product(range(N_j_ops), repeat=2):
            if i == j:
                current[i] = _data.expect_super(J_ops[i], rhoss_vec).real
                noise[j, i, k] = current[i]
            noise[j, i, k] -= (
                _data.expect_super(_data.matmul(J_ops[j], Q), X_rho[i]) +
                _data.expect_super(_data.matmul(J_ops[i], Q), X_rho[j])
            ).real
    return current, noise