Repository: e2nIEE/pandapower

Issue:
calc_sc branch calculation error
I have been investigating the calc_sc function with and without a fault impedance and met with a behavior I cannot easily explain. I know this is still in beta mode.

I used the tutorial on SC analysis from [here](https://github.com/e2nIEE/pandapower/blob/master/tutorials/shortcircuit.ipynb).

I compute the 3ph short circuit with fault impedance Z_F=0 and then again with Z_F=5+j10 Ohm. While the max fault currents at the nodes are decreasing as expected, the max fault currents on the branches are many times higher. 

As I understand, the currents are computed based on the post-fault voltages and the Ybus matrix. However, I expected when the short-circuit currents on each node are decreasing so would the max short-circuit current on the branches (roughly based on Kirchoff's law). Is this a bug or something I don't understand?

<img width="844" alt="pp_sc_issue" src="https://user-images.githubusercontent.com/10074820/99818234-a1671000-2b56-11eb-82ac-2668f9aaa92b.png">


Discussion:

Patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 6c0cd51c..c8675895 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,5 +1,7 @@
 Change Log
 =============
+- [ADDED] Factorization mode instead of inversion of Ybus in short-circuit calculation
+- [ADDED] Optimized the calculation of single/selected buses in 1ph/2ph/3ph short-circuit calculation
 
 [2.5.0]- 2021-01-08
 ----------------------
diff --git a/pandapower/auxiliary.py b/pandapower/auxiliary.py
index b537629f..b3948e43 100644
--- a/pandapower/auxiliary.py
+++ b/pandapower/auxiliary.py
@@ -577,7 +577,8 @@ def _add_opf_options(net, trafo_loading, ac, v_debug=False, **kwargs):
 
 
 def _add_sc_options(net, fault, case, lv_tol_percent, tk_s, topology, r_fault_ohm,
-                    x_fault_ohm, kappa, ip, ith, branch_results, kappa_method, return_all_currents):
+                    x_fault_ohm, kappa, ip, ith, branch_results, kappa_method, return_all_currents,
+                    inverse_y):
     """
     creates dictionary for pf, opf and short circuit calculations from input parameters.
     """
@@ -594,7 +595,8 @@ def _add_sc_options(net, fault, case, lv_tol_percent, tk_s, topology, r_fault_oh
         "ith": ith,
         "branch_results": branch_results,
         "kappa_method": kappa_method,
-        "return_all_currents": return_all_currents
+        "return_all_currents": return_all_currents,
+        "inverse_y": inverse_y
     }
     _add_options(net, options)
 
diff --git a/pandapower/results.py b/pandapower/results.py
index cf8882f7..8df714b2 100644
--- a/pandapower/results.py
+++ b/pandapower/results.py
@@ -189,14 +189,23 @@ def _ppci_other_to_ppc(result, ppc, mode):
 
 def _ppci_internal_to_ppc(result, ppc):
     for key, value in result["internal"].items():
+        # Only for sc calculation
         # if branch current matrices have been stored they need to include out of service elements
-        if key in ["branch_ikss_f", "branch_ikss_t", "branch_ip_f", "branch_ip_t", "branch_ith_f", "branch_ith_t"]:
-            n_buses = np.shape(ppc['bus'])[0]
+        if key in ["branch_ikss_f", "branch_ikss_t",
+                   "branch_ip_f", "branch_ip_t",
+                   "branch_ith_f", "branch_ith_t"]:
+
+            # n_buses = np.shape(ppc['bus'])[0]
             n_branches = np.shape(ppc['branch'])[0]
-            n_rows_result = np.shape(result['bus'])[0]
-            update_matrix = np.empty((n_branches, n_buses)) * np.nan
-            update_matrix[result["internal"]['branch_is'], :n_rows_result] = result["internal"][key]
-            ppc['internal'][key] = np.copy(update_matrix)
+            # n_rows_result = np.shape(result['bus'])[0]
+            # update_matrix = np.empty((n_branches, n_buses)) * np.nan
+            # update_matrix[result["internal"]['branch_is'], :n_rows_result] = result["internal"][key]
+
+            # # To select only required buses and pad one column of nan value for oos bus
+            update_matrix = np.empty((n_branches, value.shape[1]+1)) * 0.0
+            update_matrix[result["internal"]['branch_is'],
+                          :value.shape[1]] = result["internal"][key]
+            ppc['internal'][key] = update_matrix
         else:
             ppc["internal"][key] = value
 
diff --git a/pandapower/shortcircuit/calc_sc.py b/pandapower/shortcircuit/calc_sc.py
index b5f0c23e..5ae78cfd 100644
--- a/pandapower/shortcircuit/calc_sc.py
+++ b/pandapower/shortcircuit/calc_sc.py
@@ -10,14 +10,17 @@ except ImportError:
     import logging
 
 logger = logging.getLogger(__name__)
-# import time
+
+import numpy as np
+from scipy.sparse.linalg import factorized
 
 from pandapower.auxiliary import _clean_up, _add_ppc_options, _add_sc_options, _add_auxiliary_elements
 from pandapower.pd2ppc import _pd2ppc
 from pandapower.pd2ppc_zero import _pd2ppc_zero
 from pandapower.results import _copy_results_ppci_to_ppc
-from pandapower.shortcircuit.currents import _calc_ikss, _calc_ikss_1ph, _calc_ip, _calc_ith, _calc_branch_currents, \
-    _calc_single_bus_sc
+from pandapower.shortcircuit.currents import _calc_ikss, \
+    _calc_ikss_1ph, _calc_ip, _calc_ith, _calc_branch_currents, \
+    _calc_single_bus_sc, _calc_single_bus_sc_no_y_inv
 from pandapower.shortcircuit.impedance import _calc_zbus, _calc_ybus, _calc_rx
 from pandapower.shortcircuit.kappa import _add_kappa_to_ppc
 from pandapower.shortcircuit.results import _extract_results, _extract_single_results
@@ -26,7 +29,8 @@ from pandapower.results import init_results
 
 def calc_sc(net, fault="3ph", case='max', lv_tol_percent=10, topology="auto", ip=False,
             ith=False, tk_s=1., kappa_method="C", r_fault_ohm=0., x_fault_ohm=0.,
-            branch_results=False, check_connectivity=True, return_all_currents=False):
+            branch_results=False, check_connectivity=True, return_all_currents=False,
+            bus=None, inverse_y=True):
     """
     Calculates minimal or maximal symmetrical short-circuit currents.
     The calculation is based on the method of the equivalent voltage source
@@ -85,6 +89,10 @@ def calc_sc(net, fault="3ph", case='max', lv_tol_percent=10, topology="auto", ip
         **return_all_currents** (bool, False) applies only if branch_results=True, if True short-circuit currents for
         each (branch, bus) tuple is returned otherwise only the max/min is returned
 
+        **bus** (int, list, np.array, None) defines if short-circuit calculations should only be calculated for defined bus
+
+        **inverse_y** (bool, True) defines if complete inverse should be used instead of LU factorization, factorization version is in experiment which should be faster and memory efficienter
+
 
     OUTPUT:
 
@@ -112,6 +120,12 @@ def calc_sc(net, fault="3ph", case='max', lv_tol_percent=10, topology="auto", ip
         logger.warning("Branch results are in beta mode and might not always be reliable, "
                        "especially for transformers")
 
+    # Convert bus to numpy array for better performance
+    if isinstance(bus, int):
+        bus = np.array([bus])
+    elif isinstance(bus, list):
+        bus = np.array(bus)
+
     kappa = ith or ip
     net["_options"] = {}
     _add_ppc_options(net, calculate_voltage_angles=False, trafo_model="pi",
@@ -121,17 +135,19 @@ def calc_sc(net, fault="3ph", case='max', lv_tol_percent=10, topology="auto", ip
     _add_sc_options(net, fault=fault, case=case, lv_tol_percent=lv_tol_percent, tk_s=tk_s,
                     topology=topology, r_fault_ohm=r_fault_ohm, kappa_method=kappa_method,
                     x_fault_ohm=x_fault_ohm, kappa=kappa, ip=ip, ith=ith,
-                    branch_results=branch_results, return_all_currents=return_all_currents)
+                    branch_results=branch_results, return_all_currents=return_all_currents,
+                    inverse_y=inverse_y)
     init_results(net, "sc")
-    if fault == "3ph":
-        _calc_sc(net)
-    if fault == "2ph":
-        _calc_sc(net)
-    if fault == "1ph":
-        _calc_sc_1ph(net)
+    if fault in ("2ph", "3ph"):
+        _calc_sc(net, bus)
+    elif fault == "1ph":
+        _calc_sc_1ph(net, bus)
+    else:
+        raise ValueError("Invalid fault %s" % fault)
 
 
-def calc_single_sc(net, bus, fault="3ph", case='max', lv_tol_percent=10, check_connectivity=True):
+def calc_single_sc(net, bus, fault="3ph", case='max', lv_tol_percent=10,
+                   check_connectivity=True, inverse_y=True):
     """
     Calculates minimal or maximal symmetrical short-circuit currents.
     The calculation is based on the method of the equivalent voltage source
@@ -190,9 +206,10 @@ def calc_single_sc(net, bus, fault="3ph", case='max', lv_tol_percent=10, check_c
     _add_sc_options(net, fault=fault, case=case, lv_tol_percent=lv_tol_percent, tk_s=1.,
                     topology="auto", r_fault_ohm=0., kappa_method="C",
                     x_fault_ohm=0., kappa=False, ip=False, ith=False,
-                    branch_results=True, return_all_currents=False)
+                    branch_results=True, return_all_currents=False,
+                    inverse_y=inverse_y)
     init_results(net, "sc")
-    if fault == "3ph" or fault == "2ph":
+    if fault in ("2ph", "3ph"):
         _calc_sc_single(net, bus)
     elif fault == "1ph":
         raise NotImplementedError("1ph short-circuits are not yet implemented")
@@ -204,48 +221,64 @@ def _calc_sc_single(net, bus):
     _add_auxiliary_elements(net)
     ppc, ppci = _pd2ppc(net)
     _calc_ybus(ppci)
-    try:
-        _calc_zbus(ppci)
-    except Exception as e:
-        _clean_up(net, res=False)
-        raise (e)
-    _calc_rx(net, ppci)
-    _calc_ikss(net, ppci)
-    _calc_single_bus_sc(net, ppci, bus)
+
+    if net["_options"]["inverse_y"]:
+        _calc_zbus(net, ppci)
+        _calc_rx(net, ppci, bus=None)
+        _calc_ikss(net, ppci, bus=None)
+        _calc_single_bus_sc(net, ppci, bus)
+    else:
+        # Factorization Ybus once
+        ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])
+
+        _calc_rx(net, ppci, bus)
+        _calc_ikss(net, ppci, bus)
+        _calc_single_bus_sc_no_y_inv(net, ppci, bus)
+
+        # Delete factorization object
+        ppci["internal"].pop("ybus_fact")
+
     ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
     _extract_single_results(net, ppc)
     _clean_up(net)
 
 
-def _calc_sc(net):
+def _calc_sc(net, bus):
     _add_auxiliary_elements(net)
     ppc, ppci = _pd2ppc(net)
     _calc_ybus(ppci)
-    try:
-        _calc_zbus(ppci)
-    except Exception as e:
-        _clean_up(net, res=False)
-        raise (e)
-    _calc_rx(net, ppci)
-    _add_kappa_to_ppc(net, ppci)
-    _calc_ikss(net, ppci)
+
+    if net["_options"]["inverse_y"]:
+        _calc_zbus(net, ppci)
+    else:
+        # Factorization Ybus once
+        ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])
+
+    _calc_rx(net, ppci, bus)
+
+    # kappa required inverse of Zbus, which is optimized
+    if net["_options"]["kappa"]:
+        _add_kappa_to_ppc(net, ppci)
+    _calc_ikss(net, ppci, bus)
+
     if net["_options"]["ip"]:
         _calc_ip(net, ppci)
     if net["_options"]["ith"]:
         _calc_ith(net, ppci)
+
     if net._options["branch_results"]:
-        _calc_branch_currents(net, ppci)
+        _calc_branch_currents(net, ppci, bus)
 
     ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
-
-    if net["_options"]["return_all_currents"]:
-        _extract_results(net, ppc, ppc_0=None)
-    else:
-        _extract_results(net, ppc, ppc_0=None)
+    _extract_results(net, ppc, ppc_0=None, bus=bus)
     _clean_up(net)
 
+    if "ybus_fact" in ppci["internal"]:
+        # Delete factorization object
+        ppci["internal"].pop("ybus_fact")
 
-def _calc_sc_1ph(net):
+
+def _calc_sc_1ph(net, bus):
     """
     calculation method for single phase to ground short-circuit currents
     """
@@ -253,26 +286,28 @@ def _calc_sc_1ph(net):
     # pos. seq bus impedance
     ppc, ppci = _pd2ppc(net)
     _calc_ybus(ppci)
-    try:
-        _calc_zbus(ppci)
-    except Exception as e:
-        _clean_up(net, res=False)
-        raise (e)
-    _calc_rx(net, ppci)
-    _add_kappa_to_ppc(net, ppci)
+
     # zero seq bus impedance
     ppc_0, ppci_0 = _pd2ppc_zero(net)
     _calc_ybus(ppci_0)
