//DSBR_DV   Computes partial derivatives of branch power flows w.r.t. voltage.
//
//   The derivatives can be take with respect to polar or cartesian coordinates
//   of voltage, depending on the 5th argument.
//
//   [DSF_DVA, DSF_DVM, DST_DVA, DST_DVM, SF, ST] = DSBR_DV(BRANCH, YF, YT, V)
//   [DSF_DVA, DSF_DVM, DST_DVA, DST_DVM, SF, ST] = DSBR_DV(BRANCH, YF, YT, V, 0)
//
//   Returns four matrices containing partial derivatives of the complex
//   branch power flows at "from" and "to" ends of each branch w.r.t voltage
//   magnitude and voltage angle, respectively (for all buses).
//
//   [DSF_DVR, DSF_DVI, DST_DVR, DST_DVI, SF, ST] = DSBR_DV(BRANCH, YF, YT, V, 1)
//
//   Returns four matrices containing partial derivatives of the complex
//   branch power flows at "from" and "to" ends of each branch w.r.t real and
//   imaginary parts of voltage, respectively (for all buses).
//
//   If YF is a sparse matrix, the partial derivative matrices will be as well.
//   Optionally returns vectors containing the power flows themselves. The
//   following explains the expressions used to form the matrices:
//
//   If = Yf * V;
//   Sf = diag(Vf) * conj(If) = diag(conj(If)) * Vf
//
//   Polar coordinates:
//     Partials of V, Vf & If w.r.t. voltage angles
//       dV/dVa  = j * diag(V)
//       dVf/dVa = sparse(1:nl, f, j * V(f)) = j * sparse(1:nl, f, V(f))
//       dIf/dVa = Yf * dV/dVa = Yf * j * diag(V)
//
//     Partials of V, Vf & If w.r.t. voltage magnitudes
//       dV/dVm  = diag(V./abs(V))
//       dVf/dVm = sparse(1:nl, f, V(f)./abs(V(f))
//       dIf/dVm = Yf * dV/dVm = Yf * diag(V./abs(V))
//
//     Partials of Sf w.r.t. voltage angles
//       dSf/dVa = diag(Vf) * conj(dIf/dVa)
//                       + diag(conj(If)) * dVf/dVa
//               = diag(Vf) * conj(Yf * j * diag(V))
//                       + conj(diag(If)) * j * sparse(1:nl, f, V(f))
//               = -j * diag(Vf) * conj(Yf * diag(V))
//                       + j * conj(diag(If)) * sparse(1:nl, f, V(f))
//               = j * (conj(diag(If)) * sparse(1:nl, f, V(f))
//                       - diag(Vf) * conj(Yf * diag(V)))
//
//     Partials of Sf w.r.t. voltage magnitudes
//       dSf/dVm = diag(Vf) * conj(dIf/dVm)
//                       + diag(conj(If)) * dVf/dVm
//               = diag(Vf) * conj(Yf * diag(V./abs(V)))
//                       + conj(diag(If)) * sparse(1:nl, f, V(f)./abs(V(f)))
//
//   Cartesian coordinates:
//     Partials of V, Vf & If w.r.t. real part of complex voltage
//       dV/dVr  = diag(ones(n,1))
//       dVf/dVr = Cf
//       dIf/dVr = Yf
//     where Cf is the connection matrix for line & from buses
//
//     Partials of V, Vf & If w.r.t. imaginary part of complex voltage
//       dV/dVi  = j * diag(ones(n,1))
//       dVf/dVi = j * Cf
//       dIf/dVi = j * Yf
//
//     Partials of Sf w.r.t. real part of complex voltage
//       dSf/dVr = conj(diag(If)) * Cf + diag(Vf) * conj(Yf)
//
//     Partials of Sf w.r.t. imaginary part of complex voltage
//       dSf/dVi = j * (conj(diag(If)) * Cf - diag(Vf) * conj(Yf))
//
//   Derivations for "to" bus are similar.
//
//   Examples:
//       [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
//       [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = ...
//           dSbr_dV(branch, Yf, Yt, V);
//       [dSf_dVr, dSf_dVi, dSt_dVr, dSt_dVi, Sf, St] = ...
//           dSbr_dV(branch, Yf, Yt, V, 1);
//  计算支路功率流对电压的偏导数。

