Based on dSbr_dV.m from MATPOWER by Ray Zimmerman, developed at
        PSERC Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more
        information.

        @return: The branch power flow vectors and the partial derivatives of
                 branch power flow w.r.t voltage magnitude and voltage angle.
        @rtype: tuple

def dSbr_dV(self, Yf, Yt, V, buses=None, branches=None):
        """ Based on dSbr_dV.m from MATPOWER by Ray Zimmerman, developed at
        PSERC Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more
        information.

        @return: The branch power flow vectors and the partial derivatives of
                 branch power flow w.r.t voltage magnitude and voltage angle.
        @rtype: tuple
        """
        buses = self.buses if buses is None else buses
        branches = self.branches if branches is None else branches

        nl = len(branches)
        nb = len(V)
        il = range(nl)
        ib = range(nb)

        f = [l.from_bus._i for l in branches]
        t = [l.to_bus._i for l in branches]

        # Compute currents.
        If = Yf * V
        It = Yt * V

        Vnorm = V / abs(V)

        diagVf = csr_matrix((V[f], (il, il)))
        diagIf = csr_matrix((If, (il, il)))
        diagVt = csr_matrix((V[t], (il, il)))
        diagIt = csr_matrix((It, (il, il)))
        diagV  = csr_matrix((V, (ib, ib)))
        diagVnorm = csr_matrix((Vnorm, (ib, ib)))

        shape = (nl, nb)
        # Partial derivative of S w.r.t voltage phase angle.
        dSf_dVa = 1j * (conj(diagIf) *
            csr_matrix((V[f], (il, f)), shape) - diagVf * conj(Yf * diagV))

        dSt_dVa = 1j * (conj(diagIt) *
            csr_matrix((V[t], (il, t)), shape) - diagVt * conj(Yt * diagV))

        # Partial derivative of S w.r.t. voltage amplitude.
        dSf_dVm = diagVf * conj(Yf * diagVnorm) + conj(diagIf) * \
            csr_matrix((Vnorm[f], (il, f)), shape)

        dSt_dVm = diagVt * conj(Yt * diagVnorm) + conj(diagIt) * \
            csr_matrix((Vnorm[t], (il, t)), shape)

        # Compute power flow vectors.
        Sf = V[f] * conj(If)
        St = V[t] * conj(It)

        return dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St