"""
Pandapower Cross Validation Script

Builds the identical WSCC 9-bus network in pandapower and runs pandapower's
own AC power flow solver, for direct comparison.
"""
import sys
import os
import numpy as np
import pandapower as pp

sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "data"))
sys.path.insert(0, os.path.dirname(__file__))
import network_9bus as net_data
from network_setup import build_case
from newton_raphson import run_newton_raphson

S_BASE = net_data.S_BASE_MVA  # 100 MVA system base
V_BASE_HV = 230.0             # kV, the transmission level base


def build_pandapower_net():
    """
    Construct the WSCC 9-bus network as a pandapower network object,
    converting this project's per unit (100 MVA base) impedance data to
    the physical units (ohms, nF, %)
    """
    net = pp.create_empty_network(sn_mva=S_BASE, f_hz=60.0)

    bus_idx = {}
    for b, data in net_data.BUS_DATA.items():
        bus_idx[b] = pp.create_bus(net, vn_kv=data["base_kv"], name=f"Bus {b}")

    # Slack (external grid) at Bus 1
    pp.create_ext_grid(net, bus_idx[1], vm_pu=net_data.BUS_DATA[1]["v_set"],
                        name="Gen 1 (slack)")

    # PV generators at Buses 2 and 3
    for gbus in [2, 3]:
        gd = net_data.GEN_DATA[gbus]
        pp.create_gen(net, bus_idx[gbus], p_mw=gd["P_mw"],
                      vm_pu=net_data.BUS_DATA[gbus]["v_set"],
                      max_q_mvar=gd["Qmax"], min_q_mvar=gd["Qmin"],
                      name=f"Gen {gbus}")

    # Loads
    for lbus, ld in net_data.LOAD_DATA.items():
        pp.create_load(net, bus_idx[lbus], p_mw=ld["P_mw"], q_mvar=ld["Q_mvar"],
                        name=f"Load {lbus}")

    # Transmission lines: convert per unit R/X/B (100 MVA, 230 kV base) to
    # ohms/nF using a nominal 1 km length so the *total* per unit impedance
    # is preserved as the per km value.
    z_base = V_BASE_HV ** 2 / S_BASE  # ohms
    omega = 2 * np.pi * 60.0

    for (f, t), d in net_data.LINE_DATA.items():
        r_ohm = d["R"] * z_base
        x_ohm = d["X"] * z_base
        b_siemens = d["B"] / z_base
        c_nf = (b_siemens / omega) * 1e9
        i_max_ka = d["rate_mva"] / (np.sqrt(3) * V_BASE_HV)
        pp.create_line_from_parameters(
            net, bus_idx[int(f)], bus_idx[int(t)], length_km=1.0,
            r_ohm_per_km=r_ohm, x_ohm_per_km=x_ohm, c_nf_per_km=c_nf,
            max_i_ka=i_max_ka, name=f"Line {f}-{t}")

    # Generator step-up transformers: R = 0 (negligible), X given on the
    # 100 MVA system base -> converted to vk_percent on the transformer's
    # own MVA rating (sn_mva) as pandapower expects.
    for (f, t), d in net_data.XFMR_DATA.items():
        hv_bus, lv_bus = bus_idx[int(t)], bus_idx[int(f)]  # t=230kV side, f=gen side
        vk_percent = d["X"] * (S_BASE / d["rate_mva"]) * 100.0
        pp.create_transformer_from_parameters(
            net, hv_bus=hv_bus, lv_bus=lv_bus, sn_mva=d["rate_mva"],
            vn_hv_kv=V_BASE_HV, vn_lv_kv=net_data.BUS_DATA[int(f)]["base_kv"],
            vk_percent=vk_percent, vkr_percent=0.0, pfe_kw=0.0, i0_percent=0.0,
            shift_degree=0.0, name=f"T{f}-{t}")

    return net, bus_idx


def main():
    print("Building pandapower network model...")
    net, bus_idx = build_pandapower_net()

    print("Running pandapower AC power flow (Newton-Raphson algorithm)...")
    pp.runpp(net, algorithm="nr", calculate_voltage_angles=True)

    pp_v = np.array([net.res_bus.at[bus_idx[b], "vm_pu"] for b in range(1, 10)])
    pp_theta_deg = np.array([net.res_bus.at[bus_idx[b], "va_degree"] for b in range(1, 10)])

    print("Running this project's custom Newton-Raphson solver...")
    case = build_case(load_scale=1.0)
    res = run_newton_raphson(case["Ybus"], case["bus_types"], case["P_spec"], case["Q_spec"],
                              V0=case["V0"], theta0=case["theta0"])

    print("\nBus  |  V_custom (pu)  |  V_pandapower (pu)  |  Deviation (%)")
    max_dev = 0.0
    for i in range(9):
        dev_pct = 100 * abs(res.V[i] - pp_v[i]) / res.V[i]
        max_dev = max(max_dev, dev_pct)
        print(f"{i + 1:4d} | {res.V[i]:14.5f}  | {pp_v[i]:18.5f}  | {dev_pct:10.6f}")

    print(f"\nMaximum voltage magnitude deviation vs. pandapower: {max_dev:.6f}%")
    print("(Validation threshold for this project: 2%)")
    return max_dev


if __name__ == "__main__":
    main()