-    try:
-        _calc_zbus(ppci_0)
-    except Exception as e:
-        _clean_up(net, res=False)
-        raise (e)
-    _calc_rx(net, ppci_0)
-    _calc_ikss_1ph(net, ppci, ppci_0)
+
+    if net["_options"]["inverse_y"]:
+        _calc_zbus(net, ppci)
+        _calc_zbus(net, ppci_0)
+    else:
+        # Factorization Ybus once
+        ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])
+        ppci_0["internal"]["ybus_fact"] = factorized(ppci_0["internal"]["Ybus"])
+
+    _calc_rx(net, ppci, bus=bus)
+    _add_kappa_to_ppc(net, ppci)
+
+    _calc_rx(net, ppci_0, bus=bus)
+    _calc_ikss_1ph(net, ppci, ppci_0, bus=bus)
+
     if net._options["branch_results"]:
-        _calc_branch_currents(net, ppci)
+        _calc_branch_currents(net, ppci, bus=bus)
     ppc_0 = _copy_results_ppci_to_ppc(ppci_0, ppc_0, "sc")
     ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
-    _extract_results(net, ppc, ppc_0)
+    _extract_results(net, ppc, ppc_0, bus=bus)
     _clean_up(net)
diff --git a/pandapower/shortcircuit/currents.py b/pandapower/shortcircuit/currents.py
index f4bc9540..12a5bd60 100644
--- a/pandapower/shortcircuit/currents.py
+++ b/pandapower/shortcircuit/currents.py
@@ -12,33 +12,51 @@ from pandapower.pypower.idx_bus import BASE_KV
 from pandapower.pypower.idx_gen import GEN_BUS, MBASE
 from pandapower.shortcircuit.idx_brch import IKSS_F, IKSS_T, IP_F, IP_T, ITH_F, ITH_T
 from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX, KAPPA, R_EQUIV, IKSS1, IP, ITH, X_EQUIV, IKSS2, IKCV, M
+from pandapower.shortcircuit.impedance import _calc_zbus_diag
 
+from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
+from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci
+
+def _calc_ikss(net, ppc, bus=None):
+    # Vectorized for multiple bus
+    if bus is None:
+        # Slice(None) is equal to : select
+        bus_idx = slice(None)
+    else:
+        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)
 
-def _calc_ikss(net, ppc):
     fault = net._options["fault"]
     case = net._options["case"]
-    c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
+    c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
     ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
-    z_equiv = abs(ppc["bus"][:, R_EQUIV] + ppc["bus"][:, X_EQUIV] * 1j)
+
+    z_equiv = abs(ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j)
     if fault == "3ph":
-        ppc["bus"][:, IKSS1] = c / z_equiv / ppc["bus"][:, BASE_KV] / np.sqrt(3) * ppc["baseMVA"]
+        ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / np.sqrt(3) * ppc["baseMVA"]
     elif fault == "2ph":
-        ppc["bus"][:, IKSS1] = c / z_equiv / ppc["bus"][:, BASE_KV] / 2 * ppc["baseMVA"]
+        ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / 2 * ppc["baseMVA"]
     _current_source_current(net, ppc)
 
 
-def _calc_ikss_1ph(net, ppc, ppc_0):
-    fault = net._options["fault"]
+def _calc_ikss_1ph(net, ppc, ppc_0, bus=None):
+    # Vectorized for multiple bus
+    if bus is None:
+        # Slice(None) is equal to : select
+        bus_idx = slice(None)
+    else:
+        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)
+
     case = net._options["case"]
-    c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
+    c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
     ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
     ppc_0["internal"]["baseI"] = ppc_0["bus"][:, BASE_KV] * np.sqrt(3) / ppc_0["baseMVA"]