fn dSbr_dV(branch, Yf, Yt, V, vcart) {

    // 默认输入参数
    if is_empty(vcart) {
        vcart = 0;      // 默认为极坐标
    }

    // 定义
    f = slice(branch, [0], [F_BUS - 1, F_BUS]) - 1;     // "from" 母线列表(f索引从'0'开始)
    t = slice(branch, [0], [T_BUS - 1, T_BUS]) - 1;     // "to" 母线列表(t索引从'0'开始)
    nl = size(f, 0);
    nb = size(V, 0);

    // 计算中间值
    Yfc = conj(Yf);
    Ytc = conj(Yt);
    Vc = conj(V);
    Ifc = Yfc * Vc;     // "from" 电流的共轭
    Itc = Ytc * Vc;     // "to" 电流的共轭

    diagVf  = diag(get_multi(V, f));
    diagVt  = diag(get_multi(V, t));
    diagIfc = diag(Ifc);
    diagItc = diag(Itc);
    if vcart == 0 {  // 极坐标
        Vnorm       = V ./ abs(V);
        diagVc      = diag(Vc);
        diagVnorm   = diag(Vnorm);
        CVf  = full(sparse(range(0, nl), f, get_multi(V, f), nl, nb));
        CVnf = full(sparse(range(0, nl), f, get_multi(Vnorm, f), nl, nb));
        CVt  = full(sparse(range(0, nl), t, get_multi(V, t), nl, nb));
        CVnt = full(sparse(range(0, nl), t, get_multi(Vnorm, t), nl, nb));
    }

    if vcart {  // 直角坐标
        Cf = sparse(range(0, nl), f, ones(nl, 1), nl, nb);      // 连接矩阵 for line & from buses
        Ct = sparse(range(0, nl), t, ones(nl, 1), nl, nb);      // 连接矩阵 for line & to buses
        Af = diagIfc * Cf;
        Bf = diagVf * Yfc;
        At = diagItc * Ct;
        Bt = diagVt * Ytc;

        dSf_dV1 = Af + Bf;          // dSf_dVr
        dSf_dV2 = c(0, 1) * (Af - Bf);   // dSf_dVi
        dSt_dV1 = At + Bt;          // dSt_dVr
        dSt_dV2 = c(0, 1) * (At - Bt);   // dSt_dVi
    } else {  // 极坐标
        dSf_dV1 = c(0, 1) * (diagIfc * CVf - diagVf * Yfc * diagVc);     // dSf_dVa
        dSf_dV2 = diagVf * conj(Yf * diagVnorm) + diagIfc * CVnf;   // dSf_dVm
        dSt_dV1 = c(0, 1) * (diagItc * CVt - diagVt * Ytc * diagVc);     // dSt_dVa
        dSt_dV2 = diagVt * conj(Yt * diagVnorm) + diagItc * CVnt;   // dSt_dVm
    }

    // 返回功率流本身
    Sf = get_multi(V, f)' .* Ifc;
    St = get_multi(V, t)' .* Itc;

    /*if is_sparse(Yf) {             // 稀疏矩阵版本 (如果Yf是稀疏矩阵)
        diagVf  = sparse(range(0, nl), range(0, nl), get_multi(V, f), nl, nl);
        diagVt  = sparse(range(0, nl), range(0, nl), get_multi(V, t), nl, nl);
        diagIfc = sparse(range(0, nl), range(0, nl), Ifc, nl, nl);
        diagItc = sparse(range(0, nl), range(0, nl), Itc, nl, nl);
        if vcart == 0 {  // 极坐标
            Vnorm       = V ./ abs(V);
            diagVc      = sparse(range(0, nb), range(0, nb), Vc, nb, nb);
            diagVnorm   = sparse(range(0, nb), range(0, nb), Vnorm, nb, nb);
            CVf  = sparse(range(0, nl), f, get_multi(V, f), nl, nb);
            CVnf = sparse(range(0, nl), f, get_multi(get_multi(Vnorm, f), range(0, size(f, 0))), nl, nb);
            CVt  = sparse(range(0, nl), t, get_multi(V, t), nl, nb);
            CVnt = sparse(range(0, nl), t, get_multi(get_multi(Vnorm, t), range(0, size(t, 0))), nl, nb);
        }
    } else {                        // 稠密矩阵版本
        diagVf  = diag(get_multi(V, f));
        diagVt  = diag(get_multi(V, t));
        diagIfc = diag(Ifc);
        diagItc = diag(Itc);
        if vcart == 0 {  // 极坐标
            Vnorm       = V ./ abs(V);
            diagVc      = diag(Vc);
            diagVnorm   = diag(Vnorm);
            CVf  = full(sparse(range(0, nl), f, get_multi(V, f), nl, nb));
            CVnf = full(sparse(range(0, nl), f, get_multi(Vnorm, f), nl, nb));
            CVt  = full(sparse(range(0, nl), t, get_multi(V, t), nl, nb));
            CVnt = full(sparse(range(0, nl), t, get_multi(Vnorm, t), nl, nb));
        }
    }
    */

    return (dSf_dV1, dSf_dV2, dSt_dV1, dSt_dV2, Sf, St);
}
