// runpf - Runs an AC Newton power flow.
//
//   [V, BUS, GEN, BRANCH, SUCCESS, ITERATIONS] = RUNPF()
//   [V, BUS, GEN, BRANCH, SUCCESS, ITERATIONS] =
//       RUNPF_CASE(BASEMVA, BUS, GEN, BRANCH)
//
//   Uses the BASEMVA, BUS, GEN and BRANCH values in the script context. This
//   rspower implementation intentionally supports only the AC power-balance,
//   polar-voltage Newton method. It performs external/internal bus-number
//   conversion and returns solved matrices using the original bus numbers.

fn runpf_case(case_baseMVA, case_bus, case_gen, case_branch) {
    baseMVA = case_baseMVA;
    bus = case_bus;
    gen = case_gen;
    branch = case_branch;

    pf_branch = branch;
    if size(pf_branch, 1) < QT {
        pf_branch = horzcat(pf_branch, zeros(size(pf_branch, 0), QT - size(pf_branch, 1)));
    }

    [i2e, pf_bus, pf_gen, pf_branch] = ext2int(bus, gen, pf_branch);
    [ref, pv, pq] = bustypes(pf_bus, pf_gen);
    ref = ref';
    pv = pv';
    pq = pq';

    nb = size(pf_bus, 0);

    // generator info
    on = find(slice(pf_gen, [0], [GEN_STATUS-1, GEN_STATUS]) > 0);
    gbus = get_multi(slice(pf_gen, [0], [GEN_BUS-1, GEN_BUS]), on) - 1;

    // initial state
    V0 = slice(pf_bus, [0], [VM-1,VM]) .* exp(1j * pi/180 * slice(pf_bus, [0], [VA-1,VA]));
    vcb = set(ones(size(V0)), pq, zeros(size(pq)));
    k = find(get_multi(vcb, gbus));
    on_k = get_multi(on, k);
    gen_vg = get_multi(slice(pf_gen, [0], [VG-1,VG]), on_k);
    gbus_k = get_multi(gbus, k);
    voltage_g = get_multi(V0, gbus_k);
    V_init = set(V0, gbus_k, gen_vg ./ abs(voltage_g) .* voltage_g);

    // build admittance matrices
    [Ybus, Yf, Yt] = makeYbus(baseMVA, pf_bus, pf_branch);
    Sbus = makeSbus(baseMVA, pf_bus, pf_gen);
    [V, success, iterations] = newtonpf(Ybus, Sbus, V_init, ref, pv, pq, [1e-6, 100]);

    [pf_bus, pf_gen, pf_branch] = pfsoln(
        baseMVA, pf_bus, pf_gen, pf_branch,
        Ybus, Yf, Yt, V, ref, pv, pq, []
    );
    [result_bus, result_gen, result_branch] = int2ext(i2e, pf_bus, pf_gen, pf_branch);

    return (V, result_bus, result_gen, result_branch, success, iterations);
}

// Entry point for scripts that load a case into globals.
fn runpf() {
    [V, result_bus, result_gen, result_branch, success, iterations] =
        runpf_case(baseMVA, bus, gen, branch);
    return (V, result_bus, result_gen, result_branch, success, iterations);
}