-    z_equiv = abs((ppc["bus"][:, R_EQUIV] + ppc["bus"][:, X_EQUIV] * 1j) * 2 +
-                  (ppc_0["bus"][:, R_EQUIV] + ppc_0["bus"][:, X_EQUIV] * 1j))
-    if fault == "1ph":
-        ppc_0["bus"][:, IKSS1] = c / z_equiv / ppc_0["bus"][:, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
-        ppc["bus"][:, IKSS1] = c / z_equiv / ppc_0["bus"][:, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
-        _current_source_current(net, ppc)
+
+    z_equiv = abs((ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j) * 2 +
+                  (ppc_0["bus"][bus_idx, R_EQUIV] + ppc_0["bus"][bus_idx, X_EQUIV] * 1j))
+
+    ppc_0["bus"][bus_idx, IKSS1] = c / z_equiv / ppc_0["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
+    ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppc["baseMVA"]
+    _current_source_current(net, ppc)
 
 
 def _current_source_current(net, ppc):
@@ -56,13 +74,18 @@ def _current_source_current(net, ppc):
     baseI = ppc["internal"]["baseI"]
     sgen_buses = sgen.bus.values
     sgen_buses_ppc = bus_lookup[sgen_buses]
-    Zbus = ppc["internal"]["Zbus"]
     if not "k" in sgen:
         raise ValueError("Nominal to short-circuit current has to specified in net.sgen.k")
     i_sgen_pu = sgen.sn_mva.values / net.sn_mva * sgen.k.values
     buses, ikcv_pu, _ = _sum_by_group(sgen_buses_ppc, i_sgen_pu, i_sgen_pu)
     ppc["bus"][buses, IKCV] = ikcv_pu
-    ppc["bus"][:, IKSS2] = abs(1 / np.diag(Zbus) * np.dot(Zbus, ppc["bus"][:, IKCV] * -1j) / baseI)
+    if net["_options"]["inverse_y"]:
+        Zbus = ppc["internal"]["Zbus"]
+        ppc["bus"][:, IKSS2] = abs(1 / np.diag(Zbus) * np.dot(Zbus, ppc["bus"][:, IKCV] * -1j) / baseI)
+    else:
+        ybus_fact = ppc["internal"]["ybus_fact"]
+        diagZ = _calc_zbus_diag(net, ppc)
+        ppc["bus"][:, IKSS2] = abs(ybus_fact(ppc["bus"][:, IKCV] * -1j) / diagZ / baseI)
     ppc["bus"][buses, IKCV] /= baseI[buses]
 
 
@@ -83,84 +106,6 @@ def _calc_ith(net, ppc):
     ppc["bus"][:, ITH] = ith
 
 
-def _calc_branch_currents(net, ppc):
-    case = net._options["case"]
-    Zbus = ppc["internal"]["Zbus"]
-    Yf = ppc["internal"]["Yf"]
-    Yt = ppc["internal"]["Yt"]
-    baseI = ppc["internal"]["baseI"]
-    n = ppc["bus"].shape[0]
-    fb = np.real(ppc["branch"][:, 0]).astype(int)
-    tb = np.real(ppc["branch"][:, 1]).astype(int)
-    minmax = np.nanmin if case == "min" else np.nanmax
-
-    # calculate voltage source branch current
-    V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
-    ikss1_all_f = np.conj(Yf.dot(V_ikss))
-    ikss1_all_t = np.conj(Yt.dot(V_ikss))
-    ikss1_all_f[abs(ikss1_all_f) < 1e-10] = 0.
-    ikss1_all_t[abs(ikss1_all_t) < 1e-10] = 0.
-
-    # add current source branch current if there is one
-    current_sources = any(ppc["bus"][:, IKCV]) > 0
-    if current_sources:
-        current = np.tile(-ppc["bus"][:, IKCV], (n, 1))
-        np.fill_diagonal(current, current.diagonal() + ppc["bus"][:, IKSS2])
-        V = np.dot((current * baseI), Zbus).T
-        fb = np.real(ppc["branch"][:, 0]).astype(int)
-        tb = np.real(ppc["branch"][:, 1]).astype(int)
-        ikss2_all_f = np.conj(Yf.dot(V))
-        ikss2_all_t = np.conj(Yt.dot(V))
-        ikss_all_f = abs(ikss1_all_f + ikss2_all_f)
-        ikss_all_t = abs(ikss1_all_t + ikss2_all_t)
-    else:
-        ikss_all_f = abs(ikss1_all_f)
-        ikss_all_t = abs(ikss1_all_t)
-
-    if net._options["return_all_currents"]:
-        ppc["internal"]["branch_ikss_f"] = ikss_all_f / baseI[fb, None]
-        ppc["internal"]["branch_ikss_t"] = ikss_all_t / baseI[tb, None]
-    else:
-        ikss_all_f[ikss_all_f < 1e-10] = np.nan
-        ikss_all_t[ikss_all_t < 1e-10] = np.nan
-        ppc["branch"][:, IKSS_F] = minmax(ikss_all_f, axis=1) / baseI[fb]
-        ppc["branch"][:, IKSS_T] = minmax(ikss_all_t, axis=1) / baseI[tb]
-
-    if net._options["ip"]:
-        kappa = ppc["bus"][:, KAPPA]
-        if current_sources:
-            ip_all_f = np.sqrt(2) * (ikss1_all_f * kappa + ikss2_all_f)
-            ip_all_t = np.sqrt(2) * (ikss1_all_t * kappa + ikss2_all_t)
-        else:
-            ip_all_f = np.sqrt(2) * ikss1_all_f * kappa
-            ip_all_t = np.sqrt(2) * ikss1_all_t * kappa
-
-        if net._options["return_all_currents"]:
-            ppc["internal"]["branch_ip_f"] = abs(ip_all_f) / baseI[fb, None]
-            ppc["internal"]["branch_ip_t"] = abs(ip_all_t) / baseI[tb, None]
-        else:
-            ip_all_f[abs(ip_all_f) < 1e-10] = np.nan
-            ip_all_t[abs(ip_all_t) < 1e-10] = np.nan
-            ppc["branch"][:, IP_F] = minmax(abs(ip_all_f), axis=1) / baseI[fb]
-            ppc["branch"][:, IP_T] = minmax(abs(ip_all_t), axis=1) / baseI[tb]
-
-    if net._options["ith"]:
-        n = 1
-        m = ppc["bus"][:, M]
-        ith_all_f = ikss_all_f * np.sqrt(m + n)
-        ith_all_t = ikss_all_t * np.sqrt(m + n)
-
-        if net._options["return_all_currents"]:
-            ppc["internal"]["branch_ith_f"] = ith_all_f / baseI[fb, None]
-            ppc["internal"]["branch_ith_t"] = ith_all_t / baseI[tb, None]
-        else:
-            ppc["branch"][:, ITH_F] = minmax(ith_all_f, axis=1) / baseI[fb]
-            ppc["branch"][:, ITH_T] = minmax(ith_all_t, axis=1) / baseI[fb]
-
-
-
-
-
 def _calc_ib_generator(net, ppci):
     Zbus = ppci["internal"]["Zbus"]
     baseI = ppci["internal"]["baseI"]
@@ -230,7 +175,7 @@ def _calc_ib_generator(net, ppci):
 
 
 def _calc_single_bus_sc(net, ppc, bus):
-    case = net._options["case"]
+    # case = net._options["case"]
     bus_idx = net._pd2ppc_lookups["bus"][bus]
     n = ppc["bus"].shape[0]
     Zbus = ppc["internal"]["Zbus"]
@@ -239,7 +184,7 @@ def _calc_single_bus_sc(net, ppc, bus):
     baseI = ppc["internal"]["baseI"]
     #    fb = np.real(ppc["branch"][:, 0]).astype(int)
     #    tb = np.real(ppc["branch"][:, 1]).astype(int)
-    c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
+    # c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
 
     # calculate voltage source branch current
     V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
@@ -257,15 +202,167 @@ def _calc_single_bus_sc(net, ppc, bus):
     #    ppc["branch"][:, IKSS_T] = abs(ikss_all_t[:, bus_idx] / baseI[tb])
     calc_branch_results(net, ppc, V)
 
+def _calc_single_bus_sc_no_y_inv(net, ppc, bus):
+    # Vectorized for multiple bus
+    if bus is None:
+        # Slice(None) is equal to : select
+        bus_idx = slice(None)
+    else:
+        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)
 
-from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
-from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci
+    ybus = ppc["internal"]["Ybus"]
+    ybus_fact = ppc["internal"]["ybus_fact"]
+    # case = net._options["case"]
+    baseI = ppc["internal"]["baseI"]
+    # vqj = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX] #this is the source voltage in per unit (VQj)
+
+    # Solve Ikss from voltage source
+    n_bus = ybus.shape[0]
+
+    # ybus_sub_mask = (np.arange(ybus.shape[0]) != bus_idx)
+    # V_ikss = np.zeros(n_bus, dtype=np.complex)
+    # V_ikss[bus_idx] = vqj[bus_idx]
+
+    # # Solve Ax = b
+    # b = np.zeros(n_bus-1, dtype=np.complex) -\
+    #     (ybus[:, ~ybus_sub_mask].toarray())[ybus_sub_mask].ravel() * V_ikss[bus_idx]
+    # ybus_sub = ybus[ybus_sub_mask, :][:, ybus_sub_mask]
+    # x = spsolve(ybus_sub, b)
+
+    # V_ikss[ybus_sub_mask] = x
+    # I_ikss = np.zeros(n_bus, dtype=np.complex)
+    # I_ikss[bus_idx] = np.dot(ybus[bus_idx, :].toarray(), V_ikss)
+    # V = V_ikss
+
+    # Version 2
+    I_ikss = np.zeros(n_bus, dtype=np.complex)
+    I_ikss[bus_idx] = ppc["bus"][bus_idx, IKSS1]
+    V_ikss = ybus_fact(I_ikss * baseI)
+    V = V_ikss
+
+    #TODO include current sources
+    current_sources = any(ppc["bus"][:, IKCV]) > 0
+    if current_sources:
+        current = -ppc["bus"][:, IKCV]
+        current[bus_idx] += ppc["bus"][bus_idx, IKSS2]
+        V_source = ybus_fact(current)
+        V += V_source
+
+    calc_branch_results(net, ppc, V)
 
 
 def calc_branch_results(net, ppci, V):
     Ybus = ppci["internal"]["Ybus"]
     Yf = ppci["internal"]["Yf"]
     Yt = ppci["internal"]["Yt"]
-    baseMVA, bus, gen, branch, ref, _, pq, _, _, V0, ref_gens = _get_pf_variables_from_ppci(ppci)
+    baseMVA, bus, gen, branch, ref, _, _, _, _, _, ref_gens = _get_pf_variables_from_ppci(ppci)
     bus, gen, branch = pfsoln_pypower(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens)
     ppci["bus"], ppci["gen"], ppci["branch"] = bus, gen, branch
+
+
+def _calc_branch_currents(net, ppc, bus):
+    # Vectorized for multiple bus
+    if bus is None:
+        # Slice(None) is equal to select all
+        bus = net.bus.index
+
+    bus_idx = net._pd2ppc_lookups["bus"][bus]
+    # Select only in service bus for sc calculation
+    bus_idx = bus_idx[bus_idx < ppc['bus'].shape[0]]
+    n_sc_bus = np.shape(bus_idx)[0]
+
+    case = net._options["case"]
+
+    Yf = ppc["internal"]["Yf"]
+    Yt = ppc["internal"]["Yt"]
+    baseI = ppc["internal"]["baseI"]
+    n_bus = ppc["bus"].shape[0]
+    fb = np.real(ppc["branch"][:, 0]).astype(int)
+    tb = np.real(ppc["branch"][:, 1]).astype(int)
+    minmax = np.nanmin if case == "min" else np.nanmax
+
+    # calculate voltage source branch current
+    if net["_options"]["inverse_y"]:
+        Zbus = ppc["internal"]["Zbus"]
+        V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
+        V_ikss = V_ikss[:, bus_idx]
+    else:
+        ybus_fact = ppc["internal"]["ybus_fact"]
+        V_ikss = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
+        for ix, b in enumerate(bus_idx):
+            ikss = np.zeros(n_bus, dtype=np.complex)
+            ikss[b] = ppc["bus"][b, IKSS1] * baseI[b]
+            V_ikss[:, ix] = ybus_fact(ikss)
+
+    ikss1_all_f = np.conj(Yf.dot(V_ikss))
+    ikss1_all_t = np.conj(Yt.dot(V_ikss))
+    ikss1_all_f[abs(ikss1_all_f) < 1e-10] = 0.
+    ikss1_all_t[abs(ikss1_all_t) < 1e-10] = 0.
+
+    # add current source branch current if there is one
+    current_sources = any(ppc["bus"][:, IKCV]) > 0
+    if current_sources:
+        current = np.tile(-ppc["bus"][:, IKCV], (n_sc_bus, 1))
+        for ix, b in enumerate(bus_idx):
+            current[ix, b] += ppc["bus"][b, IKSS2]
+
+        # calculate voltage source branch current
+        if net["_options"]["inverse_y"]:
+            Zbus = ppc["internal"]["Zbus"]
+            V = np.dot((current * baseI), Zbus).T
+        else:
+            ybus_fact = ppc["internal"]["ybus_fact"]
+            V = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
+            for ix, b in enumerate(bus_idx):
+                V[:, ix] = ybus_fact(current[ix, :] * baseI[b])
+
+        fb = np.real(ppc["branch"][:, 0]).astype(int)
+        tb = np.real(ppc["branch"][:, 1]).astype(int)
+        ikss2_all_f = np.conj(Yf.dot(V))
+        ikss2_all_t = np.conj(Yt.dot(V))
+
+        ikss_all_f = abs(ikss1_all_f + ikss2_all_f)
+        ikss_all_t = abs(ikss1_all_t + ikss2_all_t)
+    else:
+        ikss_all_f = abs(ikss1_all_f)
+        ikss_all_t = abs(ikss1_all_t)
+
+    if net._options["return_all_currents"]:
+        ppc["internal"]["branch_ikss_f"] = ikss_all_f / baseI[fb, None]
+        ppc["internal"]["branch_ikss_t"] = ikss_all_t / baseI[tb, None]
+    else:
+        ikss_all_f[abs(ikss_all_f) < 1e-10] = np.nan
+        ikss_all_t[abs(ikss_all_t) < 1e-10] = np.nan
+        ppc["branch"][:, IKSS_F] = np.nan_to_num(minmax(ikss_all_f, axis=1) / baseI[fb])
+        ppc["branch"][:, IKSS_T] = np.nan_to_num(minmax(ikss_all_t, axis=1) / baseI[tb])
+
+    if net._options["ip"]:
+        kappa = ppc["bus"][:, KAPPA]
+        if current_sources:
+            ip_all_f = np.sqrt(2) * (ikss1_all_f * kappa[bus_idx] + ikss2_all_f)
+            ip_all_t = np.sqrt(2) * (ikss1_all_t * kappa[bus_idx] + ikss2_all_t)
+        else:
+            ip_all_f = np.sqrt(2) * ikss1_all_f * kappa[bus_idx]
+            ip_all_t = np.sqrt(2) * ikss1_all_t * kappa[bus_idx]
+
+        if net._options["return_all_currents"]:
+            ppc["internal"]["branch_ip_f"] = abs(ip_all_f) / baseI[fb, None]
+            ppc["internal"]["branch_ip_t"] = abs(ip_all_t) / baseI[tb, None]
+        else:
+            ip_all_f[abs(ip_all_f) < 1e-10] = np.nan
+            ip_all_t[abs(ip_all_t) < 1e-10] = np.nan
+            ppc["branch"][:, IP_F] = np.nan_to_num(minmax(abs(ip_all_f), axis=1) / baseI[fb])
+            ppc["branch"][:, IP_T] = np.nan_to_num(minmax(abs(ip_all_t), axis=1) / baseI[tb])
+
+    if net._options["ith"]:
+        n = 1
+        m = ppc["bus"][bus_idx, M]
+        ith_all_f = ikss_all_f * np.sqrt(m + n)
+        ith_all_t = ikss_all_t * np.sqrt(m + n)
+
+        if net._options["return_all_currents"]:
+            ppc["internal"]["branch_ith_f"] = ith_all_f / baseI[fb, None]
+            ppc["internal"]["branch_ith_t"] = ith_all_t / baseI[tb, None]
+        else:
+            ppc["branch"][:, ITH_F] = np.nan_to_num(minmax(ith_all_f, axis=1) / baseI[fb])
+            ppc["branch"][:, ITH_T] = np.nan_to_num(minmax(ith_all_t, axis=1) / baseI[fb])
diff --git a/pandapower/shortcircuit/impedance.py b/pandapower/shortcircuit/impedance.py
index 8b988540..f07aef33 100644
--- a/pandapower/shortcircuit/impedance.py
+++ b/pandapower/shortcircuit/impedance.py
@@ -13,23 +13,37 @@ from scipy.linalg import inv
 
 from pandapower.shortcircuit.idx_bus import R_EQUIV, X_EQUIV
 from pandapower.pypower.idx_bus import BASE_KV
+from pandapower.auxiliary import _clean_up
 try:
     from pandapower.pf.makeYbus_numba import makeYbus
 except ImportError:
     from pandapower.pypower.makeYbus import makeYbus
 
 
-def _calc_rx(net, ppc):
-    Zbus = ppc["internal"]["Zbus"]
+def _calc_rx(net, ppc, bus):
+    # Vectorized for multiple bus
+    if bus is None:
+        # Slice(None) is equal to select all
+        bus_idx = slice(None)
+    else:
+        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)
+
     r_fault = net["_options"]["r_fault_ohm"]
     x_fault = net["_options"]["x_fault_ohm"]
     if r_fault > 0 or x_fault > 0:
-        base_r = np.square(ppc["bus"][:, BASE_KV]) / ppc["baseMVA"]
+        base_r = np.square(ppc["bus"][bus_idx, BASE_KV]) / ppc["baseMVA"]
         fault_impedance = (r_fault + x_fault * 1j) / base_r
-        np.fill_diagonal(Zbus, Zbus.diagonal() + fault_impedance)
-    z_equiv = np.diag(Zbus)
-    ppc["bus"][:, R_EQUIV] = z_equiv.real
-    ppc["bus"][:, X_EQUIV] = z_equiv.imag
+    else:
+        fault_impedance = 0 + 0j
+
+    if net["_options"]["inverse_y"]:
+        Zbus = ppc["internal"]["Zbus"]
+        z_equiv = np.diag(Zbus)[bus_idx] + fault_impedance
+    else:
+        z_equiv = _calc_zbus_diag(net, ppc, bus) + fault_impedance
+    ppc["bus"][bus_idx, R_EQUIV] = z_equiv.real
+    ppc["bus"][bus_idx, X_EQUIV] = z_equiv.imag
+
 
 def _calc_ybus(ppc):
     Ybus, Yf, Yt = makeYbus(ppc["baseMVA"], ppc["bus"],  ppc["branch"])
@@ -39,12 +53,60 @@ def _calc_ybus(ppc):
     ppc["internal"]["Yt"] = Yt
     ppc["internal"]["Ybus"] = Ybus
 
-def _calc_zbus(ppc):
-    Ybus = ppc["internal"]["Ybus"]
-    sparsity = Ybus.nnz / Ybus.shape[0]**2
-    if sparsity < 0.002:
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore")
-            ppc["internal"]["Zbus"] = inv_sparse(Ybus).toarray()
+
+def _calc_zbus(net, ppc):
+    try:
+        Ybus = ppc["internal"]["Ybus"]
+        sparsity = Ybus.nnz / Ybus.shape[0]**2
+        if sparsity < 0.002:
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                ppc["internal"]["Zbus"] = inv_sparse(Ybus).toarray()
+        else:
+            ppc["internal"]["Zbus"] = inv(Ybus.toarray())
+    except Exception as e:
+        _clean_up(net, res=False)
+        raise (e)
+
+
+def _calc_zbus_diag(net, ppc, bus=None):
+    ybus_fact = ppc["internal"]["ybus_fact"]
+    n_bus = ppc["bus"].shape[0]
+
+    if bus is None:
+        diagZ = np.zeros(n_bus, dtype=np.complex)
+        for i in range(ppc["bus"].shape[0]):
+            b = np.zeros(n_bus, dtype=np.complex)
+            b[i] = 1 + 0j
+            diagZ[i] = ybus_fact(b)[i]
+        ppc["internal"]["diagZ"] = diagZ
+        return diagZ
     else:
-        ppc["internal"]["Zbus"] = inv(Ybus.toarray())
+        if isinstance(bus, int):
+            bus = np.array([bus])
+        diagZ = np.zeros(np.shape(bus)[0], dtype=np.complex)
+        for ix, b in enumerate(bus):
+            bus_idx = net._pd2ppc_lookups["bus"][b] #bus where the short-circuit is calculated (j)
+            b = np.zeros(n_bus, dtype=np.complex)
+            b[bus_idx] = 1 + 0j
+            diagZ[ix] = ybus_fact(b)[bus_idx]
+        return diagZ
+
+    # if bus is None:
+    #     bus = net.bus.index
+
+    # diagZ = np.zeros(np.shape(bus)[0], dtype=np.complex)
+    # ix = 0
+
+    # # Use windows size 32 to calculate Zbus
+    # while ix < np.shape(bus)[0]:
+    #     ix_end = min(ix+32, np.shape(bus)[0])
+    #     bus_idx = net._pd2ppc_lookups["bus"][bus[ix: ix_end]]
+    #     b = np.zeros((n_bus, (ix_end-ix)), dtype=np.complex)
+    #     for this_ix, this_bus_ix in enumerate(bus_idx):
+    #         b[this_bus_ix, this_ix] = 1 + 0j
+    #     res = ybus_fact(b)
+    #     for this_ix, this_bus_ix in enumerate(bus_idx):
+    #         diagZ[ix] = res[this_bus_ix, this_ix]
+    #     ix += 32
+    # return diagZ
\ No newline at end of file
diff --git a/pandapower/shortcircuit/kappa.py b/pandapower/shortcircuit/kappa.py
index ae03b07a..d03453b4 100644
--- a/pandapower/shortcircuit/kappa.py
+++ b/pandapower/shortcircuit/kappa.py
@@ -7,6 +7,7 @@
 import copy
 import networkx as nx
 import numpy as np
+from scipy.sparse.linalg import factorized
 
 from pandapower.pypower.idx_brch import F_BUS, T_BUS, BR_R, BR_X
 from pandapower.pypower.idx_bus import BUS_I, GS, BS, BASE_KV
@@ -29,6 +30,9 @@ def _add_kappa_to_ppc(net, ppc):
         raise ValueError("Unknown kappa method %s - specify B or C"%kappa_method)
     ppc["bus"][:, KAPPA] = kappa
 
+def _kappa(rx):
+    return 1.02 + .98 * np.exp(-3 * rx)
+
 def _kappa_method_c(net, ppc):
     if net.f_hz == 50:
         fc = 20
@@ -48,14 +52,17 @@ def _kappa_method_c(net, ppc):
     ppc_c["bus"][conductance, GS] = y_shunt.real[0]
     ppc_c["bus"][conductance, BS] = y_shunt.imag[0]
     _calc_ybus(ppc_c)
-    _calc_zbus(ppc_c)
-    _calc_rx(net, ppc_c)
+
+    if net["_options"]["inverse_y"]:
+        _calc_zbus(net, ppc_c)
+    else:
+        # Factorization Ybus once
+        ppc_c["internal"]["ybus_fact"] = factorized(ppc_c["internal"]["Ybus"])
+
+    _calc_rx(net, ppc_c, bus=None)
     rx_equiv_c = ppc_c["bus"][:, R_EQUIV] / ppc_c["bus"][:, X_EQUIV] * fc / net.f_hz
     return _kappa(rx_equiv_c)
 
-def _kappa(rx):
-    return 1.02 + .98 * np.exp(-3 * rx)
-
 def _kappa_method_b(net, ppc):
     topology = net._options["topology"]
     kappa_max = np.full(ppc["bus"].shape[0], 2.)
diff --git a/pandapower/shortcircuit/results.py b/pandapower/shortcircuit/results.py
index bc10d64f..13bf1441 100644
--- a/pandapower/shortcircuit/results.py
+++ b/pandapower/shortcircuit/results.py
@@ -9,25 +9,32 @@ import pandas as pd
 
 from pandapower.shortcircuit.idx_brch import IKSS_F, IKSS_T, IP_F, IP_T, ITH_F, ITH_T
 from pandapower.shortcircuit.idx_bus import IKSS1, IP, ITH, IKSS2
-from pandapower.pypower.idx_bus import VM, VA
+from pandapower.pypower.idx_bus import VM, VA, BUS_TYPE
 from pandapower.results_bus import _get_bus_idx, _set_buses_out_of_service
 from pandapower.results import _get_aranged_lookup, _get_branch_results
 from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX
 
+def _get_bus_ppc_idx_for_br_all_results(net, ppc, bus):
+    bus_lookup = net._pd2ppc_lookups["bus"]
+    if bus is None:
+        bus = net.bus.index
+
+    ppc_index = np.arange(np.shape(bus)[0])
+    ppc_index[ppc["bus"][bus_lookup[ppc_index], BUS_TYPE] == 4] = -1
+    return bus, ppc_index
 
-def _extract_results(net, ppc, ppc_0):
-    _get_bus_results(net, ppc, ppc_0)
+def _extract_results(net, ppc, ppc_0, bus):
+    _get_bus_results(net, ppc, ppc_0, bus)
     if net._options["branch_results"]:
         if net._options['return_all_currents']:
-            _get_line_all_results(net, ppc)
-            _get_trafo_all_results(net, ppc)
-            _get_trafo3w_all_results(net, ppc)
+            _get_line_all_results(net, ppc, bus)
+            _get_trafo_all_results(net, ppc, bus)
+            _get_trafo3w_all_results(net, ppc, bus)
         else:
             _get_line_results(net, ppc)
             _get_trafo_results(net, ppc)
             _get_trafo3w_results(net, ppc)
 
-
 def _extract_single_results(net, ppc):
     for element in ["line", "trafo"]:
         net["res_%s_sc"%element] = pd.DataFrame(np.nan, index=net[element].index,
@@ -51,7 +58,10 @@ def _get_single_bus_results(net, ppc):
     net["res_bus_sc"]["va_degree"] = ppc["bus"][bus_idx, VA]
 
 
-def _get_bus_results(net, ppc, ppc_0):
+def _get_bus_results(net, ppc, ppc_0, bus):
+    if bus is None:
+        bus = slice(None)
+
     bus_lookup = net._pd2ppc_lookups["bus"]
     ppc_index = bus_lookup[net.bus.index]
 
@@ -64,6 +74,8 @@ def _get_bus_results(net, ppc, ppc_0):
     if net._options["ith"]:
         net.res_bus_sc["ith_ka"] = ppc["bus"][ppc_index, ITH]
 
+    net.res_bus_sc = net.res_bus_sc.loc[bus, :]
+
 
 def _get_line_results(net, ppc):
     branch_lookup = net._pd2ppc_lookups["branch"]
@@ -78,13 +90,13 @@ def _get_line_results(net, ppc):
             net.res_line_sc["ith_ka"] = minmax(ppc["branch"][f:t, [ITH_F, ITH_T]].real, axis=1)
 
 
-def _get_line_all_results(net, ppc):
+def _get_line_all_results(net, ppc, bus):
     case = net._options["case"]
-    bus_lookup = net._pd2ppc_lookups["bus"]
-    ppc_index = bus_lookup[net.bus.index]
+
+    bus, ppc_index = _get_bus_ppc_idx_for_br_all_results(net, ppc, bus)
     branch_lookup = net._pd2ppc_lookups["branch"]
 
-    multindex = pd.MultiIndex.from_product([net.res_line_sc.index, net.bus.index], names=['line','bus'])
+    multindex = pd.MultiIndex.from_product([net.res_line_sc.index, bus], names=['line','bus'])
     net.res_line_sc = net.res_line_sc.reindex(multindex)
 
     if "line" in branch_lookup:
@@ -100,7 +112,6 @@ def _get_line_all_results(net, ppc):
             net.res_line_sc["ith_ka"] = minmax(ppc["internal"]["branch_ith_f"][f:t, ppc_index].real.reshape(-1, 1),
                                                ppc["internal"]["branch_ith_t"][f:t, ppc_index].real.reshape(-1, 1))
 
-
 def _get_trafo_results(net, ppc):
     branch_lookup = net._pd2ppc_lookups["branch"]
     if "trafo" in branch_lookup:
@@ -109,12 +120,11 @@ def _get_trafo_results(net, ppc):
         net.res_trafo_sc["ikss_lv_ka"] = ppc["branch"][f:t, IKSS_T].real
 
 
-def _get_trafo_all_results(net, ppc):
-    bus_lookup = net._pd2ppc_lookups["bus"]
-    ppc_index = bus_lookup[net.bus.index]
+def _get_trafo_all_results(net, ppc, bus):
+    bus, ppc_index = _get_bus_ppc_idx_for_br_all_results(net, ppc, bus)
     branch_lookup = net._pd2ppc_lookups["branch"]
 
-    multindex = pd.MultiIndex.from_product([net.res_trafo_sc.index, net.bus.index], names=['trafo', 'bus'])
+    multindex = pd.MultiIndex.from_product([net.res_trafo_sc.index, bus], names=['trafo', 'bus'])
     net.res_trafo_sc = net.res_trafo_sc.reindex(multindex)
 
     if "trafo" in branch_lookup:
@@ -135,12 +145,11 @@ def _get_trafo3w_results(net, ppc):
         net.res_trafo3w_sc["ikss_lv_ka"] = ppc["branch"][mv:lv, IKSS_T].real
 
 
-def _get_trafo3w_all_results(net, ppc):
-    bus_lookup = net._pd2ppc_lookups["bus"]
-    ppc_index = bus_lookup[net.bus.index]
+def _get_trafo3w_all_results(net, ppc, bus):
+    bus, ppc_index = _get_bus_ppc_idx_for_br_all_results(net, ppc, bus)
     branch_lookup = net._pd2ppc_lookups["branch"]
 
-    multindex = pd.MultiIndex.from_product([net.res_trafo3w_sc.index, net.bus.index], names=['trafo3w', 'bus'])
+    multindex = pd.MultiIndex.from_product([net.res_trafo3w_sc.index, bus], names=['trafo3w', 'bus'])
     net.res_trafo3w_sc = net.res_trafo3w_sc.reindex(multindex)
 
     if "trafo3w" in branch_lookup:


Test patch:
diff --git a/pandapower/test/shortcircuit/test_1ph.py b/pandapower/test/shortcircuit/test_1ph.py
index da4dab22..0bd8cd3a 100644
--- a/pandapower/test/shortcircuit/test_1ph.py
+++ b/pandapower/test/shortcircuit/test_1ph.py
@@ -89,39 +89,68 @@ def test_1ph_shortcircuit_min():
                 ,"Dd":  [0.52209346201, 0.66632662571, 0.66756160176, 0.72517293174]
                }
 
-    net = pp.create_empty_network()
-    for vc in results.keys():
-         add_network(net, vc)
-         try:
-             sc.calc_sc(net, fault="1ph", case="min")
-         except:
-             raise UserWarning("Did not converge after adding transformer with vector group %s"%vc)
-
-    for vc, result in results.items():
-        check_results(net, vc, result)
+    for inv_y in (False, True):
+        net = pp.create_empty_network()
+        for vc in results.keys():
+             add_network(net, vc)
+             try:
+                 sc.calc_sc(net, fault="1ph", case="min", inverse_y=inv_y)
+             except:
+                 raise UserWarning("Did not converge after adding transformer with vector group %s"%vc)
+    
+        for vc, result in results.items():
+            check_results(net, vc, result)
 
 def test_iec60909_example_4():
     file = os.path.join(pp.pp_dir, "test", "test_files", "IEC60909-4_example.json")
     net = pp.from_json(file)
-    sc.calc_sc(net, fault="1ph")
-    assert np.isclose(net.res_bus_sc[net.bus.name=="Q"].ikss_ka.values[0], 10.05957231)
-    assert np.isclose(net.res_bus_sc[net.bus.name=="T2LV"].ikss_ka.values[0], 34.467353142)
-    assert np.isclose(net.res_bus_sc[net.bus.name=="F1"].ikss_ka.values[0], 35.53066312)
-    assert np.isclose(net.res_bus_sc[net.bus.name=="F2"].ikss_ka.values[0], 34.89135137)
-    assert np.isclose(net.res_bus_sc[net.bus.name=="F3"].ikss_ka.values[0], 5.0321033105)
-    assert np.isclose(net.res_bus_sc[net.bus.name=="Cable/Line IC"].ikss_ka.values[0], 16.362586813)
+    for inv_y in (False, True):
+        sc.calc_sc(net, fault="1ph", inverse_y=inv_y)
+        assert np.isclose(net.res_bus_sc[net.bus.name=="Q"].ikss_ka.values[0], 10.05957231)
+        assert np.isclose(net.res_bus_sc[net.bus.name=="T2LV"].ikss_ka.values[0], 34.467353142)
+        assert np.isclose(net.res_bus_sc[net.bus.name=="F1"].ikss_ka.values[0], 35.53066312)
+        assert np.isclose(net.res_bus_sc[net.bus.name=="F2"].ikss_ka.values[0], 34.89135137)
+        assert np.isclose(net.res_bus_sc[net.bus.name=="F3"].ikss_ka.values[0], 5.0321033105)
+        assert np.isclose(net.res_bus_sc[net.bus.name=="Cable/Line IC"].ikss_ka.values[0], 16.362586813)
+
+def test_iec60909_example_4_bus_selection():
+    file = os.path.join(pp.pp_dir, "test", "test_files", "IEC60909-4_example.json")
+    net = pp.from_json(file)
+    for inv_y in (False, True):
+        sc.calc_sc(net, fault="1ph", inverse_y=inv_y, 
+                   bus=net.bus[net.bus.name.isin(("F1", "F2"))].index)
+        assert np.isclose(net.res_bus_sc.at[net.bus[net.bus.name=="F1"].index[0],
+                                            "ikss_ka"], 35.53066312)
+        assert np.isclose(net.res_bus_sc.at[net.bus[net.bus.name=="F2"].index[0],
+                                            "ikss_ka"], 34.89135137)
+
+def test_iec60909_example_4_bus_selection_br_res():
+    file = os.path.join(pp.pp_dir, "test", "test_files", "IEC60909-4_example.json")
+    net = pp.from_json(file)
+    for inv_y in (False, True):
+        sc.calc_sc(net, fault="1ph", inverse_y=inv_y,
+                   bus=net.bus[net.bus.name.isin(("F1", "F2"))].index,
+                   branch_results=True)
+        sc.calc_sc(net, fault="1ph", inverse_y=inv_y,
+                   bus=net.bus[net.bus.name.isin(("F1", "F2"))].index,
+                   branch_results=True, return_all_currents=True)
+        assert np.isclose(net.res_bus_sc.at[net.bus[net.bus.name=="F1"].index[0],
+                                            "ikss_ka"], 35.53066312)
+        assert np.isclose(net.res_bus_sc.at[net.bus[net.bus.name=="F2"].index[0],
+                                            "ikss_ka"], 34.89135137)
 
 def test_1ph_with_switches():
-    net = pp.create_empty_network()
-    vc = "Yy"
-    l1, l2, _ = add_network(net, vc)
-    sc.calc_sc(net, fault="1ph", case="max")
-    pp.create_line(net, net.line.to_bus.at[l2], net.line.from_bus.at[l1], length_km=15,
-                   std_type="unsymmetric_line_type", parallel=2.)
-    pp.add_zero_impedance_parameters(net)
-    pp.create_switch(net, bus=net.line.to_bus.at[l2], element=l2, et="l", closed=False)
-    sc.calc_sc(net, fault="1ph", case="max")
-    check_results(net, vc, [0.52209347338, 2.0620266652, 2.3255761263, 2.3066467489])
+    for inv_y in (False, True):
+        net = pp.create_empty_network()
+        vc = "Yy"
+        l1, l2, _ = add_network(net, vc)
+        sc.calc_sc(net, fault="1ph", case="max", inverse_y=inv_y)
+        pp.create_line(net, net.line.to_bus.at[l2], net.line.from_bus.at[l1], length_km=15,
+                       std_type="unsymmetric_line_type", parallel=2.)
+        pp.add_zero_impedance_parameters(net)
+        pp.create_switch(net, bus=net.line.to_bus.at[l2], element=l2, et="l", closed=False)
+        sc.calc_sc(net, fault="1ph", case="max")
+        check_results(net, vc, [0.52209347338, 2.0620266652, 2.3255761263, 2.3066467489])
 
 if __name__ == "__main__":
-    pytest.main([__file__])
\ No newline at end of file
+    pytest.main([__file__])
diff --git a/pandapower/test/shortcircuit/test_all_currents.py b/pandapower/test/shortcircuit/test_all_currents.py
index 845bda93..b73198c0 100644
--- a/pandapower/test/shortcircuit/test_all_currents.py
+++ b/pandapower/test/shortcircuit/test_all_currents.py
@@ -197,12 +197,14 @@ def test_all_currents_with_oos_elements(three_bus_example):
 
     assert np.allclose(net.res_line_sc.ikss_ka.loc[[(0, 0), (0, 1)]].values,
                        np.array([0.01259673, 0.49593036]), atol=1e-5)
-    assert all(net.res_line_sc.ikss_ka.loc[[(0, 2), (1, 0), (1, 1), (1, 2)]].isnull())
+    assert np.allclose(net.res_line_sc.ikss_ka.loc[[(0, 2), (1, 0), (1, 1), (1, 2)]].values,
+                       0, atol=1e-10)
 
     sc.calc_sc(net, case="min", branch_results=True, return_all_currents=True)
     assert np.allclose(net.res_line_sc.ikss_ka.loc[[(0, 0), (0, 1)]].values,
                        np.array([0.01259673, 0.3989686]), atol=1e-5)
-    assert all(net.res_line_sc.ikss_ka.loc[[(0, 2), (1, 0), (1, 1), (1, 2)]].isnull())
+    assert np.allclose(net.res_line_sc.ikss_ka.loc[[(0, 2), (1, 0), (1, 1), (1, 2)]].values,
+                       0, atol=1e-10)
 
 
 def test_branch_all_currents_gen(gen_three_bus_example):
@@ -245,7 +247,7 @@ def test_against_single_sc_results_line(three_bus_permuted_index):
 
 def test_against_single_sc_results_trafo(net_transformer):
     net = net_transformer
-    sc.calc_sc(net, case="max", branch_results=True, return_all_currents=True)
+    sc.calc_sc(net, case="max", branch_results=True, return_all_currents=True, inverse_y=False)
     multi_results = net.res_trafo_sc.copy()
 
     for bus in net.bus.index[net.bus.in_service]:
diff --git a/pandapower/test/shortcircuit/test_gen.py b/pandapower/test/shortcircuit/test_gen.py
index 742fcb23..733826d3 100644
--- a/pandapower/test/shortcircuit/test_gen.py
+++ b/pandapower/test/shortcircuit/test_gen.py
@@ -45,7 +45,7 @@ def gen_three_bus_example():
 
 def test_max_gen(one_line_one_generator):
     net = one_line_one_generator
-    sc.calc_sc(net, case="max")
+    sc.calc_sc(net, case="max", inverse_y=False)
     assert abs(net.res_bus_sc.ikss_ka.at[0] - 1.5395815) < 1e-7
     assert abs(net.res_bus_sc.ikss_ka.at[2] - 1.5395815) < 1e-7
     assert abs(net.res_bus_sc.ikss_ka.at[1] - 1.5083952) < 1e-7
@@ -144,7 +144,6 @@ def test_close_to_gen_simple2():
     # sc.calc_single_sc(net, b2)
     sc.calc_sc(net, tk_s=5e-2)
 
-
 def test_generator_book():
     net=pp.create_empty_network()
     b1= pp.create_bus(net, 110)
@@ -156,7 +155,6 @@ def test_generator_book():
     pp.create_gen(net, b2, 0, 1, sn_mva=25, vn_kv=6.3, xdss_pu=0.11, cos_phi=np.cos(np.arcsin(0.8)))
     sc.calc_sc(net, tk_s=2.5e-2)
 
-
 def test_shunt():
     net=pp.create_empty_network()
     b1= pp.create_bus(net, 110)
@@ -168,7 +166,6 @@ def test_shunt():
 
     sc.calc_sc(net, tk_s=2.5e-2)
 
-
 def test_power_station_unit():
     net = pp.create_empty_network()
     b1, b2, b3, b4 = pp.create_buses(net, 4, 20)
@@ -185,21 +182,5 @@ def test_power_station_unit():
 
     sc.calc_sc(net)
 
-
-
-
 if __name__ == '__main__':
     pytest.main(['test_gen.py'])
-#    net = pp.create_empty_network(sn_mva=2)
-#    b1 = pp.create_bus(net, vn_kv=10.)
-#    b2 = pp.create_bus(net, vn_kv=10.)
-#    b3 = pp.create_bus(net, vn_kv=10.)
-#    pp.create_bus(net, vn_kv=0.4, in_service=False)
-#    pp.create_gen(net, b1, vn_kv=10.5, xdss_pu=0.2, rdss_pu=0.001, cos_phi=0.8, p_mw=0.1, sn_mva=2.5)
-#    pp.create_gen(net, b1, vn_kv=10.5, xdss_pu=0.2, rdss_pu=0.001, cos_phi=0.8, p_mw=0.1, sn_mva=2.5)
-#    l = pp.create_line_from_parameters(net, b2, b1, length_km=1.0, max_i_ka=0.29,
-#                                       r_ohm_per_km=0.1548, x_ohm_per_km=0.0816814, c_nf_per_km=165)
-#    net.line.loc[l, "endtemp_degree"] = 165
-#    pp.create_switch(net, b3, b1, et="b")
-#    sc.calc_sc(net)
-#    pp.reset_results(net, "sc")
\ No newline at end of file
diff --git a/pandapower/test/shortcircuit/test_iec60909_4.py b/pandapower/test/shortcircuit/test_iec60909_4.py
new file mode 100644
index 00000000..6c391411
--- /dev/null
+++ b/pandapower/test/shortcircuit/test_iec60909_4.py
@@ -0,0 +1,249 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (c) 2016-2017 by University of Kassel and Fraunhofer Institute for Wind Energy and
+# Energy System Technology (IWES), Kassel. All rights reserved. Use of this source code is governed
+# by a BSD-style license that can be found in the LICENSE file.
+
+import numpy as np
+import pytest
+
+import pandapower as pp
+import pandapower.shortcircuit as sc
+
+@pytest.fixture
+def iec_60909_4():
+    net = pp.create_empty_network()
+
+    b1 = pp.create_bus(net, vn_kv=380.)
+    b2 = pp.create_bus(net, vn_kv=110.)
+    b3 = pp.create_bus(net, vn_kv=110.)
+    b4 = pp.create_bus(net, vn_kv=110.)
+    b5 = pp.create_bus(net, vn_kv=110.)
+    b6 = pp.create_bus(net, vn_kv=10.)
+    b7 = pp.create_bus(net, vn_kv=10.)
+    b8 = pp.create_bus(net, vn_kv=30.)
+    H = pp.create_bus(net, vn_kv=30.)
+    HG1 = pp.create_bus(net, vn_kv=21.)
+    HG2 = pp.create_bus(net, vn_kv=10.5)  # 10.5kV?
+    T_T5 = pp.create_bus(net, vn_kv=10.5)
+    T_T6 = pp.create_bus(net, vn_kv=10.5)
+
+    pp.create_ext_grid(net, b1, s_sc_max_mva=38 * 380 * np.sqrt(3), rx_max=0.1)
+    pp.create_ext_grid(net, b5, s_sc_max_mva=16 * 110 * np.sqrt(3), rx_max=0.1)
+
+    pp.create_gen(net, HG2, p_mw=0.9 * 100, vn_kv=10.5, xdss_pu=0.16, rdss_pu=0.005, cos_phi=0.9, sn_mva=100)
+    pp.create_gen(net, HG1, p_mw=0.85 * 150, vn_kv=21, xdss_pu=0.14, rdss_pu=0.002, cos_phi=0.85, sn_mva=150)
+    pp.create_gen(net, b6, p_mw=0.8 * 10, vn_kv=10.5, xdss_pu=0.1, rdss_pu=0.018, cos_phi=0.8, sn_mva=10)
+
+    pp.create_transformer_from_parameters(net, b4, HG1, sn_mva=150,
+        pfe_kw=0, i0_percent=0,  # FIXME: Optional for SC
+        vn_hv_kv=115., vn_lv_kv=21, vk_percent=16, vkr_percent=0.5)
+    pp.create_transformer_from_parameters(net, b3, HG2, sn_mva=100,
+        pfe_kw=0, i0_percent=0,
+        vn_hv_kv=120., vn_lv_kv=10.5, vk_percent=12, vkr_percent=0.5)
+
+    pp.create_transformer3w_from_parameters(net,
+        hv_bus=b1, mv_bus=b2, lv_bus=H,
+        vn_hv_kv=400, vn_mv_kv=120, vn_lv_kv=30,
+        sn_hv_mva=350, sn_mv_mva=350, sn_lv_mva=50,
+        pfe_kw=0, i0_percent=0,  # FIXME: Optional for SC
+        vk_hv_percent=21, vkr_hv_percent=.26,
+        vk_mv_percent=10, vkr_mv_percent=.16,
+        vk_lv_percent=7., vkr_lv_percent=.16)
+    pp.create_transformer3w_from_parameters(net,
+        hv_bus=b1, mv_bus=b2, lv_bus=b8,
+        vn_hv_kv=400, vn_mv_kv=120, vn_lv_kv=30,
+        sn_hv_mva=350, sn_mv_mva=350, sn_lv_mva=50,
+        pfe_kw=0, i0_percent=0,
+        vk_hv_percent=21, vkr_hv_percent=.26,
+        vk_mv_percent=10, vkr_mv_percent=.16,
+        vk_lv_percent=7., vkr_lv_percent=.16)
+
+    pp.create_transformer3w_from_parameters(net,
+        hv_bus=b5, mv_bus=b6, lv_bus=T_T5,
+        vn_hv_kv=115., vn_mv_kv=10.5, vn_lv_kv=10.5,
+        sn_hv_mva=31.5, sn_mv_mva=31.5, sn_lv_mva=31.5,
+        pfe_kw=0, i0_percent=0,  # FIXME: Optional for SC
+        vk_hv_percent=12, vkr_hv_percent=.5,
+        vk_mv_percent=12, vkr_mv_percent=.5,
+        vk_lv_percent=12, vkr_lv_percent=.5)
+    pp.create_transformer3w_from_parameters(net,
+        hv_bus=b5, mv_bus=b6, lv_bus=T_T6,
+        vn_hv_kv=115., vn_mv_kv=10.5, vn_lv_kv=10.5,
+        sn_hv_mva=31.5, sn_mv_mva=31.5, sn_lv_mva=31.5,
+        pfe_kw=0, i0_percent=0,  # FIXME: Optional for SC
+        vk_hv_percent=12, vkr_hv_percent=.5,
+        vk_mv_percent=12, vkr_mv_percent=.5,
+        vk_lv_percent=12, vkr_lv_percent=.5)
+
+    # pp.create_transformer_from_parameters(net, b5, b6, sn_mva=31.5,
+    #     pfe_kw=0, i0_percent=0,  # FIXME: Optional for SC
+    #     vn_hv_kv=115., vn_lv_kv=10.5, vk_percent=12, vkr_percent=0.5)
+    # pp.create_transformer_from_parameters(net, b5, b6, sn_mva=31.5,
+    #     pfe_kw=0, i0_percent=0,  # FIXME: Optional for SC
+    #     vn_hv_kv=115., vn_lv_kv=10.5, vk_percent=12, vkr_percent=0.5)
+
+    # First try without motor in powerfactory
+    # pp.create_motor(net, b7, p_mw=5.0, sn_mva=5.0 / 0.88,  k=5, rx=0.1)  # FIXME: R/X
+    # pp.create_motor(net, b7, p_mw=2.0, sn_mva=2.0 / 0.89, k=5.2, rx=0.1)
+    # pp.create_motor(net, b7, p_mw=2.0, sn_mva=2.0 / 0.89, k=5.2, rx=0.1)
+
+    pp.create_line_from_parameters(net, b2, b3, name="L1",
+        c_nf_per_km=0, max_i_ka=0,  # FIXME: Optional for SC
+        length_km=20, r_ohm_per_km=0.12, x_ohm_per_km=0.39)
+    pp.create_line_from_parameters(net, b3, b4, name="L2",
+        c_nf_per_km=0, max_i_ka=0,
+        length_km=10, r_ohm_per_km=0.12, x_ohm_per_km=0.39)
+    pp.create_line_from_parameters(net, b2, b5, name="L3a",
+        c_nf_per_km=0, max_i_ka=0,
+        length_km=5, r_ohm_per_km=0.12, x_ohm_per_km=0.39)
+    pp.create_line_from_parameters(net, b2, b5, name="L3b",
+        c_nf_per_km=0, max_i_ka=0,
+        length_km=5, r_ohm_per_km=0.12, x_ohm_per_km=0.39)
+    pp.create_line_from_parameters(net, b5, b3, name="L4",
+        c_nf_per_km=0, max_i_ka=0,
+        length_km=10, r_ohm_per_km=0.096, x_ohm_per_km=0.388)
+    pp.create_line_from_parameters(net, b5, b4, name="L5",
+        c_nf_per_km=0, max_i_ka=0,
+        length_km=15, r_ohm_per_km=0.12, x_ohm_per_km=0.39)
+    pp.create_line_from_parameters(net, b6, b7, name="L6",
+        c_nf_per_km=0, max_i_ka=0,
+        length_km=1, r_ohm_per_km=0.12, x_ohm_per_km=0.39)
+
+    return net
+
+@pytest.fixture
+def iec_60909_4_small():
+    net = pp.create_empty_network()
+
+    b1 = pp.create_bus(net, vn_kv=380.)
+    b2 = pp.create_bus(net, vn_kv=110.)
+    b3 = pp.create_bus(net, vn_kv=110.)
+    b5 = pp.create_bus(net, vn_kv=110.)
+    b8 = pp.create_bus(net, vn_kv=30.)
+    H = pp.create_bus(net, vn_kv=30.)
+    HG1 = pp.create_bus(net, vn_kv=21.)
+    HG2 = pp.create_bus(net, vn_kv=10)  # 10.5kV?
+    T_T5 = pp.create_bus(net, vn_kv=10.5)
+    T_T6 = pp.create_bus(net, vn_kv=10.5)
+
+    pp.create_ext_grid(net, b1, s_sc_max_mva=38 * 380 * np.sqrt(3), rx_max=0.1)
+    pp.create_ext_grid(net, b5, s_sc_max_mva=16 * 110 * np.sqrt(3), rx_max=0.1)
+
+    pp.create_gen(net, HG2, p_mw=0.9 * 100, vn_kv=10.5,
+                  xdss_pu=0.16, rdss_pu=0.005, cos_phi=0.9, sn_mva=100)
+
+    pp.create_transformer_from_parameters(net, b3, HG2, sn_mva=100,
+        pfe_kw=0, i0_percent=0,
+        vn_hv_kv=120., vn_lv_kv=10.5, vk_percent=12, vkr_percent=0.5)
+
+    pp.create_transformer3w_from_parameters(net,
+        hv_bus=b1, mv_bus=b2, lv_bus=H,
+        vn_hv_kv=400, vn_mv_kv=120, vn_lv_kv=30,
+        sn_hv_mva=350, sn_mv_mva=350, sn_lv_mva=50,
+        pfe_kw=0, i0_percent=0,  # FIXME: Optional for SC
+        vk_hv_percent=21, vkr_hv_percent=.26,
+        vk_mv_percent=7, vkr_mv_percent=.16,
+        vk_lv_percent=10., vkr_lv_percent=.16)
+    pp.create_transformer3w_from_parameters(net,
+        hv_bus=b1, mv_bus=b2, lv_bus=b8,
+        vn_hv_kv=400, vn_mv_kv=120, vn_lv_kv=30,
+        sn_hv_mva=350, sn_mv_mva=350, sn_lv_mva=50,
+        pfe_kw=0, i0_percent=0,
+        vk_hv_percent=21, vkr_hv_percent=.26,
+        vk_mv_percent=7, vkr_mv_percent=.16,
+        vk_lv_percent=10., vkr_lv_percent=.16)
+
+    pp.create_line_from_parameters(net, b2, b3, name="L1",
+        c_nf_per_km=0, max_i_ka=0,  # FIXME: Optional for SC
+        length_km=20, r_ohm_per_km=0.12, x_ohm_per_km=0.39,)
+    pp.create_line_from_parameters(net, b2, b5, name="L3a",
+        c_nf_per_km=0, max_i_ka=0,
+        length_km=5, r_ohm_per_km=0.12, x_ohm_per_km=0.39)
+    pp.create_line_from_parameters(net, b2, b5, name="L3b",
+        c_nf_per_km=0, max_i_ka=0,
+        length_km=5, r_ohm_per_km=0.12, x_ohm_per_km=0.39)
+    pp.create_line_from_parameters(net, b5, b3, name="L4",
+        c_nf_per_km=0, max_i_ka=0,
+        length_km=10, r_ohm_per_km=0.096, x_ohm_per_km=0.388)
+
+    return net
+
+
+def iec_60909_4_small_gen_only():
+    net = pp.create_empty_network()
+
+    b3 = pp.create_bus(net, vn_kv=110.)
+    HG2 = pp.create_bus(net, vn_kv=10)  # 10.5kV?
+
+    pp.create_gen(net, HG2, p_mw=0.9 * 100, vn_kv=10.5,
+                  xdss_pu=0.16, rdss_pu=0.00453, cos_phi=0.9, sn_mva=100, slack=True)
+
+    pp.create_transformer_from_parameters(net, b3, HG2, sn_mva=100,
+        pfe_kw=0, i0_percent=0,
+        vn_hv_kv=120., vn_lv_kv=10.5, vk_percent=12, vkr_percent=0.5)
+    return net
+
+
+def test_iec_60909_4_3ph_small_without_gen(iec_60909_4_small):
+    net = iec_60909_4_small
+    # Deactivate all gens
+    net.gen = net.gen.iloc[0:0, :]
+    
+    sc.calc_sc(net, fault="3ph", case="max", ip=True, tk_s=0.1, kappa_method="C")
+    ikss_pf = [40.3390,	28.4130, 14.2095, 28.7195, 13.4191]
+    ip_pf = [99.7374, 72.6580, 32.1954, 72.1443, 36.5036]
+    
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:5], np.array(ikss_pf), atol=1e-4)
+    assert np.allclose(net.res_bus_sc.ip_ka.values[:5], np.array(ip_pf), atol=1e-4)
+
+# WIP
+@pytest.mark.skip("Gens are not supported yet")
+def test_iec_60909_4_3ph_small_with_gen(iec_60909_4_small):
+    net = iec_60909_4_small
+    
+    sc.calc_sc(net, fault="3ph", case="max", ip=True, tk_s=0.1, kappa_method="C")
+    ikss_pf = [40.4754, 29.8334, 16.1684, 30.3573]
+    ip_pf = [100.1164, 76.1134, 37.3576, 76.2689]
+    ib_pf = [40.4754, 29.7337, 15.9593, 30.2245]
+    
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:4], np.array(ikss_pf), atol=1e-4)
+    assert np.allclose(net.res_bus_sc.ip_ka.values[:4], np.array(ip_pf), atol=1e-4)
+
+@pytest.mark.skip("Gens are not supported yet")
+def test_iec_60909_4_3ph_small_gen_only(iec_60909_4_small_gen_only):
+    net = iec_60909_4_small_gen_only
+    
+    sc.calc_sc(net, fault="3ph", case="max", ip=True, tk_s=0.1, kappa_method="C")
+    ikss_pf = [1.9755, 39.5042]
+    ip_pf = [5.2316, 104.1085]
+    ib_pf = [1.6071, 27.3470]
+
+    assert np.allclose(net.res_bus_sc.ikss_ka[:2].values, np.array(ikss_pf), atol=1e-4)
+    assert np.allclose(net.res_bus_sc.ip_ka[:2].values, np.array(ip_pf), atol=1e-4)
+
+@pytest.mark.skip("Gens are not supported yet")
+def test_iec_60909_4_3ph_without_motor(iec_60909_4):
+    net = iec_60909_4
+    sc.calc_sc(net, fault="3ph", case="max", ip=True, tk_s=0.1, kappa_method="C")
+
+    ikss_pf = [40.6347, 31.6635, 19.6231, 16.1956, 32.9971, 34.3559, 22.2762, 13.5726]
+    ip_pf = [100.5427, 80.3509, 45.7157, 36.7855, 82,9406, 90.6143, 43.3826, 36,9103]
+    
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:8], np.array(ikss_pf), atol=1e-4)
+
+# def test_iec_60909_4_3ph(iec_60909_4):
+#     net = iec_60909_4
+#     sc.calc_sc(net, fault="3ph", case="max", ip=True, tk_s=0.1, kappa_method="C")
+
+#     ikss = [40.6447, 31.7831, 19.6730, 16.2277, 33.1894, 37.5629, 25.5895, 13.5778]
+#     ip = [100.5766, 80.8249, 45.8249, 36.8041, 83.6266, 99.1910, 51.3864, 36.9201]
+#     ib = [40.645, 31.570, 19.388, 16.017, 32.795, 34.028, 23.212, 13.578]
+
+#     assert abs(net.res_bus_sc.ikss_ka.at[0] - ikss[0]) < 1e-4
+
+#     assert np.allclose(net.res_bus_sc.ikss_ka.values[:8], np.array(ikss), atol=1e-4)
+#     assert np.allclose(net.res_bus_sc.ip.values[:8], np.array(ip), rtol=1e-4)
+
+if __name__ == '__main__':
+    pytest.main(["test_iec60909_4.py"])
diff --git a/pandapower/test/shortcircuit/test_meshing_detection.py b/pandapower/test/shortcircuit/test_meshing_detection.py
index 0571b1f7..272670e2 100644
--- a/pandapower/test/shortcircuit/test_meshing_detection.py
+++ b/pandapower/test/shortcircuit/test_meshing_detection.py
@@ -5,7 +5,7 @@
 
 
 import os
