# -*- coding: utf-8 -*-

# Copyright (c) 2016-2022 by University of Kassel and Fraunhofer Institute for Energy Economics
# and Energy System Technology (IEE), Kassel. All rights reserved.

import math
import numpy as np
from itertools import product

import pandapower.auxiliary as aux
import pandas as pd
from pandapower.build_bus import _build_bus_ppc
from pandapower.build_gen import _build_gen_ppc
#from pandapower.pd2ppc import _ppc2ppci, _init_ppc
from pandapower.pypower.idx_brch import BR_B, BR_R, BR_X, F_BUS, T_BUS, branch_cols, BR_STATUS, SHIFT, TAP
from pandapower.pypower.idx_bus import BASE_KV, BS, GS, BUS_TYPE
from pandapower.shortcircuit.idx_brch import branch_cols_sc
from pandapower.shortcircuit.idx_bus import C_MAX, C_MIN
from pandapower.build_branch import _calc_tap_from_dataframe, _transformer_correction_factor, _calc_nominal_ratio_from_dataframe,\
     get_trafo_values, _trafo_df_from_trafo3w, _calc_branch_values_from_trafo_df, _calc_switch_parameter
from pandapower.build_branch import _switch_branches, _branches_with_oos_buses, _initialize_branch_lookup, _end_temperature_correction_factor
from pandapower.pd2ppc import _ppc2ppci, _init_ppc


def _pd2ppc_zero(net, k_st, sequence=0):
    """
    Builds the ppc data structure for zero impedance system. Includes the impedance values of
    lines and transformers, but no load or generation data.

    For short-circuit calculation, the short-circuit impedance of external grids is also considered.
    """
    # select elements in service (time consuming, so we do it once)
    net["_is_elements"] = aux._select_is_elements_numba(net, sequence=sequence)

    ppc = _init_ppc(net, sequence)

    _build_bus_ppc(net, ppc)
    _build_gen_ppc(net, ppc)
    _add_gen_sc_impedance_zero(net, ppc)
    _add_ext_grid_sc_impedance_zero(net, ppc)
    _build_branch_ppc_zero(net, ppc, k_st)

    # adds auxilary buses for open switches at branches
    _switch_branches(net, ppc)

    # add auxilary buses for out of service buses at in service lines.
    # Also sets lines out of service if they are connected to two out of service buses
    _branches_with_oos_buses(net, ppc)
    if hasattr(net, "_isolated_buses"):
        ppc["bus"][net._isolated_buses, BUS_TYPE] = 4.

    # generates "internal" ppci format (for powerflow calc) from "external" ppc format and updates the bus lookup
    # Note: Also reorders buses and gens in ppc
    ppci = _ppc2ppci(ppc, net)
    #net._ppc0 = ppc    <--Obsolete. now covered in _init_ppc
    return ppc, ppci

def _build_branch_ppc_zero(net, ppc, k_st=None):
    """
    Takes the empty ppc network and fills it with the zero imepdance branch values. The branch
    datatype will be np.complex 128 afterwards.

    .. note:: The order of branches in the ppc is:
            1. Lines
            2. Transformers

    **INPUT**:
        **net** -The pandapower format network

        **ppc** - The PYPOWER format network to fill in values

    """
    length = _initialize_branch_lookup(net)
    lookup = net._pd2ppc_lookups["branch"]
    mode = net._options["mode"]
    ppc["branch"] = np.zeros(shape=(length, branch_cols), dtype=np.complex128)
    if mode == "sc":
        branch_sc = np.empty(shape=(length, branch_cols_sc), dtype=float)
        branch_sc.fill(np.nan)
        ppc["branch"] = np.hstack((ppc["branch"], branch_sc ))
    ppc["branch"][:, :13] = np.array([0, 0, 0, 0, 0, 250, 250, 250, 1, 0, 1, -360, 360])

    _add_line_sc_impedance_zero(net, ppc)
    _add_trafo_sc_impedance_zero(net, ppc, k_st=k_st)
    if "switch" in lookup:
        _calc_switch_parameter(net, ppc)
    if mode == "sc":
        _add_trafo3w_sc_impedance_zero(net, ppc)
    else:
        if "trafo3w" in lookup:
            raise NotImplementedError("Three winding transformers are not implemented for unbalanced calculations")


