Computes the branch power flow vector and the partial derivative of
        branch power flow w.r.t voltage.

def dSbr_dV(Yf, Yt, V, buses, branches):
    """ Computes the branch power flow vector and the partial derivative of
        branch power flow w.r.t voltage.
    """
    nl = len(branches)
    nb = len(V)

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

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

    Vnorm = div(V, abs(V))

    diagVf = spdiag(V[f])
    diagIf = spdiag(If)
    diagVt = spdiag(V[t])
    diagIt = spdiag(It)
    diagV = spdiag(V)
    diagVnorm = spdiag(Vnorm)

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

    dSt_dVa = 1j * (conj(diagIt) *
        spmatrix(V[t], ibr, t, size) - diagVt * conj(Yt * diagV))

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

    dSt_dVm = diagVt * conj(Yt * diagVnorm) + conj(diagIt) * \
        spmatrix(Vnorm[t], ibr, t, size)

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

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