-
+import numpy as np
 import pytest
 
 import pandapower as pp
@@ -24,148 +24,65 @@ def meshed_grid():
 def test_max_10_meshed_grid(meshed_grid):
     net = meshed_grid
     sc.calc_sc(net, case='max', ip=True, ith=True, lv_tol_percent= 10., kappa_method="B")
-    assert (abs(net.res_bus_sc.ikss_ka.at[0] - 5.773503) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[1] - 14.82619) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[2] - 4.606440) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[3] - 4.068637) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[4] - 13.61509) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[5] - 2.812111) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[6] - 1.212288) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[7] - 1.525655) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[8] - 1.781087) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[9] - 1.568337) <1e-5)
-
-    assert (abs(net.res_bus_sc.ip_ka.at[0] - 14.256050) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[1] - 33.751300) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[2] - 6.759302) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[3] - 6.359403) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[4] - 26.49241) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[5] - 4.726619) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[6] - 2.015958) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[7] - 2.538654) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[8] - 2.576375) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[9] - 2.608065) <1e-5)
-
-    assert (abs(net.res_bus_sc.ith_ka.at[0] - 5.871191) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[1] - 14.97527) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[2] - 4.613454) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[3] - 4.077662) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[4] - 13.68449) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[5] - 2.820525) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[6] - 1.215770) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[7] - 1.530048) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[8] - 1.783442) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[9] - 1.572843) <1e-5)
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:10],
+                       [5.773503, 14.82619, 4.606440, 4.068637, 13.61509,
+                        2.812111, 1.212288, 1.525655, 1.781087, 1.568337], atol=1e-5)
+
+    assert np.allclose(net.res_bus_sc.ip_ka.values[:10],
+                       [14.256050, 33.751300, 6.759302, 6.359403, 26.49241,
+                        4.726619, 2.015958, 2.538654, 2.576375, 2.608065], atol=1e-5)
+
+    assert np.allclose(net.res_bus_sc.ith_ka.values[:10],
+                       [5.871191, 14.97527, 4.613454, 4.077662, 13.68449,
+                        2.820525, 1.215770, 1.530048, 1.783442, 1.572843], atol=1e-5)
 
 
 def test_max_6_meshed_grid(meshed_grid):
     net = meshed_grid
     sc.calc_sc(net, case='max', ip=True, ith=True, lv_tol_percent = 6., kappa_method="B")
