//NEWTONPF  Solves power flow using full Newton's method (power/polar)
//   [V, CONVERGED, I] = NEWTONPF(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
//
//   Solves for bus voltages using a full Newton-Raphson method, using nodal
//   power balance equations and polar coordinate representation of
//   voltages, given the following inputs:
//       YBUS  - full system admittance matrix (for all buses)
//       SBUS  - handle to function that returns the complex bus power
//               injection vector (for all buses), given the bus voltage
//               magnitude vector (for all buses)
//       V0    - initial vector of complex bus voltages
//       REF   - bus index of reference bus (voltage ang reference & gen slack)
//       PV    - vector of bus indices for PV buses
//       PQ    - vector of bus indices for PQ buses
//       MPOPT - (optional) MATPOWER option struct, used to set the
//               termination tolerance, maximum number of iterations, and
//               output options (see MPOPTION for details).
//
//   The bus voltage vector contains the set point for generator
//   (including ref bus) buses, and the reference angle of the swing
//   bus, as well as an initial guess for remaining magnitudes and
//   angles.
//
//   Returns the final complex voltages, a flag which indicates whether it
//   converged or not, and the number of iterations performed.
//
//   See also RUNPF, NEWTONPF_S_CART, NEWTONPF_I_POLAR, NEWTONPF_I_CART.

//   MATPOWER
//   Copyright (c) 1996-2019, Power Systems Engineering Research Center (PSERC)
//   by Ray Zimmerman, PSERC Cornell
//
//   This file is part of MATPOWER.
//   Covered by the 3-clause BSD License (see LICENSE file for details).
//   See https://matpower.org for more info.

fn newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) {
    // options
    tol         = get(mpopt, 0);
    max_it      = get(mpopt, 1);

    // initialize
    converged = 0;
    i = 0;
    V = V0;
    Va = angle(V);
    Vm = abs(V);

    // set up indexing for updating V
    nb = length(V);
    npv = length(pv);
    npq = length(pq);
    j1 = 0;         j2 = npv;       // j1:j2 - V angle of pv buses
    j3 = j2;    j4 = j2 + npq;      // j3:j4 - V angle of pq buses
    j5 = j4;    j6 = j4 + npq;      // j5:j6 - V mag of pq buses
    pv_pq_i = horzcat(pv, pq);

    // evaluate F(x0)
    mis = V .* conj(Ybus * V) - Sbus;
    F = vertcat(real(get_multi(mis, pv_pq_i))', imag(get_multi(mis, pq))');

    // check tolerance
    normF = norm_max(F);
    if normF < tol {
        converged = 1;
    }

    // do Newton iterations
    while ~~converged && i < max_it {
        // update iteration counter
        i = i + 1;
        // evaluate Jacobian 因迭代中V变化，因此jac需要从dSbus_dV重新计算
        [dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V, 0);

        j11 = real(select(dSbus_dVa, pv_pq_i, pv_pq_i));
        j12 = real(select(dSbus_dVm, pv_pq_i, pq));
        j21 = imag(select(dSbus_dVa, pq, pv_pq_i));
        j22 = imag(select(dSbus_dVm, pq, pq));
        J = blockcat(2, 2, j11, j12, j21, j22);
        // compute update step
        negF = 0 - F;
        dx = J \ negF;
        // update voltage vector
        if npv {
            Va = set2(Va, pv, slice(dx, [j1,j2], [0]));
        }
        if npq {
            Va = set2(Va, pq, slice(dx, [j3,j4], [0]));
            Vm = set2(Vm, pq, slice(dx, [j5,j6], [0]));
        }
        V = Vm .* exp(1j * Va);
        Vm = abs(V);            // update Vm and Va again in case
        Va = angle(V);          // we wrapped around with a negative Vm

        mis = V .* conj(Ybus * V) - Sbus;
        F = vertcat(real(get_multi(mis, pv_pq_i))', imag(get_multi(mis, pq))');
        // check for convergence
        normF = norm_max(F);
        if normF < tol {
            converged = 1;
        }
    }
    return (V, converged, i);
}