def _add_trafo_sc_impedance_zero(net, ppc, trafo_df=None, k_st=None):
    if trafo_df is None:
        trafo_df = net["trafo"]
    if k_st is None:
        k_st = np.ones(len(ppc['branch']))
    if "xn_ohm" not in trafo_df.columns:
        trafo_df["xn_ohm"] = 0.
    branch_lookup = net["_pd2ppc_lookups"]["branch"]
    if not "trafo" in branch_lookup:
        return
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
    mode = net["_options"]["mode"]
    trafo_model = net["_options"]["trafo_model"]
    f, t = branch_lookup["trafo"]
    trafo_df["_ppc_idx"] = range(f, t)
    trafo_df['k_st'] = k_st[trafo_df["_ppc_idx"].values].real
    bus_lookup = net["_pd2ppc_lookups"]["bus"]

    hv_bus = get_trafo_values(trafo_df, "hv_bus").astype(int)
    lv_bus = get_trafo_values(trafo_df, "lv_bus").astype(int)
    in_service = get_trafo_values(trafo_df, "in_service").astype(int)
    ppc["branch"][f:t, F_BUS] = bus_lookup[hv_bus]
    ppc["branch"][f:t, T_BUS] = bus_lookup[lv_bus]
    buses_all, gs_all, bs_all = np.array([], dtype=int), np.array([]), \
                                np.array([])
    BIG_NUMBER = 1e20 * ppc["baseMVA"]
    if mode == "sc":
        # Should be considered as connected for all in_service branches
        ppc["branch"][f:t, BR_X] = BIG_NUMBER
        ppc["branch"][f:t, BR_R] = BIG_NUMBER
        ppc["branch"][f:t, BR_B] = 0
        ppc["branch"][f:t, BR_STATUS] = in_service
    else:
        ppc["branch"][f:t, BR_STATUS] = 0

    if not "vector_group" in trafo_df:
        raise ValueError("Vector Group of transformer needs to be specified for zero \
                         sequence modelling \n Try : net.trafo[\"vector_group\"] = 'Dyn'" )

    for vector_group, trafos in trafo_df.groupby("vector_group"):
        # TODO Roman: check this/expand this
        ppc_idx = trafos["_ppc_idx"].values.astype(int)

        if vector_group.lower() in ["yy", "yd", "dy", "dd"]:
            continue

        vk_percent = trafos["vk_percent"].values.astype(float)
        vkr_percent = trafos["vkr_percent"].values.astype(float)
        sn_trafo_mva = trafos["sn_mva"].values.astype(float)
        # Just put pos seq parameter if zero seq parameter is zero
        if not "vk0_percent" in trafos:
            raise ValueError("Short circuit voltage of transformer Vk0 needs to be specified for zero \
                             sequence modelling \n Try : net.trafo[\"vk0_percent\"] = net.trafo[\"vk_percent\"]" )
        vk0_percent = trafos["vk0_percent"].values.astype(float) if \
            trafos["vk0_percent"].values.astype(float).all() != 0. else \
            trafos["vk_percent"].values.astype(float)
        # Just put pos seq parameter if zero seq parameter is zero
        if not "vkr0_percent" in trafos:
            raise ValueError("Real part of short circuit voltage Vk0(Real) needs to be specified for transformer \
                             modelling \n Try : net.trafo[\"vkr0_percent\"] = net.trafo[\"vkr_percent\"]" )
        vkr0_percent = trafos["vkr0_percent"].values.astype(float) if \
            trafos["vkr0_percent"].values.astype(float).all() != 0. else \
            trafos["vkr_percent"].values.astype(float)
        lv_buses = trafos["lv_bus"].values.astype(int)
        hv_buses = trafos["hv_bus"].values.astype(int)
        lv_buses_ppc = bus_lookup[lv_buses]
        hv_buses_ppc = bus_lookup[hv_buses]
        if not "mag0_percent" in trafos:
            # For Shell Type transformers vk0 = vk * 1
            #        and mag0_percent = 10 ... 100  Zm0/ Zsc0
            # --pg 50 DigSilent Power Factory Transformer manual
            raise ValueError("Magnetizing impedance to vk0 ratio needs to be specified for transformer \
                             modelling  \n Try : net.trafo[\"mag0_percent\"] = 100" )
        mag0_ratio = trafos.mag0_percent.values.astype(float)
        if not "mag0_rx" in trafos:
            raise ValueError("Magnetizing impedance R/X ratio needs to be specified for transformer \
                             modelling \n Try : net.trafo[\"mag0_rx\"] = 0 " )
        mag0_rx = trafos["mag0_rx"].values.astype(float)
        if not "si0_hv_partial" in trafos:
            raise ValueError("Zero sequence short circuit impedance partition towards HV side needs to be specified for transformer \
                             modelling \n Try : net.trafo[\"si0_hv_partial\"] = 0.9 " )
        si0_hv_partial = trafos.si0_hv_partial.values.astype(float)
        parallel = trafos.parallel.values.astype(float)
        if "power_station_unit" in trafos.columns:
            power_station_unit = trafos.power_station_unit.fillna(False).astype(bool)
        else:
            power_station_unit = np.zeros(len(trafos), dtype=bool)
        in_service = trafos["in_service"].astype(int)

        ppc["branch"][ppc_idx, F_BUS] = hv_buses_ppc
        ppc["branch"][ppc_idx, T_BUS] = lv_buses_ppc

        vn_trafo_hv, vn_trafo_lv, shift = _calc_tap_from_dataframe(net, trafos)
        vn_bus_lv = ppc["bus"][lv_buses_ppc, BASE_KV]
        vn_bus_hv = ppc["bus"][hv_buses_ppc, BASE_KV]
        ratio = _calc_nominal_ratio_from_dataframe(ppc, trafos, vn_trafo_hv, \
                                                   vn_trafo_lv, bus_lookup)
        ppc["branch"][ppc_idx, TAP] = ratio
        ppc["branch"][ppc_idx, SHIFT] = shift

        # zero seq. transformer impedance
        tap_lv = np.square(vn_trafo_lv / vn_bus_lv) * net.sn_mva
        tap_hv = np.square(vn_trafo_hv / vn_bus_hv) * net.sn_mva
        if mode == 'pf_3ph':
            if vector_group.lower() not in ["ynyn", "dyn", "yzn"]:
                raise NotImplementedError("Calculation of 3-phase power flow is only implemented for the transformer "
                                          "vector groups 'YNyn', 'Dyn', 'Yzn'")
            # =============================================================================
            #     Changing base from transformer base to Network base to get Zpu(Net)
            #     Zbase = (kV).squared/S_mva
            #     Zpu(Net)={Zpu(trafo) * Zb(trafo)} / {Zb(Net)}
            #        Note:
            #             Network base voltage is Line-Neutral voltage in each phase
            #             Line-Neutral voltage= Line-Line Voltage(vn_bus_lv) divided by sq.root(3)
            # =============================================================================
            tap_lv = np.square(vn_trafo_lv / vn_bus_lv) * (3 * net.sn_mva)
            tap_hv = np.square(vn_trafo_hv / vn_bus_hv) * (3 * net.sn_mva)

        tap_corr = tap_hv if vector_group.lower() in ("ynd", "yny") else tap_lv
        z_sc = vk0_percent / 100. / sn_trafo_mva * tap_corr
        r_sc = vkr0_percent / 100. / sn_trafo_mva * tap_corr
        z_sc = z_sc.astype(float)
        r_sc = r_sc.astype(float)
        x_sc = np.sign(z_sc) * np.sqrt(z_sc ** 2 - r_sc ** 2)
        # TODO: This equation needs to be checked!
        # z0_k = (r_sc + x_sc * 1j) / parallel  * max(1, ratio) **2
        # z0_k = (r_sc + x_sc * 1j) / parallel * vn_trafo_hv / vn_bus_hv
        # z0_k = (r_sc + x_sc * 1j) / parallel * tap_hv
        z0_k = (r_sc + x_sc * 1j) / parallel
        z_n_ohm = trafos["xn_ohm"].fillna(0).values
        k_st_tr = trafos["k_st"].fillna(1).values

        if mode == "sc":# or trafo_model == "pi":
            cmax = net._ppc["bus"][lv_buses_ppc, C_MAX]
            kt = _transformer_correction_factor(trafos, vk_percent, vkr_percent, sn_trafo_mva, cmax)
            z0_k *= kt

            # different formula must be applied for power station unit transformers:
            # z_0THV is for power station block unit transformer -> page 20 of IEC60909-4:2021 (example 4.4.2):
            # todo: check if sn_mva must be included here?
            vkx0_percent = np.sqrt(np.square(vk0_percent) - np.square(vkr0_percent))
            z_0THV = (vkr0_percent / 100 + 1j * vkx0_percent / 100) * (np.square(vn_trafo_hv) / sn_trafo_mva) / parallel
            z0_k_psu = (z_0THV * k_st_tr + 3j * z_n_ohm) / ((vn_bus_hv ** 2) / net.sn_mva)
            z0_k = np.where(power_station_unit, z0_k_psu, z0_k)

        y0_k = 1 / z0_k  # adding admittance for "pi" model
        # y0_k = 1 / (z0_k * k_st_tr + 3j * z_n_ohm)  # adding admittance for "pi" model

        # =============================================================================
        #       Transformer magnetising impedance for zero sequence
        # =============================================================================
        z_m = z_sc * mag0_ratio
        x_m = z_m / np.sqrt(mag0_rx ** 2 + 1)
        r_m = x_m * mag0_rx
        r0_trafo_mag = r_m / parallel
        x0_trafo_mag = x_m / parallel
        z0_mag = r0_trafo_mag + x0_trafo_mag * 1j
        # =============================================================================
        #         Star - Delta conversion ( T model to Pi Model)
        #      ----------- |__zc=ZAB__|-----------------
        #            _|                   _|
        #     za=ZAN|_|                  |_| zb=ZBN
        #            |                    |
        # =============================================================================
        z1 = si0_hv_partial * z0_k
        z2 = (1 - si0_hv_partial) * z0_k
        z3 = z0_mag
        z_temp = z1 * z2 + z2 * z3 + z1 * z3
        za = z_temp / z2