-    assert (abs(net.res_bus_sc.ikss_ka.at[0] - 5.773503) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[1] - 14.75419) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[2] - 4.437882) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[3] - 4.068637) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[4] - 13.53425) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[5] - 2.701411) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[6] - 1.159945) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[7] - 1.460757) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[8] - 1.705172) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[9] - 1.501673) <1e-5)
-
-    assert (abs(net.res_bus_sc.ip_ka.at[0] - 14.25605) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[1] - 33.59996) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[2] - 6.50406) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[3] - 6.359403) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[4] - 26.28476) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[5] - 4.537759) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[6] - 1.928734) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[7] - 2.430331) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[8] - 2.466185) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[9] - 2.496901) <1e-5)
-
-    assert (abs(net.res_bus_sc.ith_ka.at[0] - 5.871191) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[1] - 14.90284) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[2] - 4.44457) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[3] - 4.077662) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[4] - 13.60275) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[5] - 2.709475) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[6] - 1.163276) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[7] - 1.464961) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[8] - 1.707423) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[9] - 1.505985) <1e-5)
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:10],
+                       [5.773503, 14.75419, 4.437882, 4.068637, 13.53425,
+                        2.701411, 1.159945, 1.460757, 1.705172, 1.501673], atol=1e-5)
+
+    assert np.allclose(net.res_bus_sc.ip_ka.values[:10],
+                       [14.25605, 33.59996, 6.50406, 6.359403, 26.28476,
+                        4.537759, 1.928734, 2.430331, 2.466185, 2.496901], atol=1e-5)
+
+    assert np.allclose(net.res_bus_sc.ith_ka.values[:10],
+                       [5.871191, 14.90284, 4.44457, 4.077662, 13.60275,
+                        2.709475, 1.163276, 1.464961, 1.707423, 1.505985], atol=1e-5)
 
 
 def test_min_10_meshed_grid(meshed_grid):
     net = meshed_grid
     sc.calc_sc(net, case='min', ip=True, ith=True, lv_tol_percent= 10., kappa_method="B")
