fn make_y_bus(baseMVA, bus, branch) {
    // constants
    nb = size(bus, 0);          // number of buses
    nl = size(branch, 0);       // number of lines

    stat = slice(branch, [0], BR_STATUS);                    // ones at in-service branches
    Ys = stat ./ (slice(branch, [0], BR_R) + c(0,1) * slice(branch, [0], BR_X));  // series admittance
    Bc = stat .* slice(branch, [0], BR_B);                           // line charging susceptance
    tap = ones(nl, 1);                              // default tap ratio = 1

    //i = find(branch(:, TAP));                       // indices of non-zero tap ratios
    //tap(i) = slice(branch, [0], TAP);               // assign non-zero tap ratios
    tap = tap .* exp(c(0,1)*pi/180 * slice(branch, [0], SHIFT)); // add phase shifters
    Ytt = Ys + c(0,1) * Bc/2;
    Yff = Ytt ./ (tap .* conj(tap));
    Yft = - Ys ./ conj(tap);
    Ytf = - Ys ./ tap;

    // compute shunt admittance
    // if Psh is the real power consumed by the shunt at V = 1.0 p.u.
    // and Qsh is the reactive power injected by the shunt at V = 1.0 p.u.
    // then Psh - j Qsh = V * conj(Ysh * V) = conj(Ysh) = Gs - j Bs,
    // i.e. Ysh = Psh + j Qsh, so ...
    //Ysh = (slice(bus, [0], GS) + c(0,1) * slice(bus, [0], BS)) / baseMVA; // vector of shunt admittances

    // bus indices
    //f = slice(branch, [0], F_BUS);                           // list of "from" buses
    //t = slice(branch, [0], T_BUS);                           // list of "to" buses

    // for best performance, choose method based on MATLAB vs Octave and size
    if nb < 300 {  // small case
        // build Yf and Yt such that Yf * V is the vector of complex branch currents injected
        // at each branch's "from" bus, and Yt is the same for the "to" bus end
        //i = [1:nl 1:nl]';                           // double set of row indices
        //Yf = sparse(i, [f, t], [Yff, Yft], nl, nb);
        //Yt = sparse(i, [f, t], [Ytf, Ytt], nl, nb);

        // build Ybus
        //Ybus = sparse([f,f,t,t], [f,t,f,t], [Yff,Yft,Ytf,Ytt], nb, nb) + // branch admittances
        //        sparse(1:nb, 1:nb, Ysh, nb, nb);        // shunt admittance
    } else {  // large case running on MATLAB
        // build connection matrices
        //Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      // connection matrix for line & from buses
        //Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      // connection matrix for line & to buses

        // build Yf and Yt such that Yf * V is the vector of complex branch currents injected
        // at each branch's "from" bus, and Yt is the same for the "to" bus end
        //Yf = sparse(1:nl, 1:nl, Yff, nl, nl) * Cf + sparse(1:nl, 1:nl, Yft, nl, nl) * Ct;
        //Yt = sparse(1:nl, 1:nl, Ytf, nl, nl) * Cf + sparse(1:nl, 1:nl, Ytt, nl, nl) * Ct;

        // build Ybus
        //Ybus = Cf' * Yf + Ct' * Yt + ...            // branch admittances
        //        sparse(1:nb, 1:nb, Ysh, nb, nb);    // shunt admittance
    }

    return Ytt;
}