#        za = z_temp / (z2+z3)
        zb = z_temp / z1
#        zb = z_temp / (z1+z3)
        zc = z_temp / z3 # ZAB  Transfer impedance
#        zc = z_temp / (z1+z2)  # ZAB  Transfer impedance
        YAB = 1 / zc.astype(complex)
        YAN = 1 / za.astype(complex)
        YBN = 1 / zb.astype(complex)

#        YAB_AN = (zc + za) /(zc * za).astype(complex)  # Series conn YAB and YAN
#        YAB_BN = (zc + zb) / (zc * zb).astype(complex)  # Series conn YAB and YBN

        YAB_AN = 1 / (zc + za).astype(complex)  # Series conn YAB and YAN
        YAB_BN = 1 / (zc + zb).astype(complex)  # Series conn YAB and YBN

        # y0_k = 1 / z0_k #adding admittance for "pi" model
        if vector_group.lower() == "dyn":
            buses_all = np.hstack([buses_all, lv_buses_ppc])
            if trafo_model == "pi":
                y = y0_k * ppc["baseMVA"] # pi model
            else:
                y = (YAB + YBN).astype(complex) * ppc["baseMVA"]  # T model
            gs_all = np.hstack([gs_all, y.real * in_service])
            bs_all = np.hstack([bs_all, y.imag * in_service])

        elif vector_group.lower() == "ynd":
            buses_all = np.hstack([buses_all, hv_buses_ppc])
            if trafo_model == "pi":
                y = y0_k * ppc["baseMVA"] # pi model
                # y = 1/0.99598 * 1 / (1/(y0_k * ppc["baseMVA"]) + 1/0.99598 * (1j * 3 * 22 /( (110 ** 2) / 1))) # pi model
                # y = 1/0.99598 * 1 / (1/(y0_k * ppc["baseMVA"]) + 1/0.99598 * (1j * 3 * 22 /( (110 ** 2) / 1))) # pi model

                # z0_k_k = z0_k * 0.99598 + 1j * 3 * 22 /( (110 ** 2) / 1)
                # print(z0_k_k)
                # y = 1 / z0_k_k # pi model
            else:
                y = (YAB_BN + YAN).astype(complex) * ppc["baseMVA"] #T model
            gs_all = np.hstack([gs_all, y.real * in_service])
            bs_all = np.hstack([bs_all, y.imag * in_service])

        elif vector_group.lower() == "yyn":
            buses_all = np.hstack([buses_all, lv_buses_ppc])
            if trafo_model == "pi":
                y = 1/(z0_mag+z0_k).astype(complex) * ppc["baseMVA"] #pi model
            else:
#                y = (YAB_AN + YBN).astype(complex) #T model
                y = (YAB + YAB_BN + YBN).astype(complex) * ppc["baseMVA"] # T model

            gs_all = np.hstack([gs_all, y.real * in_service])
            bs_all = np.hstack([bs_all, y.imag * in_service])

        elif vector_group.lower() == "ynyn":
            ppc["branch"][ppc_idx, BR_STATUS] = in_service
            # zc = ZAB
            ppc["branch"][ppc_idx, BR_R] = zc.real
            ppc["branch"][ppc_idx, BR_X] = zc.imag

            buses_all = np.hstack([buses_all, hv_buses_ppc])
            gs_all = np.hstack([gs_all, YAN.real * in_service * ppc["baseMVA"] * tap_lv / tap_hv])
            bs_all = np.hstack([bs_all, YAN.imag * in_service * ppc["baseMVA"] * tap_lv / tap_hv])

            buses_all = np.hstack([buses_all, lv_buses_ppc])
            gs_all = np.hstack([gs_all, YBN.real * in_service * ppc["baseMVA"]])
            bs_all = np.hstack([bs_all, YBN.imag * in_service * ppc["baseMVA"]])

        elif vector_group.lower() == "yny":
            buses_all = np.hstack([buses_all, hv_buses_ppc])
            if trafo_model == "pi":
                y = 1/(z0_mag+z0_k).astype(complex) * ppc["baseMVA"]#pi model
            else:
                y = (YAB_BN + YAN).astype(complex) * ppc["baseMVA"]  #T model
            gs_all = np.hstack([gs_all, y.real * in_service])
            bs_all = np.hstack([bs_all, y.imag * in_service])

        elif vector_group.lower() == "yzn":
            buses_all = np.hstack([buses_all, lv_buses_ppc])
            #            y = 1/(z0_mag+z0_k).astype(complex)* int(ppc["baseMVA"])#T model
            #            y= (za+zb+zc)/((za+zc)*zb).astype(complex)* int(ppc["baseMVA"])#pi model
            y = (YAB_AN + YBN).astype(complex) * ppc["baseMVA"] ** 2  #T model # why sn_mva squared here?
            gs_all = np.hstack([gs_all, (1.1547) * y.real * in_service])  # what's the 1.1547 value?
            bs_all = np.hstack([bs_all, (1.1547) * y.imag * in_service])

        elif vector_group[-1].isdigit():
            raise ValueError("Unknown transformer vector group %s -\
                             please specify vector group without \
                             phase shift number. Phase shift can be \
                             specified in net.trafo.shift_degree" % vector_group)
        else:
            raise ValueError("Transformer vector group %s is unknown\
                    / not implemented for three phase load flow" % vector_group)

    buses, gs, bs = aux._sum_by_group(buses_all, gs_all, bs_all)
    ppc["bus"][buses, GS] += gs
    ppc["bus"][buses, BS] += bs
    del net.trafo["_ppc_idx"]