-    assert (abs(net.res_bus_sc.ikss_ka.at[0] - 2.309401) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[1] - 11.3267) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[2] - 2.879343) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[3] - 1.884323) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[4] - 10.40083) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[5] - 1.693922) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[6] - 0.7107017) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[7] - 0.9000445) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[8] - 1.055881) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[9] - 0.928488) <1e-5)
-
-    assert (abs(net.res_bus_sc.ip_ka.at[0] - 5.702418) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[1] - 26.01655) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[2] - 4.166047) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[3] - 3.124163) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[4] - 20.04053) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[5] - 2.813883) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[6] - 1.179085) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[7] - 1.493293) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[8] - 1.523338) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[9] - 1.540432) <1e-5)
-
-    assert (abs(net.res_bus_sc.ith_ka.at[0] - 2.348476) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[1] - 11.44622) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[2] - 2.883161) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[3] - 1.889675) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[4] - 10.45195) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[5] - 1.698768) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[6] - 0.712725) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[7] - 0.9026074) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[8] - 1.057233) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[9] - 0.9311316) <1e-5)
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:10],
+                       [2.309401, 11.3267, 2.879343, 1.884323, 10.40083,
+                        1.693922, 0.7107017, 0.9000445, 1.055881, 0.928488], atol=1e-5)
+
+    assert np.allclose(net.res_bus_sc.ip_ka.values[:10],
+                       [5.702418, 26.01655, 4.166047, 3.124163, 20.04053,
+                        2.813883, 1.179085, 1.493293, 1.523338, 1.540432], atol=1e-5)
+
+    assert np.allclose(net.res_bus_sc.ith_ka.values[:10],
+                       [2.348476, 11.44622, 2.883161, 1.889675, 10.45195,
+                        1.698768, 0.712725, 0.9026074, 1.057233, 0.9311316], atol=1e-5)
+
 
 def test_min_6_meshed_grid(meshed_grid):
     net = meshed_grid
     sc.calc_sc(net, case='min', ip=True, ith=True, lv_tol_percent = 6., kappa_method="B")
-    assert (abs(net.res_bus_sc.ikss_ka.at[0] - 2.309401) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[1] - 11.75072) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[2] - 2.895465) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[3] - 1.884323) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[4] - 10.77961) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[5] - 1.700202) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[6] - 0.7116519) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[7] - 0.9016006) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[8] - 1.0576) <1e-5)
-    assert (abs(net.res_bus_sc.ikss_ka.at[9] - 0.9301236) <1e-5)
-
-    assert (abs(net.res_bus_sc.ip_ka.at[0] - 5.702418) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[1] - 27.00861) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[2] - 4.18812) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[3] - 3.124163) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[4] - 20.72881) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[5] - 2.824028) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[6] - 1.180654) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[7] - 1.495858) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[8] - 1.525799) <1e-5)
-    assert (abs(net.res_bus_sc.ip_ka.at[9] - 1.543131) <1e-5)
-
-    assert (abs(net.res_bus_sc.ith_ka.at[0] - 2.348476) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[1] - 11.87518) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[2] - 2.899291) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[3] - 1.889675) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[4] - 10.8322) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[5] - 1.705064) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[6] - 0.7136779) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[7] - 0.9041679) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[8] - 1.058954) <1e-5)
-    assert (abs(net.res_bus_sc.ith_ka.at[9] - 0.9327717) <1e-5)
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:10],
+                       [2.309401, 11.75072, 2.895465, 1.884323, 10.77961,
+                        1.700202, 0.7116519, 0.9016006, 1.0576, 0.9301236], atol=1e-5)
+
+    assert np.allclose(net.res_bus_sc.ip_ka.values[:10],
+                       [5.702418, 27.00861, 4.18812, 3.124163, 20.72881,
+                        2.824028, 1.180654, 1.495858, 1.525799, 1.543131], atol=1e-5)
+
+    assert np.allclose(net.res_bus_sc.ith_ka.values[:10],
+                       [2.348476, 11.87518, 2.899291, 1.889675, 10.8322,
+                        1.705064, 0.7136779, 0.9041679, 1.058954, 0.9327717], atol=1e-5)
 
 if __name__ == '__main__':
     pytest.main(["test_meshing_detection.py", '-xs'])
\ No newline at end of file
diff --git a/pandapower/test/shortcircuit/test_motor.py b/pandapower/test/shortcircuit/test_motor.py
index 0956bc89..10745274 100644
--- a/pandapower/test/shortcircuit/test_motor.py
+++ b/pandapower/test/shortcircuit/test_motor.py
@@ -33,16 +33,14 @@ def motor_net():
 def test_motor_min(motor_net):
     net = motor_net
     sc.calc_sc(net, case="min")
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[0], 11.547005315, rtol=1e-4)
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[1], 0.53709235574, rtol=1e-4)
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[2], 0.18070949061, rtol=1e-4)
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:3],
+                       [11.547005315, 0.53709235574, 0.18070949061], rtol=1e-4)
 
 def test_motor_max(motor_net):
     net = motor_net
     sc.calc_sc(net, case="max")
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[0], 14.743809197, rtol=1e-4)
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[1], 5.626994278, rtol=1e-4)
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[2], 0.370730612, rtol=1e-4)
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:3],
+                       [14.743809197, 5.626994278, 0.370730612], rtol=1e-4)
 
 def test_motor_max_branch(motor_net):
     net = motor_net
@@ -52,9 +50,8 @@ def test_motor_max_branch(motor_net):
 
     net.motor.in_service = True
     sc.calc_sc(net, case="max", branch_results=True)
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[0], 14.743809197, rtol=1e-4)
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[1], 5.626994278, rtol=1e-4)
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[2], 0.370730612, rtol=1e-4)
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:3],
+                       [14.743809197, 5.626994278, 0.370730612], rtol=1e-4)
     # The maximum current through the first branch is the short-circuit current
     # at the second bus without the motor contribution, which does not flow
     # through the line
@@ -77,9 +74,8 @@ def test_large_motor(motor_net):
     net = motor_net
     net.motor.pn_mech_mw = 10
     sc.calc_sc(net, case="max")
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[0], 14.695869025, rtol=1e-4)
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[1], 103.16722971, rtol=1e-4)
-    assert np.isclose(net.res_bus_sc.ikss_ka.at[2], 0.38693418116, rtol=1e-4)
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:3],
+                       [14.695869025, 103.16722971, 0.38693418116], rtol=1e-4)
 
 if __name__ == '__main__':
     pytest.main([__file__])
\ No newline at end of file
diff --git a/pandapower/test/shortcircuit/test_ring.py b/pandapower/test/shortcircuit/test_ring.py
index cfe21ad3..15011e69 100644
--- a/pandapower/test/shortcircuit/test_ring.py
+++ b/pandapower/test/shortcircuit/test_ring.py
@@ -30,9 +30,19 @@ def ring_network():
 
 def test_branch_results_open_ring(ring_network):
     net = ring_network
-    sc.calc_sc(net, branch_results=True)
+    sc.calc_sc(net, branch_results=True, inverse_y=False)
     assert np.allclose(net.res_trafo_sc.ikss_lv_ka.values, [0.47705988])
-    assert np.allclose(net.res_line_sc.ikss_ka.values, [0.45294928, 0.44514686, 0.47125418])
+    assert np.allclose(net.res_line_sc.ikss_ka.values, [0.45294928, 0.0, 0.47125418])
+    
+def test_branch_results_open_ring_with_impedance(ring_network):
+    net = ring_network
+    sc.calc_sc(net, branch_results=True, inverse_y=False)
+    res_line_no_imp = net.res_line_sc.ikss_ka.values.copy()
+
+    # Make sure that with fault impedance, the total current should be smaller
+    sc.calc_sc(net, branch_results=True, inverse_y=False, r_fault_ohm=1, x_fault_ohm=5)
+    non_null_flag = np.abs(res_line_no_imp) > 1e-10
+    assert np.all(net.res_line_sc.ikss_ka.values[non_null_flag] < res_line_no_imp[non_null_flag])
 
 def test_branch_results_closed_ring(ring_network):
     net = ring_network
@@ -45,7 +55,7 @@ def test_branch_results_closed_ring(ring_network):
 def test_kappa_methods(ring_network):
     net = ring_network
     net.switch.closed = True
-    sc.calc_sc(net, kappa_method="B", ip=True)
+    sc.calc_sc(net, kappa_method="B", ip=True, inverse_y=False)
     assert np.allclose(net.res_bus_sc.ip_ka.values,
                        [0.48810547956, 0.91192962511, 1.0264898716, 1.0360554521])
     sc.calc_sc(net, kappa_method="C", ip=True, topology="auto")
@@ -53,7 +63,4 @@ def test_kappa_methods(ring_network):
                        [0.48810547956, 0.91192962511, 0.89331396461, 0.90103415924])
 
 if __name__ == '__main__':
-#    net = ring_network()
-#    test_branch_results_open_ring(net)
     pytest.main(["test_ring.py"])