def _add_gen_sc_impedance_zero(net, ppc):
    mode = net["_options"]["mode"]
    if mode == 'pf_3ph':
        return

    eg = net["gen"][net._is_elements["gen"]]
    if len(eg) == 0:
        return
    eg_buses = eg.bus.values
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
    eg_buses_ppc = bus_lookup[eg_buses]

    y0_gen = 1 / (1e3 + 1e3*1j)
    # buses, gs, bs = aux._sum_by_group(eg_buses_ppc, y0_gen.real, y0_gen.imag)
    ppc["bus"][eg_buses_ppc, GS] += y0_gen.real
    ppc["bus"][eg_buses_ppc, BS] += y0_gen.imag

def _add_ext_grid_sc_impedance_zero(net, ppc):
    mode = net["_options"]["mode"]

    if mode == "sc":
        case = net._options["case"]
    else:
        case = "max"
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
    eg = net["ext_grid"][net._is_elements["ext_grid"]]
    if len(eg) == 0:
        return
    eg_buses = eg.bus.values
    eg_buses_ppc = bus_lookup[eg_buses]

    if mode == "sc":
        c = ppc["bus"][eg_buses_ppc, C_MAX] if case == "max" else ppc["bus"][eg_buses_ppc, C_MIN]
    else:
        c = 1.1
    if not "s_sc_%s_mva" % case in eg:
        raise ValueError("short circuit apparent power s_sc_%s_mva needs to be specified for "% case +
                         "external grid" )
    s_sc = eg["s_sc_%s_mva" % case].values
    if not "rx_%s" % case in eg:
        raise ValueError("short circuit R/X rate rx_%s needs to be specified for external grid" %
                         case)
    rx = eg["rx_%s" % case].values

    z_grid = c / s_sc
    if mode == 'pf_3ph':
        z_grid = c / (s_sc/3)
    x_grid = z_grid / np.sqrt(rx ** 2 + 1)
    r_grid = rx * x_grid
    eg["r"] = r_grid
    eg["x"] = x_grid

    # ext_grid zero sequence impedance
    if case == "max":
        x0_grid = net.ext_grid["x0x_%s" % case].values * x_grid
        r0_grid = net.ext_grid["r0x0_%s" % case].values * x0_grid
    elif case == "min":
        x0_grid = net.ext_grid["x0x_%s" % case].values * x_grid
        r0_grid = net.ext_grid["r0x0_%s" % case].values * x0_grid
    y0_grid = 1 / (r0_grid + x0_grid*1j)

    buses, gs, bs = aux._sum_by_group(eg_buses_ppc, y0_grid.real, y0_grid.imag)
    ppc["bus"][buses, GS] = gs
    ppc["bus"][buses, BS] = bs


def _add_line_sc_impedance_zero(net, ppc):
    branch_lookup = net["_pd2ppc_lookups"]["branch"]
    mode = net["_options"]["mode"]
    if not "line" in branch_lookup:
        return
    line = net["line"]
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
    length = line["length_km"].values
    parallel = line["parallel"].values

    fb = bus_lookup[line["from_bus"].values]
    tb = bus_lookup[line["to_bus"].values]
    baseR = np.square(ppc["bus"][fb, BASE_KV]) / net.sn_mva
    if mode == 'pf_3ph':
        baseR = np.square(ppc["bus"][fb, BASE_KV]) / (3*net.sn_mva)
    f, t = branch_lookup["line"]
    # line zero sequence impedance
    ppc["branch"][f:t, F_BUS] = fb
    ppc["branch"][f:t, T_BUS] = tb
    ppc["branch"][f:t, BR_R] = line["r0_ohm_per_km"].values * length / baseR / parallel
    if mode == "sc":
        # temperature correction
        if net["_options"]["case"] == "min":
            ppc["branch"][f:t, BR_R] *= _end_temperature_correction_factor(net, short_circuit=True)
    ppc["branch"][f:t, BR_X] = line["x0_ohm_per_km"].values * length / baseR / parallel
    ppc["branch"][f:t, BR_B] = (2 * net["f_hz"] * math.pi * line["c0_nf_per_km"].values * 1e-9 * baseR * length * parallel)
    ppc["branch"][f:t, BR_STATUS] = line["in_service"].astype(int)