-
diff --git a/pandapower/test/shortcircuit/test_sc_multi_bus.py b/pandapower/test/shortcircuit/test_sc_multi_bus.py
new file mode 100644
index 00000000..1a3c20d1
--- /dev/null
+++ b/pandapower/test/shortcircuit/test_sc_multi_bus.py
@@ -0,0 +1,100 @@
+ # -*- coding: utf-8 -*-
+
+# Copyright (c) 2016-2020 by University of Kassel and Fraunhofer Institute for Energy Economics
+# and Energy System Technology (IEE), Kassel. All rights reserved.
+
+
+import os
+
+import pytest
+
+import numpy as np
+import pandapower as pp
+import pandapower.shortcircuit as sc
+from pandapower.test.shortcircuit.test_meshing_detection import meshed_grid
+
+@pytest.fixture
+def radial_grid():
+    net = pp.create_empty_network(sn_mva=2.)
+    b0 = pp.create_bus(net, 220)
+    b1 = pp.create_bus(net, 110)
+    b2 = pp.create_bus(net, 110)
+    b3 = pp.create_bus(net, 110)
+    pp.create_ext_grid(net, b0, s_sc_max_mva=100., s_sc_min_mva=80., rx_min=0.4, rx_max=0.4)
+    pp.create_transformer(net, b0, b1, "100 MVA 220/110 kV")
+    pp.create_line(net, b1, b2, std_type="305-AL1/39-ST1A 110.0" , length_km=20.)
+    pp.create_line(net, b2, b3, std_type="N2XS(FL)2Y 1x185 RM/35 64/110 kV" , length_km=15.)
+    return net
+
+@pytest.fixture
+def three_bus_big_sgen_example():
+     net = pp.create_empty_network()
+     b1 = pp.create_bus(net, 110)
+     b2 = pp.create_bus(net, 110)
+     b3 = pp.create_bus(net, 110)
+
+     pp.create_ext_grid(net, b1, s_sc_max_mva=100., s_sc_min_mva=80., rx_min=0.4, rx_max=0.4)
+     pp.create_line(net, b1, b2, std_type="305-AL1/39-ST1A 110.0" , length_km=20.)
+     pp.create_line(net, b2, b3, std_type="N2XS(FL)2Y 1x185 RM/35 64/110 kV" , length_km=15.)
+     net.line["endtemp_degree"] = 80
+
+     pp.create_sgen(net, b2, sn_mva=200., p_mw=0, k=1.2)
+     return net
+
+def test_radial_network(radial_grid):   
+    net = radial_grid
+    sc_bus = 3
+    sc.calc_sc(net)
+    ik = net.res_bus_sc.ikss_ka.at[sc_bus]
+    sc.calc_sc(net, bus=sc_bus, inverse_y=False, branch_results=True, return_all_currents=True)
+    assert np.isclose(net.res_bus_sc.ikss_ka.at[sc_bus], ik)
+    assert np.isclose(net.res_line_sc.ikss_ka.loc[(1, sc_bus)], ik)
+    assert np.isclose(net.res_line_sc.ikss_ka.loc[(0, sc_bus)], ik)
+    assert np.isclose(net.res_trafo_sc.ikss_lv_ka.loc[(0, sc_bus)], ik)
+    trafo_ratio = net.trafo.vn_lv_kv.values / net.trafo.vn_hv_kv.values
+    assert np.isclose(net.res_trafo_sc.ikss_hv_ka.loc[(0, sc_bus)], ik*trafo_ratio)
+
+    sc_bus = 2
+    sc.calc_sc(net)
+    ik = net.res_bus_sc.ikss_ka.at[sc_bus]
+    sc.calc_sc(net, bus=sc_bus, inverse_y=False, branch_results=True, return_all_currents=True)
+    assert np.isclose(net.res_bus_sc.ikss_ka.at[sc_bus], ik)
+    assert np.isclose(net.res_line_sc.ikss_ka.loc[(1, sc_bus)], 0)
+    assert np.isclose(net.res_line_sc.ikss_ka.loc[(0, sc_bus)], ik)
+    assert np.isclose(net.res_trafo_sc.ikss_lv_ka.loc[(0, sc_bus)], ik)
+    trafo_ratio = net.trafo.vn_lv_kv.values / net.trafo.vn_hv_kv.values
+    assert np.isclose(net.res_trafo_sc.ikss_hv_ka.loc[(0, sc_bus)], ik*trafo_ratio)
+
+def test_meshed_network(meshed_grid):
+    net = meshed_grid
+    sc.calc_sc(net)
+    sc_bus = 5
+    ik = net.res_bus_sc.ikss_ka.at[sc_bus]
+    
+    sc.calc_sc(net, bus=sc_bus, inverse_y=False, branch_results=True, return_all_currents=True)
+    assert np.isclose(net.res_bus_sc.ikss_ka.at[sc_bus], ik)
+    line_ix = net.line[(net.line.to_bus==sc_bus) | (net.line.from_bus==sc_bus)].index
+    line_flow_into_sc = net.res_line_sc.loc[(line_ix, sc_bus), "ikss_ka"].sum()
+    assert np.isclose(line_flow_into_sc, ik, atol=2e-3)
+
+def test_big_gen_network(three_bus_big_sgen_example):    
+    net = three_bus_big_sgen_example
+    sc_bus = [0, 1, 2]
+    sc.calc_sc(net, bus=sc_bus, branch_results=True, return_all_currents=True, inverse_y=False)
+    assert np.isclose(net.res_line_sc.loc[(0, 0),"ikss_ka"], 1.25967331, atol=1e-3)
+    assert np.isclose(net.res_line_sc.loc[(1, 0),"ikss_ka"], 0., atol=2e-3)
+    assert np.isclose(net.res_line_sc.loc[(0, 2),"ikss_ka"], 0.46221808, atol=1e-3)
+    assert np.isclose(net.res_line_sc.loc[(1, 2),"ikss_ka"], 1.72233192, atol=1e-3)
+
+def test_big_gen_network_v2(three_bus_big_sgen_example):    
+    net = three_bus_big_sgen_example
+    sc_bus = [0, 2]
+    sc.calc_sc(net, bus=sc_bus, branch_results=True, return_all_currents=True, inverse_y=False)
+    assert np.isclose(net.res_line_sc.loc[(0, 0),"ikss_ka"], 1.25967331, atol=1e-3)
+    assert np.isclose(net.res_line_sc.loc[(1, 0),"ikss_ka"], 0., atol=2e-3)
+    assert np.isclose(net.res_line_sc.loc[(0, 2),"ikss_ka"], 0.46221808, atol=1e-3)
+    assert np.isclose(net.res_line_sc.loc[(1, 2),"ikss_ka"], 1.72233192, atol=1e-3)
+
+
+if __name__ == '__main__':
+    pytest.main([__file__])
\ No newline at end of file
diff --git a/pandapower/test/shortcircuit/test_sc_single_bus.py b/pandapower/test/shortcircuit/test_sc_single_bus.py
index 7c1e543f..e10900f6 100644
--- a/pandapower/test/shortcircuit/test_sc_single_bus.py
+++ b/pandapower/test/shortcircuit/test_sc_single_bus.py
@@ -8,7 +8,7 @@ import os
 
 import pytest
 
-from numpy import isclose
+import numpy as np
 import pandapower as pp
 import pandapower.shortcircuit as sc
 from pandapower.test.shortcircuit.test_meshing_detection import meshed_grid
@@ -47,23 +47,24 @@ def test_radial_network(radial_grid):
     sc.calc_sc(net)
     ik = net.res_bus_sc.ikss_ka.at[sc_bus]
     sc.calc_single_sc(net, bus=sc_bus)
-    assert isclose(net.res_bus_sc.vm_pu.at[sc_bus], 0)
-    assert isclose(net.res_line_sc.i_ka.at[1], ik)
-    assert isclose(net.res_line_sc.i_ka.at[0], ik)
-    assert isclose(net.res_trafo_sc.i_lv_ka.at[0], ik)
+    assert np.isclose(net.res_bus_sc.vm_pu.at[sc_bus], 0)
+    assert np.isclose(net.res_line_sc.i_ka.at[1], ik)
+    assert np.isclose(net.res_line_sc.i_ka.at[0], ik)
+    assert np.isclose(net.res_trafo_sc.i_lv_ka.at[0], ik)
     trafo_ratio = net.trafo.vn_lv_kv.values / net.trafo.vn_hv_kv.values
-    assert isclose(net.res_trafo_sc.i_hv_ka.at[0], ik*trafo_ratio)
+    assert np.isclose(net.res_trafo_sc.i_hv_ka.at[0], ik*trafo_ratio)
     
-    sc_bus = 2
-    sc.calc_sc(net)
-    ik = net.res_bus_sc.ikss_ka.at[sc_bus]
-    sc.calc_single_sc(net, bus=sc_bus)
-    assert isclose(net.res_bus_sc.vm_pu.at[sc_bus], 0)
-    assert isclose(net.res_line_sc.i_ka.at[1], 0)
-    assert isclose(net.res_line_sc.i_ka.at[0], ik)
-    assert isclose(net.res_trafo_sc.i_lv_ka.at[0], ik)
-    trafo_ratio = net.trafo.vn_lv_kv.values / net.trafo.vn_hv_kv.values
-    assert isclose(net.res_trafo_sc.i_hv_ka.at[0], ik*trafo_ratio)
+    for inv_y in (False, True):
+        sc_bus = 2
+        sc.calc_sc(net, inverse_y=inv_y)
+        ik = net.res_bus_sc.ikss_ka.at[sc_bus]
+        sc.calc_single_sc(net, bus=sc_bus)
+        assert np.isclose(net.res_bus_sc.vm_pu.at[sc_bus], 0)
+        assert np.isclose(net.res_line_sc.i_ka.at[1], 0)
+        assert np.isclose(net.res_line_sc.i_ka.at[0], ik)
+        assert np.isclose(net.res_trafo_sc.i_lv_ka.at[0], ik)
+        trafo_ratio = net.trafo.vn_lv_kv.values / net.trafo.vn_hv_kv.values
+        assert np.isclose(net.res_trafo_sc.i_hv_ka.at[0], ik*trafo_ratio)
 
 
 def test_meshed_network(meshed_grid):
@@ -72,25 +73,60 @@ def test_meshed_network(meshed_grid):
     sc_bus = 5
     ik = net.res_bus_sc.ikss_ka.at[sc_bus]
     
-    sc.calc_single_sc(net, bus=sc_bus)
-    assert isclose(net.res_bus_sc.vm_pu.at[sc_bus], 0)
-    line_flow_into_sc = net.res_line_sc.i_ka[(net.line.to_bus==sc_bus) | (net.line.from_bus==sc_bus)].sum()
-    assert isclose(line_flow_into_sc, ik, atol=2e-3)
-
+    for inv_y in (False, True):
+        sc.calc_single_sc(net, bus=sc_bus, inverse_y=inv_y)
+        assert np.isclose(net.res_bus_sc.vm_pu.at[sc_bus], 0)
+        line_flow_into_sc = net.res_line_sc.i_ka[(net.line.to_bus==sc_bus) | (net.line.from_bus==sc_bus)].sum()
+        assert np.isclose(line_flow_into_sc, ik, atol=2e-3)
 
 def test_big_gen_network(three_bus_big_sgen_example):
     net = three_bus_big_sgen_example
+
     sc_bus = 0
-    sc.calc_single_sc(net, sc_bus)
+    sc.calc_single_sc(net, sc_bus, inverse_y=True)
+    assert np.isclose(net.res_line_sc.i_ka.at[0], 1.25967331, atol=1e-3)
+    assert np.isclose(net.res_line_sc.i_ka.at[1], 0., atol=1e-3)
 
-    assert isclose(net.res_line_sc.i_ka.at[0], 1.25967331, atol=1e-3)
-    assert isclose(net.res_line_sc.i_ka.at[1], 0., atol=2e-3)
+    sc_bus = 2
+    sc.calc_single_sc(net, sc_bus, inverse_y=True)
+    assert np.isclose(net.res_line_sc.i_ka.at[0], 0.46221808, atol=1e-3)
+    assert np.isclose(net.res_line_sc.i_ka.at[1], 1.72233192, atol=1e-3)
 
+@pytest.mark.xfail
+def test_big_gen_network_no_y_inv(three_bus_big_sgen_example):
+    # ("TODO: Fix the bug in the single sc mode without Yinv")
     net = three_bus_big_sgen_example
+
+    sc_bus = 0
+    sc.calc_single_sc(net, sc_bus, inverse_y=False)
+
+    assert np.isclose(net.res_line_sc.i_ka.at[0], 1.25967331, atol=1e-3)
+    assert np.isclose(net.res_line_sc.i_ka.at[1], 0., atol=1e-3)
+
     sc_bus = 2
-    sc.calc_single_sc(net, sc_bus)
-    assert isclose(net.res_line_sc.i_ka.at[0], 0.46221808, atol=1e-3)
-    assert isclose(net.res_line_sc.i_ka.at[1], 1.72233192, atol=1e-3)
+    sc.calc_single_sc(net, sc_bus, inverse_y=False)
+    assert np.isclose(net.res_line_sc.i_ka.at[0], 0.46221808, atol=1e-3)
+    assert np.isclose(net.res_line_sc.i_ka.at[1], 1.72233192, atol=1e-3)
+
+
+def test_big_gen_network_calc_sc(three_bus_big_sgen_example):
+    net = three_bus_big_sgen_example
+       
+    for inv_y in (False, True):
+        sc_bus = 0
+        sc.calc_sc(net, bus=sc_bus,
+                   branch_results=True, inverse_y=inv_y)
+    
+        assert np.isclose(net.res_line_sc.ikss_ka.at[0], 1.25967331, atol=1e-3)
+        assert np.isclose(net.res_line_sc.ikss_ka.at[1], 0., atol=1e-3)
+    
+        sc_bus = 2
+        sc.calc_sc(net, bus=sc_bus,
+                   branch_results=True, inverse_y=inv_y)
+        assert np.isclose(net.res_line_sc.ikss_ka.at[0], 0.46221808, atol=1e-3)
+        assert np.isclose(net.res_line_sc.ikss_ka.at[1], 1.72233192, atol=1e-3)
+
 
 if __name__ == '__main__':
-    pytest.main([__file__])
\ No newline at end of file
+    pytest.main([__file__])
+     
\ No newline at end of file
diff --git a/pandapower/test/shortcircuit/test_sgen.py b/pandapower/test/shortcircuit/test_sgen.py
index fc6b7f98..c56f82a8 100644
--- a/pandapower/test/shortcircuit/test_sgen.py
+++ b/pandapower/test/shortcircuit/test_sgen.py
@@ -82,8 +82,8 @@ def test_max_3ph_branch_small_sgen(three_bus_example):
     net = three_bus_example
     sc.calc_sc(net, case="max", ip=True, ith=True, branch_results=True)
     assert np.allclose(net.res_bus_sc.ikss_ka.values, np.array([0.53746061, 0.50852707, 0.4988896]))
-    assert np.allclose(net.res_line_sc.ikss_ka.values, np.array([ 0.49593034,  0.4988896 ]))
-    assert np.allclose(net.res_line_sc.ip_ka.values, np.array([ 0.92787443,  0.9251165 ]))
+    assert np.allclose(net.res_line_sc.ikss_ka.values, np.array([ 0.49593034,  0.4988896]))
+    assert np.allclose(net.res_line_sc.ip_ka.values, np.array([ 0.92787443,  0.9251165]))
     assert np.allclose(net.res_line_sc.ith_ka.values, np.array([ 0.49811957,  0.50106881]))
 
 
@@ -91,7 +91,7 @@ def test_max_3ph_branch_big_sgen(big_sgen_three_bus_example):
     net = big_sgen_three_bus_example
     sc.calc_sc(net, case="max", ip=True, ith=True, branch_results=True)
     assert np.allclose(net.res_bus_sc.ikss_ka.values, np.array([1.78453722, 1.75560368, 1.72233192]))
-    assert np.allclose(net.res_line_sc.ikss_ka.values, np.array([1.25967331, 1.72233192 ]))
+    assert np.allclose(net.res_line_sc.ikss_ka.values, np.array([1.25967331, 1.72233192]))
     assert np.allclose(net.res_line_sc.ip_ka.values, np.array([1.78144709, 2.65532524]))
     assert np.allclose(net.res_line_sc.ith_ka.values, np.array([1.26511638, 1.7298553]))