def _add_trafo3w_sc_impedance_zero(net, ppc):
    # TODO Roman: check this/expand this
    branch_lookup = net["_pd2ppc_lookups"]["branch"]
    if not "trafo3w" in branch_lookup:
        return
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
    branch = ppc["branch"]
    f, t = net["_pd2ppc_lookups"]["branch"]["trafo3w"]
    trafo_df = _trafo_df_from_trafo3w(net, sequence=0)
    hv_bus = get_trafo_values(trafo_df, "hv_bus").astype(int)
    lv_bus = get_trafo_values(trafo_df, "lv_bus").astype(int)
    in_service = get_trafo_values(trafo_df, "in_service").astype(int)
    branch[f:t, F_BUS] = bus_lookup[hv_bus]
    branch[f:t, T_BUS] = bus_lookup[lv_bus]

    r, x, _, ratio, shift = _calc_branch_values_from_trafo_df(net, ppc, trafo_df, sequence=0)

    # Y0y0d5,  YN0y0d5,  Y0yn0d5,  YN0yn0d5, Y0y0y0, Y0d5d5,
    # YN0d5d5,  Y0d5y0,  Y0y0d11  und  D0d0d0

    #
    # already implemented:
    #   YNyd
    #   YYnd
    #   YNynd
    #   YNdd
    # not relevant (for 1ph):
    #   Yyy
    #   Ydd
    #   Ydy
    #   Yyd
    #   Ddd

    BIG_NUMBER = 1e20 * ppc["baseMVA"]

    n_t3 = net.trafo3w.shape[0]
    for t3_ix in np.arange(n_t3):
        t3 = net.trafo3w.iloc[t3_ix, :]

        if t3.vector_group.lower() in set(map(lambda vg: "".join(vg), product("dy", repeat=3))):
            x[[t3_ix, t3_ix+n_t3, t3_ix+n_t3*2]] = BIG_NUMBER
            r[[t3_ix, t3_ix+n_t3, t3_ix+n_t3*2]] = BIG_NUMBER
        elif t3.vector_group.lower() == "ynyd":
            # Correction for YNyd
            # z3->y3
            ys = ppc["baseMVA"] / ((x[t3_ix+n_t3*2] * 1j + r[t3_ix+n_t3*2]) * ratio[t3_ix+n_t3*2] ** 2)
            aux_bus = bus_lookup[lv_bus[t3_ix]]
            ppc["bus"][aux_bus, BS] += ys.imag
            ppc["bus"][aux_bus, GS] += ys.real

            # Set y2/y3 to almost 0 to avoid isolated bus
            x[[t3_ix+n_t3, t3_ix+n_t3*2]] = BIG_NUMBER
            r[[t3_ix+n_t3, t3_ix+n_t3*2]] = BIG_NUMBER
        elif t3.vector_group.lower() == "yndy":
            # Correction for YNyd
            # z3->y3
            ys = ppc["baseMVA"] / ((x[t3_ix + n_t3] * 1j + r[t3_ix + n_t3]) * ratio[t3_ix + n_t3] ** 2)
            aux_bus = bus_lookup[lv_bus[t3_ix]]
            ppc["bus"][aux_bus, BS] += ys.imag
            ppc["bus"][aux_bus, GS] += ys.real

            # Set y2/y3 to almost 0 to avoid isolated bus
            x[[t3_ix+n_t3, t3_ix+n_t3*2]] = BIG_NUMBER
            r[[t3_ix+n_t3, t3_ix+n_t3*2]] = BIG_NUMBER
        elif t3.vector_group.lower() == "yynd":
            # z3->y3
            ys = ppc["baseMVA"] / ((x[t3_ix+n_t3*2] * 1j + r[t3_ix+n_t3*2]) * ratio[t3_ix+n_t3*2] ** 2)
            aux_bus = bus_lookup[lv_bus[t3_ix]]
            ppc["bus"][aux_bus, BS] += ys.imag
            ppc["bus"][aux_bus, GS] += ys.real

            # Set y1/y3 to almost 0 to avoid isolated bus
            x[[t3_ix, t3_ix+n_t3*2]] = BIG_NUMBER
            r[[t3_ix, t3_ix+n_t3*2]] = BIG_NUMBER
        elif t3.vector_group.lower() == "ydyn":
            # z3->y3
            ys = ppc["baseMVA"] / ((x[t3_ix+n_t3] * 1j + r[t3_ix+n_t3]) * ratio[t3_ix+n_t3] ** 2)
            aux_bus = bus_lookup[lv_bus[t3_ix]]
            ppc["bus"][aux_bus, BS] += ys.imag
            ppc["bus"][aux_bus, GS] += ys.real

            # Set y1/y3 to almost 0 to avoid isolated bus
            x[[t3_ix, t3_ix+n_t3]] = BIG_NUMBER
            r[[t3_ix, t3_ix+n_t3]] = BIG_NUMBER
        elif t3.vector_group.lower() == "ynynd":
            # z3->y3
            ys = ppc["baseMVA"] / ((x[t3_ix+n_t3*2] * 1j + r[t3_ix+n_t3*2]) * ratio[t3_ix+n_t3*2] ** 2)
            aux_bus = bus_lookup[lv_bus[t3_ix]]
            ppc["bus"][aux_bus, BS] += ys.imag
            ppc["bus"][aux_bus, GS] += ys.real

            # Set y3 to almost 0 to avoid isolated bus
            x[t3_ix+n_t3*2] = BIG_NUMBER
            r[t3_ix+n_t3*2] = BIG_NUMBER
        elif t3.vector_group.lower() == "yndyn":
            # z3->y3
            ys = ppc["baseMVA"] / ((x[t3_ix + n_t3] * 1j + r[t3_ix + n_t3]) * ratio[t3_ix + n_t3] ** 2)
            aux_bus = bus_lookup[lv_bus[t3_ix]]
            ppc["bus"][aux_bus, BS] += ys.imag
            ppc["bus"][aux_bus, GS] += ys.real

            # Set y3 to almost 0 to avoid isolated bus
            x[t3_ix + n_t3] = BIG_NUMBER
            r[t3_ix + n_t3] = BIG_NUMBER
        elif t3.vector_group.lower() == "yndd":
            # Correction for YNdd
            # z3->y3
            # we need a shunt impedance for both "delta" windings -> ys1, ys2
            ys1 = ppc["baseMVA"] / ((x[t3_ix + n_t3] * 1j + r[t3_ix + n_t3]) * ratio[t3_ix + n_t3] ** 2)
            ys2 = ppc["baseMVA"] / ((x[t3_ix + n_t3 * 2] * 1j + r[t3_ix + n_t3 * 2]) * ratio[t3_ix + n_t3 * 2] ** 2)
            aux_bus = bus_lookup[lv_bus[t3_ix]]
            ppc["bus"][aux_bus, BS] += ys1.imag + ys2.imag
            ppc["bus"][aux_bus, GS] += ys1.real + ys2.real

            # Set y2/y3 to almost 0 to avoid isolated bus
            x[[t3_ix + n_t3, t3_ix + n_t3 * 2]] = BIG_NUMBER
            r[[t3_ix + n_t3, t3_ix + n_t3 * 2]] = BIG_NUMBER
        elif t3.vector_group.lower() == "ynyy":
            # Correction for YNyy
            x[[t3_ix, t3_ix + n_t3, t3_ix + n_t3 * 2]] = BIG_NUMBER
            r[[t3_ix, t3_ix + n_t3, t3_ix + n_t3 * 2]] = BIG_NUMBER
        else:
            raise UserWarning(f"{t3.vector_group} not supported yet for trafo3w!")

    branch[f:t, BR_R] = r
    branch[f:t, BR_X] = x
    branch[f:t, BR_B] = 0
    branch[f:t, TAP] = ratio
    branch[f:t, SHIFT] = shift
    branch[f:t, BR_STATUS] = in_service
