Skip to content

Analysis routines

Higher-level analysis workflows on top of a fully built network: the inverse rho problem at the bus-grounding side, the conductor thermal-limit check that assesses whether a grounding conductor survives the fault current it carries, and its node-side counterpart, the node thermal-limit check for earthing conductors and earth electrodes.

Physical / modelling context

A fully built Network solves a forward problem: given the bus grounding impedances \(Z_{\text{B},i}(\rho_E, f)\), find the EPR at the fault bus. The natural inverse question — "how poor can the soil get before the EPR exceeds a touch-voltage limit \(u_{\max}\)?" — is recurring in safety-engineering workflows: it sets the worst-case specific earth resistance the network can tolerate while still satisfying the relevant standard (e.g. EN 50522) at the assumed fault current.

find_max_rho_scaling solves that problem by log-bisection over a uniform scaling factor \(c\) of \(\rho_E\) on a user-supplied list of selected buses. At every trial \(c\), the bus-grounding impedance is re-evaluated through the bus's own BusType.impedance_formula and run_fault is invoked. The bracket converges on the largest \(c\) for which \(|U_\text{EPR}(f)|_{\text{RMS}} \le u_{\max}\) at the fault bus. The sister routine find_max_rho_f_scaling (and the diagnostic helper evaluate_max_epr_under_k) extend the same idea to a frequency- dependent rho-f characteristic: rho is scaled while a separate factor \(k\) tunes the imaginary, frequency-coupling part of the formula, so the inverse problem can be posed against a parametric rho-f curve rather than a single scalar.

The shape of the rho-f curve at each bus is controlled entirely by the existing BusType.impedance_formula; the analysis routines only vary the scalar factors, so any user-defined parametric form (linear, square root, frequency-dependent, …) is supported transparently. Original rho values are restored via a finally block, so the network is left untouched even if the bisection fails.

Example

import groundinsight as gi

# Assume `net` is a built network with at least one fault and one
# source. Pick the buses whose specific_earth_resistance should be
# scaled jointly — typically all buses sharing a soil environment.
result = gi.find_max_rho_scaling(
    network=net,
    fault_name="fault1",
    bus_names=["bus_substation", "bus_fault"],
    u_max=200.0,            # touch-voltage limit in volts (RMS)
    c_bounds=(0.1, 100.0),  # search bracket on the scaling factor
    tol_rel=1e-3,
    max_iter=40,
)

# `result` is a dict with keys: c_max, u_epr_rms_at_c_max,
# rho_max (per-bus dict), iterations, status, converged,
# c_bracket, bracket_rel_width.
#
# Read `converged` before `c_max`: c_max is always a scaling factor
# whose EPR was measured and found admissible, but only a converged
# search has also shown that nothing meaningfully larger is.
if not result["converged"]:
    print(f"not a maximum: {result['status']}, "
          f"bracket {result['c_bracket']}")

print(f"c_max = {result['c_max']:.3f}")
print(f"EPR at c_max = {result['u_epr_rms_at_c_max']:.1f} V")
for bus, rho_max in result["rho_max"].items():
    print(f"  {bus}: rho_max = {rho_max:.0f} Ω·m")

Reading the result

c_max carries a one-directional guarantee: the EPR at that scaling factor was computed and found to satisfy the limit. Whether it is also the largest such factor depends on how the search ended, and that is what status reports.

status converged what c_max means
"converged" True The bracket closed to within tol_rel. The true threshold lies inside c_bracket.
"bracket_within_tol_on_entry" True The supplied bracket was already narrower than tol_rel, so no step was needed. Same guarantee as above.
"bracket_fully_admissible" False Every factor in the bracket satisfied the limit, so c_max is the upper bound — a lower bound on the true maximum, since nothing above it was evaluated. c_bracket is (c_hi, inf), so math.isfinite(result["c_bracket"][1]) tests for this case. Widen c_bounds.
"max_iter_reached" False The step cap was hit before the bracket closed. c_max is admissible but can be far below the true threshold — on a saturating test network, a cap of 3 instead of 60 came out 81 % low. Raise max_iter, or narrow c_bounds.

iterations on its own does not identify a case: two of the four rows above can produce iterations == 0, and their c_max comes from opposite ends of the bracket. Branch on converged or status, never on the step count.

bracket_rel_width is (c_hi - c_lo) / c_lo for the reported bracket — directly comparable against tol_rel, and inf when the bracket is open upwards.

The arguments are validated strictly, because a limit search that returns a number nobody can tell apart from an answer is worse than one that raises. u_max and tol_rel must be finite and strictly positive — NaN passes a plain <= 0 check and then makes every later comparison False — c_bounds must be finite, and max_iter must be an int of at least 1.

For the rho-f variant — useful when the bus impedance carries a frequency-coupling term whose magnitude should also be probed — see :func:find_max_rho_f_scaling. It reports the same four diagnostic keys with the same meaning.

Conductor thermal limits (IEC 60949 / IEC 60909-0)

The inverse-rho routines answer a person-safety question (does the EPR stay below a touch-voltage limit). check_conductor_limits answers the complementary equipment-integrity question: does the grounding conductor survive the fault thermally? It compares the thermally equivalent short-time current

\[ I_{th} = I_{s,\text{RMS}} \sqrt{m + n} \]

against the adiabatic limit of IEC 60949,

\[ I_{adm} = \frac{k \, S}{\sqrt{t_k}}, \qquad k = K \sqrt{\ln\!\frac{\theta_f + \beta}{\theta_i + \beta}} , \]

per grounding branch, and reports the utilisation and a pass/fail flag. A branch is checked only when its BranchType carries both conductor_material and cross_section_mm2.

What is superposed and what is not

This is the modelling rule the whole short-circuit side is built on. The frequency-domain solve superposes the linear AC-RMS currents as always. The IEC 60909 factors \(\kappa\) and \(m\) are non-linear in the fault-loop \(R/X\), so they are applied once, to the already aggregated branch current — \(i_p\) and \(I_{th}\) are never superposed directly.

With several infeeds the DC components add, so the largest possible peak of the total current is the sum of the individual peaks. Written as a single factor on the aggregate that is exactly the current-weighted mean

\[ \kappa_{\text{eff}} = \frac{\sum_i \kappa_i I_i}{\sum_i I_i} , \]

which resolve_fault_sc_characteristics uses by default (aggregation="weighted"); it reproduces the sum of the individual peaks identically. aggregation="max" is the strictly conservative variant. Reusing one source's \(\kappa\) for all of them is simply wrong and errs in either direction. Where the simultaneous-peak assumption is too crude — strongly mixed \(R/X\) infeeds — the transient solver remains the exact fallback, because it integrates the actual waveforms instead of applying a standard factor.

import groundinsight as gi

# The characteristics can be set by hand ...
net.sources["src"].r_to_x = 0.1
net.faults["F"].t_k_s = 0.5

# ... or imported from a solved pandapower case, see the I/O page:
# gi.apply_shortcircuit_characteristics(net, net_pp, "F")

gi.run_fault(net, "F")
gi.check_conductor_limits(net, "F").select(
    "branch_name", "I_s_rms_A", "i_p_A", "I_th_A",
    "I_admissible_A", "utilization", "within_limit",
)

Explicit t_k, kappa / r_to_x and n arguments override whatever is stored on the model, so sensitivity studies stay possible.

Node thermal limits — earthing conductor vs earth electrode

check_conductor_limits assesses the shield / earth wire between buses. check_node_limits assesses the two grounding elements at a bus, which EN 50522 / IEC 61936-1 keep strictly apart because they carry different currents. Confusing them is the classic sizing error, and in a meshed system it is an order-of-magnitude error in both directions.

element German current column
earthing conductor Erdungsleiter full injected earth-fault current ResultBus.i_inj
earth electrode Erder share dissipated into the soil, \(u_\text{EPR}/Z_B\) ResultBus.ia

Three physically distinct currents meet at a grounding bus: the lumped injection i_inj, the electrode current ia, and the branch shield currents reported per branch on ResultBranch. The nodal balance that ties them together is

\[ i_a = i_\text{vector} + \sum_\text{branches} (u_\text{other} - u_\text{self})\, Y_\text{self} . \]

i_inj deliberately excludes the mutual Norton-equivalent injections of the inductively coupled branches: those model a distributed induced EMF along the line, not a current entering the node through a lumped conductor. It is non-zero only at source buses (the infeed) and at the fault bus (the total fault current).

Data model and current split

A bus is assessed per element, and only once its BusType carries both the material and the cross-section for that element. The two elements are independent — declaring one does not imply the other.

bt = gi.BusType(
    name="tower",
    system_type="Tower",
    voltage_level=110.0,
    impedance_formula="rho * 0.1",
    # Erdungsleiter — sized for the full earth-fault current
    earthing_conductor_material="Cu",
    earthing_conductor_cross_section_mm2=50.0,
    earthing_conductor_theta_final_C=gi.final_temperature("Cu", "bare"),
    # Erder — sized for what actually reaches the soil, four legs
    electrode_material="Steel",
    electrode_cross_section_mm2=95.0,
    electrode_current_split=0.25,
)

Each element carries a free factor current_split in \((0, 1]\), applied as \(I_\text{conductor} = I_\text{RMS}\cdot\texttt{current\_split}\). It expresses how the bus current divides between physically parallel paths that the nodal model lumps into one: 1.0 (default) for a single conductor, 1/N for \(N\) parallel legs, 0.5 for a ring fed at one point, or an IEEE Std 80 division factor. It is deliberately not derived automatically — the split depends on geometry the nodal model does not carry — and a value above 1 is rejected at the model level, because that would not be a split but an error.

gi.run_fault(net, "F")
gi.check_node_limits(net, "F", t_k=0.5, r_to_x=0.1).select(
    "bus_name", "element", "I_rms_A", "current_split",
    "I_th_A", "I_admissible_A", "utilization", "within_limit",
)

The IEC 60909 excitation (t_k, kappa, m, n) is resolved by the same helper the branch check uses, so both views of one fault always agree. Buses whose type declares neither element still appear in the frame, with the currents filled in and within_limit = None — enough to size by hand without re-running anything.

Final temperatures

FINAL_TEMPERATURES and final_temperature(material, covering) provide \(\theta_f\) values with the source named inline per entry (National Grid ETS Table 5a for bare buried conductors, IEC 60364-5-54 Table 54.2 for PVC / XLPE). The catalog is deliberately incomplete rather than filled with plausible-looking numbers; final_temperature raises for anything missing and names EN 50522 Table 2 as the source to consult.

Steel default

The IEC60949_MATERIALS["Steel"] default of 400 °C is higher — i.e. more permissive, the unsafe direction for a limit check — than the 300 °C the National Grid table gives for bare buried steel. It is left unchanged because moving a default silently would move every existing study. Pass theta_final_C explicitly until the value has been checked against EN 50522 Table 2.

Results computed before this feature

ResultBus.i_inj defaults to 0.0 on results stored by an earlier version, so the earthing-conductor rows would read as unstressed. Re-run run_fault after upgrading.

Incomplete results are an error, not a pass

Both checks build their frame from the stored result of the fault, so a missing entry would produce no row — and a missing row is indistinguishable from a passing one. They therefore verify first that the stored result covers every branch and every active bus, and raise ValueError otherwise rather than reporting the missing elements as free of violations.

This matters because the gap is reachable through ordinary use. ElectricalNetwork.solve_network() replaces network.results[fault] with a result carrying bus rows only — compute_branch_currents() fills in the branches, which is why run_fault calls both — so solving a hand-built ElectricalNetwork to look at \(Y\), \(i\) or \(u\) discards the branch results. Adding a branch after run_fault leaves the same gap. To inspect the nodal system without disturbing anything, build the ElectricalNetwork (its constructor has no side effects) and read the voltages back from network.results[fault].buses instead of re-solving.

Both frames also carry an explicit schema, so a network genuinely without branches returns an empty frame that is still selectable and filterable, rather than a schema-less (0, 0) frame on which pl.col(...) raises.

API reference

analysis

Analysis subpackage.

Contains higher-level analysis routines on top of :mod:network_operations, such as inverse problems for the bus-grounding rho-f characteristic. These functions assume that the network is already fully built (buses, branches, faults, sources, paths) and orchestrate repeated :func:run_fault calls.

BusResponse

Bases: BaseModel

Closed-form response of the network to the electrode at one bus.

Built by :func:bus_response. Evaluating it costs no solve.

Attributes:

Name Type Description
fault str

Fault the response was built for.
bus str

Bus whose electrode is the free parameter.
fault_bus str

Bus the fault sits on. Used to anchor which group of buses is named as feeding the soil when the earth-return current is split.
frequencies list of float

Frequencies covered.
z_network dict of float to complex

Driving-point impedance at bus with the local electrode removed -- the parallel impedance the rest of the network offers at that location. Independent of anything installed at bus.
u_open dict of float to dict of str to complex

Nodal voltages with the local electrode removed.
z_column dict of float to dict of str to complex

Voltage at every bus per ampere injected at bus, source-free.
r_epr dict of float to float or None

The EPR-based reduction factor, carried through because it is constant along the whole curve.
i_fault dict of float to complex

Fault current, the negated sum of source injections. Constant.

admittance 

admittance(z_bus: ElectrodeSpec, freq: float) -> complex

Translate an electrode spelling into a shunt admittance.

None and an infinite impedance both mean no electrode (Y_B = 0); a zero impedance means an ideal one and is returned as inf, which :meth:voltages handles as a limit rather than a division.

Source code in src/groundinsight/analysis/response.py 
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def admittance(self, z_bus: ElectrodeSpec, freq: float) -> complex:
    """
    Translate an electrode spelling into a shunt admittance.

    ``None`` and an infinite impedance both mean *no electrode*
    (``Y_B = 0``); a zero impedance means an *ideal* one and is returned as
    ``inf``, which :meth:`voltages` handles as a limit rather than a
    division.
    """
    if z_bus is None:
        return 0.0 + 0.0j
    z = complex(z_bus)
    if np.isinf(z.real) or np.isinf(z.imag):
        return 0.0 + 0.0j
    if z == 0:
        return complex(np.inf, 0.0)
    return 1.0 / z

driving_point 

driving_point(z_bus: ElectrodeSpec, freq: float) -> complex

Z_dp = 1 / (Y_B + 1/Z_net) -- the impedance seen at bus.

Source code in src/groundinsight/analysis/response.py 
199
200
201
202
203
204
205
206
207
208
209
210
def driving_point(self, z_bus: ElectrodeSpec, freq: float) -> complex:
    """``Z_dp = 1 / (Y_B + 1/Z_net)`` -- the impedance seen at ``bus``."""
    y = self.admittance(z_bus, freq)
    if np.isinf(y.real):
        return 0.0 + 0.0j
    z_net = self.z_network[freq]
    if z_net == 0:
        return 0.0 + 0.0j
    total = y + 1.0 / z_net
    if total == 0:
        return complex(np.inf, 0.0)
    return 1.0 / total

evaluate 

evaluate(
    z_bus: ElectrodeSpec, *, label: Optional[str] = None
) -> pl.DataFrame

Everything the location does for one electrode, as one frame.

Parameters:

Name Type Description Default
z_bus (complex, float or None)

The electrode. None or an infinite value means none installed, 0 means ideal. Both are evaluated as exact limits.

 required
label str

Value of the case column. Defaults to a readable rendering of z_bus.

 None

Returns:

Type Description
DataFrame

One row per frequency: case, bus, fault, frequency_Hz, Z_bus_Ohm, Z_network_Ohm, Z_driving_point_Ohm, EPR_V, EPR_deg, I_electrode_A, r_epr, r_current, Z_G_Ohm, plus EPR_<bus>_V for every other bus.
Source code in src/groundinsight/analysis/response.py 
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
def evaluate(
    self, z_bus: ElectrodeSpec, *, label: Optional[str] = None
) -> pl.DataFrame:
    """
    Everything the location does for one electrode, as one frame.

    Parameters
    ----------
    z_bus : complex, float or None
        The electrode. ``None`` or an infinite value means *none installed*,
        ``0`` means *ideal*. Both are evaluated as exact limits.
    label : str, optional
        Value of the ``case`` column. Defaults to a readable rendering of
        ``z_bus``.

    Returns
    -------
    pl.DataFrame
        One row per frequency: ``case``, ``bus``, ``fault``,
        ``frequency_Hz``, ``Z_bus_Ohm``, ``Z_network_Ohm``,
        ``Z_driving_point_Ohm``, ``EPR_V``, ``EPR_deg``, ``I_electrode_A``,
        ``r_epr``, ``r_current``, ``Z_G_Ohm``, plus ``EPR_<bus>_V`` for
        every other bus.
    """
    rows = []
    for freq in self.frequencies:
        u = self.voltages(z_bus, freq)
        z_dp = self.driving_point(z_bus, freq)
        u_b = u[self.bus]
        r_epr = self.r_epr.get(freq)
        i_fault = self.i_fault.get(freq)
        z_g = (
            u_b / (r_epr * i_fault)
            if r_epr not in (None, 0) and i_fault not in (None, 0)
            else None
        )
        row = {
            "case": label if label is not None else _render(z_bus),
            "fault": self.fault,
            "bus": self.bus,
            "frequency_Hz": float(freq),
            "Z_bus_Ohm": _magnitude(z_bus),
            "Z_network_Ohm": float(abs(self.z_network[freq])),
            "Z_driving_point_Ohm": float(abs(z_dp)),
            "EPR_V": float(abs(u_b)),
            "EPR_deg": float(np.degrees(np.angle(u_b))),
            "I_electrode_A": float(
                abs(self._electrode_current(z_bus, freq, u))
            ),
            "r_epr": r_epr,
            "r_current": self._earth_return_factor(z_bus, freq, u),
            "Z_G_Ohm": None if z_g is None else float(abs(z_g)),
        }
        for name in self.bus_names:
            if name != self.bus:
                row[f"EPR_{name}_V"] = float(abs(u[name]))
        rows.append(row)
    # Columns that are legitimately all-null for one case ("open" has no
    # finite electrode impedance) would otherwise come out as Null dtype and
    # refuse to stack with the other cases.
    nullable = ["Z_bus_Ohm", "r_epr", "r_current", "Z_G_Ohm"]
    return pl.DataFrame(rows).with_columns(
        [pl.col(name).cast(pl.Float64) for name in nullable]
    )

extremes 

extremes() -> pl.DataFrame

The bracket: no electrode, ideal electrode, and the passive worst case.

Three rows per frequency, labelled "open", "ideal" and "worst_passive". The first two are the endpoints of the curve; the third is the reactive electrode that maximises the driving-point magnitude (see :meth:worst_case_electrode).

Source code in src/groundinsight/analysis/response.py 
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
def extremes(self) -> pl.DataFrame:
    """
    The bracket: no electrode, ideal electrode, and the passive worst case.

    Three rows per frequency, labelled ``"open"``, ``"ideal"`` and
    ``"worst_passive"``. The first two are the endpoints of the curve; the
    third is the reactive electrode that maximises the driving-point
    magnitude (see :meth:`worst_case_electrode`).
    """
    frames = [
        self.evaluate(None, label="open"),
        self.evaluate(0.0, label="ideal"),
    ]
    worst_rows = []
    for freq in self.frequencies:
        worst = self.worst_case_electrode(freq)
        frame = self.evaluate(worst["z_bus"], label="worst_passive")
        worst_rows.append(frame.filter(pl.col("frequency_Hz") == float(freq)))
    frames.append(pl.concat(worst_rows, how="diagonal"))
    return pl.concat(frames, how="diagonal")

sweep 

sweep(
    z_values: Sequence[ElectrodeSpec],
    *,
    labels: Optional[Sequence[str]] = None
) -> pl.DataFrame

Evaluate many electrodes at once. No solve, whatever the length.

Parameters:

Name Type Description Default
z_values sequence

Electrodes, in the spellings :meth:evaluate accepts.

 required
labels sequence of str

One label per value. Defaults to a rendering of each value.

 None
Source code in src/groundinsight/analysis/response.py 
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
def sweep(
    self,
    z_values: Sequence[ElectrodeSpec],
    *,
    labels: Optional[Sequence[str]] = None,
) -> pl.DataFrame:
    """
    Evaluate many electrodes at once. No solve, whatever the length.

    Parameters
    ----------
    z_values : sequence
        Electrodes, in the spellings :meth:`evaluate` accepts.
    labels : sequence of str, optional
        One label per value. Defaults to a rendering of each value.
    """
    if labels is not None and len(labels) != len(z_values):
        raise ValueError(
            f"{len(z_values)} electrode(s) were given but {len(labels)} "
            f"label(s); they have to correspond one to one."
        )
    frames = [
        self.evaluate(z, label=None if labels is None else labels[i])
        for i, z in enumerate(z_values)
    ]
    return pl.concat(frames, how="diagonal")

voltages 

voltages(
    z_bus: ElectrodeSpec, freq: float
) -> Dict[str, complex]

Nodal voltages for one electrode at one frequency.

Both endpoints are evaluated as limits, so an ideal electrode gives exactly zero at bus instead of a very small number.

Source code in src/groundinsight/analysis/response.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
def voltages(self, z_bus: ElectrodeSpec, freq: float) -> Dict[str, complex]:
    """
    Nodal voltages for one electrode at one frequency.

    Both endpoints are evaluated as limits, so an ideal electrode gives
    exactly zero at ``bus`` instead of a very small number.
    """
    u0 = self.u_open[freq]
    z_col = self.z_column[freq]
    u0_b = u0[self.bus]
    z_net = self.z_network[freq]
    y = self.admittance(z_bus, freq)

    if np.isinf(y.real):
        # Ideal electrode: the Y_B -> inf limit of the Möbius form.
        if z_net == 0:
            return dict(u0)
        return {
            name: u0[name] - z_col[name] * u0_b / z_net
            for name in self.bus_names
        }
    if y == 0:
        return dict(u0)
    denominator = 1.0 + y * z_net
    if denominator == 0:
        # A reactive electrode exactly at the network's own pole. Physically
        # unreachable with a passive electrode (it needs Re(Y_B) < 0), but
        # named rather than returned as inf.
        raise ValueError(
            f"The electrode {z_bus!r} at bus '{self.bus}' cancels the "
            f"network admittance exactly at {freq} Hz (1 + Y_B*Z_net = 0), "
            f"so the response has a pole there. That requires a negative "
            f"conductance, which no passive electrode has."
        )
    factor = y * u0_b / denominator
    return {
        name: u0[name] - z_col[name] * factor for name in self.bus_names
    }

worst_case_electrode 

worst_case_electrode(freq: float) -> Dict[str, complex]

The passive electrode that maximises |Z_dp|, and the value it gives.

Over the closed right half-plane of Y_B the magnitude of 1/(Y_B + Y_net) is largest where the imaginary parts cancel and the real part is as small as it can be, i.e. at Y_B = -j*Im(Y_net). The result exceeds |Z_net| only slightly in a cable network, but it is the true bound rather than an assumed one.

Returns:

Type Description
dict

z_bus (the electrode, purely reactive) and z_driving_point.
Source code in src/groundinsight/analysis/response.py 
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
def worst_case_electrode(self, freq: float) -> Dict[str, complex]:
    """
    The passive electrode that maximises ``|Z_dp|``, and the value it gives.

    Over the closed right half-plane of ``Y_B`` the magnitude of
    ``1/(Y_B + Y_net)`` is largest where the imaginary parts cancel and the
    real part is as small as it can be, i.e. at ``Y_B = -j*Im(Y_net)``. The
    result exceeds ``|Z_net|`` only slightly in a cable network, but it is
    the true bound rather than an assumed one.

    Returns
    -------
    dict
        ``z_bus`` (the electrode, purely reactive) and ``z_driving_point``.
    """
    z_net = self.z_network[freq]
    if z_net == 0 or not np.isfinite(z_net):
        return {"z_bus": complex(np.inf, 0.0), "z_driving_point": z_net}
    y_net = 1.0 / z_net
    y_b = complex(0.0, -y_net.imag)
    if y_b == 0:
        return {"z_bus": complex(np.inf, 0.0), "z_driving_point": z_net}
    return {"z_bus": 1.0 / y_b, "z_driving_point": 1.0 / y_net.real}

Cut

Bases: BaseModel

A named set of branches separating one side of the network from the fault.

Attributes:

Name Type Description
name str

Label used in the result frame. Must not collide with the reserved name "rest", which is used for the implicit remainder side.
branches list of str

Names of the branches forming the cut. All of them must exist, be active, carry a grounding conductor and be incident to the fault bus.
description (str, optional)

Free text carried through to the result frame.

CutAnalysis

Bases: BaseModel

Result of :func:analyze_cuts for one fault.

Attributes:

Name Type Description
fault str

Name of the analysed fault.
fault_bus str

Bus the fault sits on -- the reference point of the decomposition.
frequencies list of float

Frequencies the analysis was run at.
sides dict of str to list of str

Buses visible from the fault bus when looking out through that cut, with the other cuts' branches removed. In a ring the entries overlap.
sides_are_disjoint bool

True when no two directions reach the same bus. The far-side current quantities are only defined in that case.
branches dict of str to list of str

The branches of each cut, including the implicit "rest" side.
z_local dict of float to complex

Impedance of the fault bus's own electrode per frequency.
z_side dict of str to dict of float to complex

Impedance of each direction, from source-free current division.
z_parallel dict of float to complex

z_local in parallel with every z_side.
z_driving_point dict of float to complex

Driving-point impedance of the complete network at the fault bus, computed independently. Equal to z_parallel up to identity_residual.
identity_residual dict of float to float

|z_parallel - z_driving_point| / |z_driving_point| per frequency.
i_shield dict of str to dict of float to complex

Shield current leaving the fault bus in each direction. Empty when no solved result was available.
i_fault, i_local dict of float to complex

Injection at the fault bus and the part of it taken by the fault bus's own electrode.
current_share dict of str to dict of float to float or None

|i_shield| / |i_fault| -- the share of the fault current leaving metallically in that direction. Always defined.
i_earth, i_total dict of str to dict of float to complex or None

Soil and total current crossing the cut, summed over the far side. None per entry when the directions overlap.
r_side dict of str to dict of float to float or None

|i_earth| / |i_total| per direction -- the earth-return share of what crosses that cut. None when the far side carries no injection or when the directions overlap.
kcl_residual dict of float to float

|i_fault - i_local - sum(i_shield)| / |i_fault|. Zero to machine precision confirms the split accounts for the whole fault current.

has_currents property 

has_currents: bool

True when a solved fault result was available.

to_polars 

to_polars() -> pl.DataFrame

Return the analysis as one long-format frame, one row per direction and frequency.

Columns

fault, fault_bus, cut, side_buses, n_side_buses, branches, sides_are_disjoint, frequency_Hz, Z_side_Ohm, Z_side_deg, Z_local_Ohm, Z_parallel_Ohm, Z_driving_point_Ohm, identity_residual, I_fault_A, I_local_A, I_shield_A, I_shield_deg, current_share, I_earth_A, I_total_A, r_side, kcl_residual

Returns:

Type Description
DataFrame

One row per (cut, frequency).
Source code in src/groundinsight/analysis/decomposition.py 
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
def to_polars(self) -> pl.DataFrame:
    """
    Return the analysis as one long-format frame, one row per direction and
    frequency.

    Columns
    -------
    fault, fault_bus, cut, side_buses, n_side_buses, branches,
    sides_are_disjoint, frequency_Hz, Z_side_Ohm, Z_side_deg, Z_local_Ohm,
    Z_parallel_Ohm, Z_driving_point_Ohm, identity_residual, I_fault_A,
    I_local_A, I_shield_A, I_shield_deg, current_share, I_earth_A,
    I_total_A, r_side, kcl_residual

    Returns
    -------
    pl.DataFrame
        One row per (cut, frequency).
    """

    def _mag(value):
        return None if value is None else float(abs(value))

    rows = []
    for cut_name in self.z_side:
        for freq in self.frequencies:
            z = self.z_side[cut_name][freq]
            row = {
                "fault": self.fault,
                "fault_bus": self.fault_bus,
                "cut": cut_name,
                "side_buses": ", ".join(self.sides[cut_name]),
                "n_side_buses": len(self.sides[cut_name]),
                "branches": ", ".join(self.branches[cut_name]),
                "sides_are_disjoint": self.sides_are_disjoint,
                "frequency_Hz": float(freq),
                "Z_side_Ohm": float(abs(z)),
                "Z_side_deg": float(np.degrees(np.angle(z))),
                "Z_local_Ohm": float(abs(self.z_local[freq])),
                "Z_parallel_Ohm": float(abs(self.z_parallel[freq])),
                "Z_driving_point_Ohm": float(abs(self.z_driving_point[freq])),
                "identity_residual": float(self.identity_residual[freq]),
            }
            if self.has_currents:
                i_sh = self.i_shield[cut_name][freq]
                row.update(
                    {
                        "I_fault_A": _mag(self.i_fault[freq]),
                        "I_local_A": _mag(self.i_local[freq]),
                        "I_shield_A": float(abs(i_sh)),
                        "I_shield_deg": float(np.degrees(np.angle(i_sh))),
                        "current_share": self.current_share[cut_name][freq],
                        "I_earth_A": _mag(self.i_earth[cut_name][freq]),
                        "I_total_A": _mag(self.i_total[cut_name][freq]),
                        "r_side": self.r_side[cut_name][freq],
                        "kcl_residual": float(self.kcl_residual[freq]),
                    }
                )
            else:
                row.update(
                    {
                        "I_fault_A": None,
                        "I_local_A": None,
                        "I_shield_A": None,
                        "I_shield_deg": None,
                        "current_share": None,
                        "I_earth_A": None,
                        "I_total_A": None,
                        "r_side": None,
                        "kcl_residual": None,
                    }
                )
            rows.append(row)
    return pl.DataFrame(rows)

FaultShortCircuitData

Bases: BaseModel

Resolved IEC 60909 characteristics of one fault.

The result of aggregating the per-source data (:attr:Source.i_k_a <groundinsight.models.core_models.Source.i_k_a>, r_to_x, kappa) with the fault-level data (:attr:Fault.t_k_s <groundinsight.models.core_models.Fault.t_k_s>, n_factor) into the single set of numbers the non-linear 60909 factors are evaluated with.

Attributes:

Name Type Description
fault_name str

Name of the fault these characteristics belong to.
frequency float

Frequency in Hz at which the source currents were weighted.
kappa (float, optional)

Effective peak factor, or None when no participating source carries kappa / r_to_x.
r_to_x (float, optional)

Effective R/X, reported only when it is unambiguous (all contributing sources share one value).
t_k_s (float, optional)

Fault duration T_k in seconds, taken from the fault.
n_factor float

AC-decay heat factor n taken from the fault.
m (float, optional)

DC heat-effect factor at kappa, frequency and t_k_s. None when either kappa or t_k_s is unknown.
i_k_a float

Arithmetic sum of the participating source injection magnitudes at frequency (in amperes), i.e. the weight base of the aggregation. This is not the branch current -- that comes from the solve.
i_p_a (float, optional)

kappa * sqrt(2) * i_k_a, reported for reference.
aggregation str

"weighted" or "max".
homogeneous bool

True when every contributing source shares the same kappa (within 1e-9), so the aggregation is exact rather than an interpolation.
sources list of str

Names of the sources that contributed a kappa.
sources_without_kappa list of str

Names of sources that inject current at frequency but carry no 60909 data; they are part of the linear solve but excluded from the kappa aggregation.

ReferenceCase 

ReferenceCase(
    name: str,
    quantity: str,
    conditions: str,
    tolerance: float,
    run: Callable[[], Tuple[float, float]],
)

One closed-form case and the model run that has to match it.

Attributes:

Name Type Description
name str

Short identifier, used as the row key.
quantity str

What is being compared, with its unit.
conditions str

The boundary conditions under which the closed form holds. A deviation outside tolerance means either the model is wrong or a condition here was not met -- the second is the more common finding.
tolerance float

Relative deviation still counted as agreement.
Source code in src/groundinsight/analysis/reference.py 
87
88
89
90
91
92
93
94
95
96
97
98
99
def __init__(
    self,
    name: str,
    quantity: str,
    conditions: str,
    tolerance: float,
    run: Callable[[], Tuple[float, float]],
):
    self.name = name
    self.quantity = quantity
    self.conditions = conditions
    self.tolerance = tolerance
    self._run = run

evaluate 

evaluate() -> Dict[str, object]

Run the case and return its row.

Source code in src/groundinsight/analysis/reference.py 
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def evaluate(self) -> Dict[str, object]:
    """Run the case and return its row."""
    closed_form, model = self._run()
    reference = abs(closed_form)
    deviation = (
        abs(model - closed_form) / reference if reference > 0 else float("nan")
    )
    return {
        "case": self.name,
        "quantity": self.quantity,
        "conditions": self.conditions,
        "closed_form": float(closed_form),
        "model": float(model),
        "rel_deviation": float(deviation),
        "tolerance": float(self.tolerance),
        "agrees": bool(deviation <= self.tolerance),
    }

admissible_short_circuit_current 

admissible_short_circuit_current(
    k: float, cross_section_mm2: float, t_k: float
) -> float

Adiabatic admissible short-circuit current I_adm = k*S/sqrt(t_k).

Parameters:

Name Type Description Default
k float

Material constant from :func:iec60949_k, in A·s^0.5/mm².

 required
cross_section_mm2 float

Conductor cross-section S in mm². Must be strictly positive.

 required
t_k float

Fault duration in seconds. Must be strictly positive.

 required

Returns:

Type Description
float

The admissible short-circuit current in amperes.

Raises:

Type Description
ValueError

If cross_section_mm2 or t_k is not strictly positive.
Source code in src/groundinsight/analysis/thermal.py 
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
def admissible_short_circuit_current(
    k: float, cross_section_mm2: float, t_k: float
) -> float:
    """
    Adiabatic admissible short-circuit current ``I_adm = k*S/sqrt(t_k)``.

    Parameters
    ----------
    k : float
        Material constant from :func:`iec60949_k`, in A·s^0.5/mm².
    cross_section_mm2 : float
        Conductor cross-section ``S`` in mm². Must be strictly positive.
    t_k : float
        Fault duration in seconds. Must be strictly positive.

    Returns
    -------
    float
        The admissible short-circuit current in amperes.

    Raises
    ------
    ValueError
        If ``cross_section_mm2`` or ``t_k`` is not strictly positive.
    """
    if cross_section_mm2 <= 0:
        raise ValueError(
            f"cross_section_mm2 must be strictly positive, got {cross_section_mm2!r}."
        )
    if t_k <= 0:
        raise ValueError(f"t_k must be strictly positive, got {t_k!r}.")
    return k * cross_section_mm2 / math.sqrt(t_k)

analyze_cuts 

analyze_cuts(
    network: Network,
    *,
    fault: str,
    cuts: Sequence[Cut],
    include_currents: bool = True
) -> CutAnalysis

Split the network at the fault bus and quantify each direction.

Parameters:

Name Type Description Default
network Network

The network. Bus and branch impedances must have been evaluated at network.frequencies (they are, after add_bus / add_branch).

 required
fault str

Name of the fault whose bus is the reference of the decomposition.

 required
cuts sequence of Cut

Named sets of branches, each incident to the fault bus and disjoint from the others. Incident branches that no cut claims form the implicit side "rest".

 required
include_currents bool

If True (default) the current split is read from network.results[fault]. When no such result exists the impedances are still returned and the current fields stay empty -- the impedance half of the analysis needs no solve.

 True

Returns:

Type Description
CutAnalysis

Impedances, currents and the residual of the parallel identity.

Raises:

Type Description
ValueError

If the fault is unknown, if a cut names a branch that is not an active grounding branch at the fault bus, or if two cuts claim the same branch.

Examples:

>>> import groundinsight as gi
>>> gi.run_fault(net, "F1")
>>> analysis = gi.analyze_cuts(
...     net,
...     fault="F1",
...     cuts=[gi.Cut(name="left", branches=["L12"]),
...           gi.Cut(name="right", branches=["L23"])],
... )
>>> analysis.to_polars()
Source code in src/groundinsight/analysis/decomposition.py 
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
def analyze_cuts(
    network: Network,
    *,
    fault: str,
    cuts: Sequence[Cut],
    include_currents: bool = True,
) -> CutAnalysis:
    """
    Split the network at the fault bus and quantify each direction.

    Parameters
    ----------
    network : Network
        The network. Bus and branch impedances must have been evaluated at
        ``network.frequencies`` (they are, after ``add_bus`` / ``add_branch``).
    fault : str
        Name of the fault whose bus is the reference of the decomposition.
    cuts : sequence of Cut
        Named sets of branches, each incident to the fault bus and disjoint from
        the others. Incident branches that no cut claims form the implicit side
        ``"rest"``.
    include_currents : bool, optional
        If ``True`` (default) the current split is read from
        ``network.results[fault]``. When no such result exists the impedances
        are still returned and the current fields stay empty -- the impedance
        half of the analysis needs no solve.

    Returns
    -------
    CutAnalysis
        Impedances, currents and the residual of the parallel identity.

    Raises
    ------
    ValueError
        If the fault is unknown, if a cut names a branch that is not an active
        grounding branch at the fault bus, or if two cuts claim the same branch.

    Examples
    --------
    >>> import groundinsight as gi  # doctest: +SKIP
    >>> gi.run_fault(net, "F1")  # doctest: +SKIP
    >>> analysis = gi.analyze_cuts(  # doctest: +SKIP
    ...     net,
    ...     fault="F1",
    ...     cuts=[gi.Cut(name="left", branches=["L12"]),
    ...           gi.Cut(name="right", branches=["L23"])],
    ... )
    >>> analysis.to_polars()  # doctest: +SKIP
    """
    if fault not in network.faults:
        raise ValueError(
            f"Fault '{fault}' does not exist in network '{network.name}'. "
            f"Available: {sorted(network.faults)}."
        )
    fault_bus = network.faults[fault].bus
    resolved, incident = _validate_cuts(network, cuts, fault_bus)
    frequencies = [float(f) for f in network.frequencies]

    directions = {
        name: _direction_buses(network, set(branches), incident, fault_bus)
        for name, branches in resolved.items()
    }
    names = list(resolved)
    disjoint = all(
        not (set(directions[a]) & set(directions[b]))
        for i, a in enumerate(names)
        for b in names[i + 1 :]
    )

    z_local: Dict[float, complex] = {}
    z_side: Dict[str, Dict[float, complex]] = {name: {} for name in resolved}
    z_parallel: Dict[float, complex] = {}
    z_driving_point: Dict[float, complex] = {}
    identity_residual: Dict[float, float] = {}

    infinite = complex(np.inf, 0.0)
    for freq in frequencies:
        local = _finite_impedance(network.buses[fault_bus].impedance.get(freq))
        z_local[freq] = local if local is not None else infinite

        u = _unit_injection(network, fault_bus, freq)
        if u is None:
            z_driving_point[freq] = infinite
            for name in resolved:
                z_side[name][freq] = infinite
            z_parallel[freq] = infinite
            identity_residual[freq] = float("nan")
            continue

        u_fault = u[fault_bus]
        z_driving_point[freq] = u_fault

        admittance = 0.0 + 0.0j if local is None else 1.0 / local
        for name, branches in resolved.items():
            i_cut = sum(
                (
                    _branch_current_out(network, b, fault_bus, u, freq)
                    for b in branches
                ),
                0.0 + 0.0j,
            )
            if i_cut == 0 or u_fault == 0:
                z_side[name][freq] = infinite
            else:
                z = u_fault / i_cut
                z_side[name][freq] = z
                admittance += 1.0 / z

        z_parallel[freq] = 1.0 / admittance if admittance != 0 else infinite
        reference = abs(z_driving_point[freq])
        identity_residual[freq] = (
            abs(z_parallel[freq] - z_driving_point[freq]) / reference
            if reference > 0 and np.isfinite(reference)
            else float("nan")
        )

    analysis_kwargs = dict(
        fault=fault,
        fault_bus=fault_bus,
        frequencies=frequencies,
        sides=directions,
        sides_are_disjoint=disjoint,
        branches=resolved,
        z_local=z_local,
        z_side=z_side,
        z_parallel=z_parallel,
        z_driving_point=z_driving_point,
        identity_residual=identity_residual,
    )

    if include_currents and fault in network.results:
        result = network.results[fault]
        branch_results = {b.name: b for b in result.branches}
        bus_results = {b.name: b for b in result.buses}

        i_shield: Dict[str, Dict[float, complex]] = {n: {} for n in resolved}
        i_earth: Dict[str, Dict[float, complex]] = {n: {} for n in resolved}
        i_total: Dict[str, Dict[float, complex]] = {n: {} for n in resolved}
        r_side: Dict[str, Dict[float, Optional[float]]] = {n: {} for n in resolved}
        current_share: Dict[str, Dict[float, Optional[float]]] = {
            n: {} for n in resolved
        }
        i_fault: Dict[float, complex] = {}
        i_local: Dict[float, complex] = {}
        kcl_residual: Dict[float, float] = {}

        for freq in frequencies:
            injection = complex(bus_results[fault_bus].i_inj_freq[freq])
            local_current = complex(bus_results[fault_bus].ia_freq[freq])
            i_fault[freq] = injection
            i_local[freq] = local_current
            leaving = 0.0 + 0.0j

            for name, branches in resolved.items():
                total_out = 0.0 + 0.0j
                for branch_name in branches:
                    branch = network.branches[branch_name]
                    stored = complex(branch_results[branch_name].i_s_freq[freq])
                    # ``i_s`` is oriented to_bus -> from_bus, so it already is
                    # the current leaving the fault bus when the fault bus is
                    # the to_bus, and the negative of it otherwise.
                    sign = 1.0 if branch.to_bus == fault_bus else -1.0
                    total_out += sign * stored
                i_shield[name][freq] = total_out
                leaving += total_out
                current_share[name][freq] = (
                    abs(total_out) / abs(injection) if injection != 0 else None
                )

                if disjoint:
                    far = directions[name]
                    earth = sum(
                        (complex(bus_results[b].ia_freq[freq]) for b in far),
                        0.0 + 0.0j,
                    )
                    total = sum(
                        (complex(bus_results[b].i_inj_freq[freq]) for b in far),
                        0.0 + 0.0j,
                    )
                    i_earth[name][freq] = earth
                    i_total[name][freq] = total
                    r_side[name][freq] = (
                        abs(earth) / abs(total) if abs(total) > 0 else None
                    )
                else:
                    i_earth[name][freq] = None
                    i_total[name][freq] = None
                    r_side[name][freq] = None

            reference = abs(injection)
            kcl_residual[freq] = (
                abs(injection - local_current - leaving) / reference
                if reference > 0
                else 0.0
            )

        if not disjoint:
            logger.info(
                "The directions out of fault bus '%s' overlap -- %s reach each "
                "other around a ring -- so there is no far side belonging to "
                "one cut alone. The impedances are unaffected (they come from "
                "current division), but I_earth, I_total and r_side are left "
                "empty; use current_share instead.",
                fault_bus,
                " and ".join(f"'{n}'" for n in names),
            )

        analysis_kwargs.update(
            i_shield=i_shield,
            i_earth=i_earth,
            i_total=i_total,
            r_side=r_side,
            current_share=current_share,
            i_fault=i_fault,
            i_local=i_local,
            kcl_residual=kcl_residual,
        )
    elif include_currents:
        logger.info(
            "No solved result for fault '%s' on network '%s', so only the "
            "source-free impedances are reported. Call run_fault first if you "
            "also want the current split.",
            fault,
            network.name,
        )

    return CutAnalysis(**analysis_kwargs)

bus_response 

bus_response(
    network: Network,
    *,
    fault: str,
    bus: Optional[str] = None
) -> BusResponse

Build the closed-form response of the network to the electrode at one bus.

Requires a solved fault: the assembled nodal system is reused, the bus's own shunt is taken back out of it, and two systems are solved once each -- the fault with the electrode removed, and a unit injection at the bus. Neither depends on the electrode, which is why the result covers every electrode.

Parameters:

Name Type Description Default
network Network

A network with network.results[fault] present, i.e. after :func:~groundinsight.network_operations.run_fault.

 required
fault str

The solved fault.

 required
bus str

Bus whose electrode is the free parameter. Defaults to the fault bus.

 None

Returns:

Type Description
BusResponse

Raises:

Type Description
ValueError

If the fault has not been solved, if the bus is unknown or inactive, or if the network without that electrode has no path to reference earth -- in which case there is nothing to characterise, because the location's behaviour is its own electrode.

Examples:

>>> gi.run_fault(net, "F1")
>>> response = gi.bus_response(net, fault="F1")
>>> response.extremes()
>>> response.evaluate(7.5)   # any electrode, no solve
Source code in src/groundinsight/analysis/response.py 
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
def bus_response(
    network: Network,
    *,
    fault: str,
    bus: Optional[str] = None,
) -> BusResponse:
    """
    Build the closed-form response of the network to the electrode at one bus.

    Requires a solved fault: the assembled nodal system is reused, the bus's own
    shunt is taken back out of it, and two systems are solved once each -- the
    fault with the electrode removed, and a unit injection at the bus. Neither
    depends on the electrode, which is why the result covers every electrode.

    Parameters
    ----------
    network : Network
        A network with ``network.results[fault]`` present, i.e. after
        :func:`~groundinsight.network_operations.run_fault`.
    fault : str
        The solved fault.
    bus : str, optional
        Bus whose electrode is the free parameter. Defaults to the fault bus.

    Returns
    -------
    BusResponse

    Raises
    ------
    ValueError
        If the fault has not been solved, if the bus is unknown or inactive, or
        if the network without that electrode has no path to reference earth --
        in which case there is nothing to characterise, because the location's
        behaviour *is* its own electrode.

    Examples
    --------
    >>> gi.run_fault(net, "F1")  # doctest: +SKIP
    >>> response = gi.bus_response(net, fault="F1")  # doctest: +SKIP
    >>> response.extremes()  # doctest: +SKIP
    >>> response.evaluate(7.5)   # any electrode, no solve  # doctest: +SKIP
    """
    if fault not in network.results:
        raise ValueError(
            f"Fault '{fault}' has no result on network '{network.name}'. The "
            f"response is read off the assembled nodal system, so call "
            f"run_fault first."
        )
    electrical = network.electrical_network
    if electrical is None or not electrical.Y_matrices:
        raise ValueError(
            f"Network '{network.name}' carries no assembled electrical network. "
            f"Call run_fault before asking for a bus response."
        )

    target = bus if bus is not None else network.faults[fault].bus
    if target not in electrical.bus_indices:
        known = sorted(electrical.bus_indices)
        raise ValueError(
            f"Bus '{target}' is not part of the solved system of network "
            f"'{network.name}'. Active buses: {known}."
        )

    index = electrical.bus_indices
    order = sorted(index, key=lambda name: index[name])
    k = index[target]
    frequencies = [float(f) for f in network.frequencies]

    result = network.results[fault]
    reduction = result.reduction_factor
    r_epr = dict(reduction.value) if reduction is not None else {}

    z_network: Dict[float, complex] = {}
    u_open: Dict[float, Dict[str, complex]] = {}
    z_column: Dict[float, Dict[str, complex]] = {}
    other_impedance: Dict[float, Dict[str, complex]] = {}
    i_fault: Dict[float, complex] = {}

    for freq in frequencies:
        Y = electrical.Y_matrices.get(freq)
        i_vector = electrical.i_vectors.get(freq)
        if Y is None or i_vector is None:
            continue
        Y0 = np.array(Y, dtype=complex, copy=True)

        # Take the bus's own shunt back out. The same rules as the assembly:
        # a missing or infinite impedance never contributed one in the first
        # place.
        z_bus_now = electrical._resolved_impedance(
            network.buses[target].impedance.get(freq), freq
        )
        if z_bus_now is not None and np.isfinite(z_bus_now) and z_bus_now != 0:
            Y0[k, k] -= 1.0 / z_bus_now

        # The check is structural, not a tolerance: at least one *other* bus has
        # to carry a finite, non-zero grounding impedance once the target's
        # shunt is gone. Leaving it to the linear solver is not enough --
        # measured on the exactly singular case, numpy returns a finite but
        # meaningless -2.25e15 instead of raising, the same inconsistency the
        # DC work found in scipy's splu.
        earthed_elsewhere = any(
            (
                z := electrical._resolved_impedance(
                    network.buses[name].impedance.get(freq), freq
                )
            )
            is not None
            and np.isfinite(z)
            and z != 0
            for name in order
            if name != target
        )
        if not earthed_elsewhere:
            raise ValueError(
                f"With the electrode at bus '{target}' removed, no other bus of "
                f"network '{network.name}' has a finite grounding impedance at "
                f"{freq} Hz, so there is no path to reference earth. There is "
                f"nothing to characterise independently of that electrode: at "
                f"this location the network's behaviour *is* the electrode. "
                f"Pick a bus that has at least one other earthed bus behind it."
            )

        unit = np.zeros(len(order), dtype=complex)
        unit[k] = 1.0
        try:
            u0 = np.linalg.solve(Y0, i_vector)
            z_col = np.linalg.solve(Y0, unit)
        except np.linalg.LinAlgError as exc:
            raise ValueError(
                f"The nodal system of network '{network.name}' is singular at "
                f"{freq} Hz once the electrode at bus '{target}' is removed, "
                f"even though other buses are earthed. Check for a subnetwork "
                f"that reaches earth only through '{target}'."
            ) from exc
        if not (np.all(np.isfinite(u0)) and np.all(np.isfinite(z_col))):
            raise ValueError(
                f"The response of bus '{target}' at {freq} Hz came out "
                f"non-finite. That points at an infinite or NaN impedance "
                f"somewhere in network '{network.name}' rather than at the bus "
                f"itself."
            )

        z_network[freq] = complex(z_col[k])
        u_open[freq] = {name: complex(u0[index[name]]) for name in order}
        z_column[freq] = {name: complex(z_col[index[name]]) for name in order}

        impedances: Dict[str, complex] = {}
        for name in order:
            z = electrical._resolved_impedance(
                network.buses[name].impedance.get(freq), freq
            )
            if z is not None and np.isfinite(z) and z != 0:
                impedances[name] = z
        other_impedance[freq] = impedances

        total_source = electrical.total_source_currents.get(freq)
        if total_source is not None:
            i_fault[freq] = -complex(total_source)

    if not z_network:
        raise ValueError(
            f"No frequency of fault '{fault}' carried an assembled system, so "
            f"there is no response to build."
        )

    return BusResponse(
        fault=fault,
        bus=target,
        fault_bus=network.faults[fault].bus,
        frequencies=[f for f in frequencies if f in z_network],
        bus_names=order,
        z_network=z_network,
        u_open=u_open,
        z_column=z_column,
        other_impedance=other_impedance,
        r_epr=r_epr,
        i_fault=i_fault,
    )

check_conductor_limits 

check_conductor_limits(
    network: Network,
    fault_name: str,
    t_k: Optional[float] = None,
    *,
    kappa: Optional[float] = None,
    r_to_x: Optional[float] = None,
    n: Optional[float] = None,
    f: Optional[float] = None,
    aggregation: str = "weighted"
) -> pl.DataFrame

Check every grounding branch against its adiabatic thermal limit.

For the given fault the RMS shield current of every branch (from the most recent :func:groundinsight.run_fault) is converted into the thermally equivalent short-time current I_th = I_s_rms * sqrt(m + n) (IEC 60909-0) and compared against the admissible adiabatic current I_adm = k * S / sqrt(t_k) (IEC 60949) of the branch conductor. A branch is checked only if its :class:BranchType carries both conductor_material and cross_section_mm2; the others are reported with within_limit = None.

.. note:: The thermal data is opt-in. A network that declares none of it still solves and is still reported here — every branch appears with its currents (I_s_rms_A, i_p_A, I_th_A) and only the judgement columns are None. Declaring half of it (a material without a cross-section or the reverse) is the one case that is announced with a logging.WARNING, because such a row is indistinguishable from a passing one at a glance.

All three IEC 60909 inputs — t_k, kappa and n — can either be passed explicitly or be left to the network data:

  • t_k falls back to Fault.t_k_s,
  • n falls back to Fault.n_factor (default 1.0),
  • kappa falls back to the current-weighted aggregation over the feeding sources' kappa / r_to_x, see :func:~groundinsight.analysis.shortcircuit.resolve_fault_sc_characteristics.

An explicit argument always wins over the stored value, so a study can override a single scenario without editing the network.

Parameters:

Name Type Description Default
network Network

A network that already has results for fault_name (call :func:run_fault first).

 required
fault_name str

Name of the solved fault whose branch currents are checked.

 required
t_k float

Fault duration T_k in seconds. Must be strictly positive. Defaults to Fault.t_k_s of the named fault.

 None
kappa float

Peak factor kappa in (1, 2]. Mutually exclusive with r_to_x. If neither is given, kappa is resolved from the sources' IEC 60909 data.

 None
r_to_x float

R/X ratio at the fault, converted to kappa internally. Mutually exclusive with kappa.

 None
n float

AC-decay heat factor n (IEC 60909-0). Defaults to Fault.n_factor, i.e. 1.0 (far-from-generator; I_k'' = I_k).

 None
f float

System frequency in Hz used in the m factor. Defaults to the lowest positive frequency of network (usually 50), or 50.0.

 None
aggregation (weighted, max)

Only relevant when kappa is resolved from the sources. "weighted" reproduces the sum of the individual source peaks exactly; "max" takes the largest source kappa and is conservative. Ignored when kappa or r_to_x is given.

 'weighted'

Returns:

Type Description
DataFrame

One row per branch with the columns branch_name, I_s_rms_A, i_p_A, kappa, m, n, t_k_s, I_th_factor, I_th_A, material, cross_section_mm2, k, I_admissible_A, utilization (I_th / I_adm) and within_limit (bool or None when the branch has no thermal parameters). i_p_A = kappa * sqrt(2) * I_s_rms is the peak current of that branch and is the input for a later electrodynamic-force check.

Raises:

Type Description
ValueError

If t_k is neither given nor stored and positive, if both kappa and r_to_x are given, if kappa can be resolved from neither the arguments nor the sources, or if network has no results for fault_name. Also if the stored result does not cover every branch of the network -- an incomplete result would silently report the missing branches as free of violations. Re-run :func:groundinsight.run_fault in that case.

Examples:

>>> gi.check_conductor_limits(net, "F1", t_k=0.5, r_to_x=0.1)
>>> gi.check_conductor_limits(net, "F1")  # all inputs from the model
Source code in src/groundinsight/analysis/thermal.py 
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
def check_conductor_limits(
    network: Network,
    fault_name: str,
    t_k: Optional[float] = None,
    *,
    kappa: Optional[float] = None,
    r_to_x: Optional[float] = None,
    n: Optional[float] = None,
    f: Optional[float] = None,
    aggregation: str = "weighted",
) -> pl.DataFrame:
    """
    Check every grounding branch against its adiabatic thermal limit.

    For the given fault the RMS shield current of every branch (from the
    most recent :func:`groundinsight.run_fault`) is converted into the
    thermally equivalent short-time current ``I_th = I_s_rms * sqrt(m + n)``
    (IEC 60909-0) and compared against the admissible adiabatic current
    ``I_adm = k * S / sqrt(t_k)`` (IEC 60949) of the branch conductor. A
    branch is checked only if its :class:`BranchType` carries both
    ``conductor_material`` and ``cross_section_mm2``; the others are
    reported with ``within_limit = None``.

    .. note::
       The thermal data is **opt-in**. A network that declares none of it
       still solves and is still reported here — every branch appears with
       its currents (``I_s_rms_A``, ``i_p_A``, ``I_th_A``) and only the
       judgement columns are ``None``. Declaring *half* of it
       (a material without a cross-section or the reverse) is the one case
       that is announced with a ``logging.WARNING``, because such a row is
       indistinguishable from a passing one at a glance.

    All three IEC 60909 inputs — ``t_k``, ``kappa`` and ``n`` — can either
    be passed explicitly or be left to the network data:

    * ``t_k`` falls back to ``Fault.t_k_s``,
    * ``n`` falls back to ``Fault.n_factor`` (default ``1.0``),
    * ``kappa`` falls back to the current-weighted aggregation over the
      feeding sources' ``kappa`` / ``r_to_x``, see
      :func:`~groundinsight.analysis.shortcircuit.resolve_fault_sc_characteristics`.

    An explicit argument always wins over the stored value, so a study can
    override a single scenario without editing the network.

    Parameters
    ----------
    network : Network
        A network that already has results for ``fault_name`` (call
        :func:`run_fault` first).
    fault_name : str
        Name of the solved fault whose branch currents are checked.
    t_k : float, optional
        Fault duration ``T_k`` in seconds. Must be strictly positive.
        Defaults to ``Fault.t_k_s`` of the named fault.
    kappa : float, optional
        Peak factor ``kappa`` in ``(1, 2]``. Mutually exclusive with
        ``r_to_x``. If neither is given, ``kappa`` is resolved from the
        sources' IEC 60909 data.
    r_to_x : float, optional
        ``R/X`` ratio at the fault, converted to ``kappa`` internally.
        Mutually exclusive with ``kappa``.
    n : float, optional
        AC-decay heat factor ``n`` (IEC 60909-0). Defaults to
        ``Fault.n_factor``, i.e. ``1.0`` (far-from-generator;
        ``I_k'' = I_k``).
    f : float, optional
        System frequency in Hz used in the ``m`` factor. Defaults to the
        lowest positive frequency of ``network`` (usually 50), or ``50.0``.
    aggregation : {'weighted', 'max'}, default 'weighted'
        Only relevant when ``kappa`` is resolved from the sources.
        ``"weighted"`` reproduces the sum of the individual source peaks
        exactly; ``"max"`` takes the largest source ``kappa`` and is
        conservative. Ignored when ``kappa`` or ``r_to_x`` is given.

    Returns
    -------
    polars.DataFrame
        One row per branch with the columns ``branch_name``,
        ``I_s_rms_A``, ``i_p_A``, ``kappa``, ``m``, ``n``, ``t_k_s``,
        ``I_th_factor``, ``I_th_A``, ``material``, ``cross_section_mm2``,
        ``k``, ``I_admissible_A``, ``utilization`` (``I_th / I_adm``) and
        ``within_limit`` (``bool`` or ``None`` when the branch has no
        thermal parameters). ``i_p_A = kappa * sqrt(2) * I_s_rms`` is the
        peak current of that branch and is the input for a later
        electrodynamic-force check.

    Raises
    ------
    ValueError
        If ``t_k`` is neither given nor stored and positive, if both
        ``kappa`` and ``r_to_x`` are given, if ``kappa`` can be resolved
        from neither the arguments nor the sources, or if ``network`` has
        no results for ``fault_name``. Also if the stored result does not
        cover every branch of the network -- an incomplete result would
        silently report the missing branches as free of violations. Re-run
        :func:`groundinsight.run_fault` in that case.

    Examples
    --------
    >>> gi.check_conductor_limits(net, "F1", t_k=0.5, r_to_x=0.1)  # doctest: +SKIP
    >>> gi.check_conductor_limits(net, "F1")  # all inputs from the model  # doctest: +SKIP
    """
    sc = _resolve_sc_inputs(
        network, fault_name, t_k, kappa, r_to_x, n, f, aggregation
    )
    t_k, kappa, m, n, i_th_factor = (
        sc.t_k, sc.kappa, sc.m, sc.n, sc.i_th_factor
    )

    result = network.results[fault_name]
    # ``compute_branch_currents`` emits a row for *every* branch, including
    # inactive and open ones, so any missing name means the stored result is
    # stale or half-built. See _require_complete_results.
    _require_complete_results(
        fault_name,
        list(network.branches),
        [rb.name for rb in result.branches],
        "branch",
    )
    rows: List[Dict[str, Any]] = []
    incomplete: List[Tuple[str, str]] = []
    for result_branch in result.branches:
        name = result_branch.name
        i_s_rms = float(result_branch.i_s)
        i_th = i_s_rms * i_th_factor

        branch = network.branches.get(name)
        material = None
        cross_section = None
        k_val = None
        i_adm = None
        utilization = None
        within = None
        if branch is not None:
            btype = branch.type
            material = getattr(btype, "conductor_material", None)
            cross_section = getattr(btype, "cross_section_mm2", None)
            if material is not None and cross_section is not None:
                k_val = iec60949_k(
                    material,
                    theta_initial_C=getattr(btype, "theta_initial_C", 20.0),
                    theta_final_C=getattr(btype, "theta_final_C", None),
                )
                i_adm = admissible_short_circuit_current(k_val, cross_section, t_k)
                utilization = i_th / i_adm if i_adm > 0 else None
                within = bool(i_th <= i_adm)
            elif material is not None:
                incomplete.append((name, "cross_section_mm2"))
            elif cross_section is not None:
                incomplete.append((name, "conductor_material"))

        rows.append(
            {
                "branch_name": name,
                "I_s_rms_A": i_s_rms,
                "i_p_A": peak_short_circuit_current(i_s_rms, kappa),
                "kappa": float(kappa),
                "m": float(m),
                "n": float(n),
                "t_k_s": float(t_k),
                "I_th_factor": float(i_th_factor),
                "I_th_A": float(i_th),
                "material": material,
                "cross_section_mm2": cross_section,
                "k": k_val,
                "I_admissible_A": i_adm,
                "utilization": utilization,
                "within_limit": within,
            }
        )

    _warn_half_declared(incomplete, fault_name, "branch(es)")

    violations = [r["branch_name"] for r in rows if r["within_limit"] is False]
    if violations:
        logger.warning(
            "Thermal limit exceeded on %d branch(es) for fault '%s' "
            "(t_k=%.3g s, kappa=%.3f): %s",
            len(violations), fault_name, t_k, kappa, violations,
        )

    return pl.DataFrame(rows, schema=_BRANCH_SCHEMA)

check_node_limits 

check_node_limits(
    network: Network,
    fault_name: str,
    t_k: Optional[float] = None,
    *,
    kappa: Optional[float] = None,
    r_to_x: Optional[float] = None,
    n: Optional[float] = None,
    f: Optional[float] = None,
    aggregation: str = "weighted",
    elements: Sequence[str] = (
        "earthing_conductor",
        "electrode",
    )
) -> pl.DataFrame

Check every bus against the adiabatic thermal limits of its grounding.

The node-side counterpart of :func:check_conductor_limits. Where the branch check assesses the shield / earth wire between buses, this one assesses the two elements at a bus, which EN 50522 / IEC 61936-1 keep strictly apart because they carry different currents:

earthing_conductor The earthing conductor (Erdungsleiter): the lumped connection that brings the earth-fault current from the installation into the grounding system. It carries the full injected current ResultBus.i_inj — at a source bus the source infeed, at the fault bus the total fault current — and is therefore the more heavily stressed of the two.

electrode The earth electrode (Erder): the buried part, which only carries the share actually dissipated into the soil at this bus, ResultBus.ia = u_EPR / Z_B. In a well-meshed system that is a small fraction of the injection, so sizing the electrode for the full fault current is wasteful — and sizing the earthing conductor for the electrode current is dangerous.

Both are reported separately, one row each, in a long-format frame.

An element is assessed only if its :class:BusType carries both the *_material and the *_cross_section_mm2 field of that prefix; otherwise the row is reported with within_limit = None.

.. note:: The thermal data is opt-in. A network that declares none of it still solves and is still reported here — every active bus appears with its currents and only the judgement columns are None. Declaring half of it (a material without a cross-section or the reverse) is the one case that is announced with a logging.WARNING, because such a row is indistinguishable from a passing one at a glance.

Current split

Each element carries a free factor current_split in (0, 1] (BusType.earthing_conductor_current_split / BusType.electrode_current_split), applied as I_conductor = I_rms * current_split. It expresses how the bus current divides between physically parallel paths that the nodal model lumps into one:

  • 1.0 (default) — one conductor / one electrode carries everything.
  • 1/N — N equal parallel paths, e.g. 0.25 for four down-conductors from a substation steelwork to the grid.
  • 0.5 — a ring electrode fed at a single point: the current splits into the two halves of the ring.
  • any IEEE Std 80 style split / division factor, entered directly.

The factor is deliberately not derived automatically: the split depends on the physical arrangement inside the substation, which the nodal grounding model does not resolve. Values above 1 are rejected at the model level — more current than the bus carries is an error, not a split.

Superposition

Identical to the branch check: the frequency-domain solve superposes the linear AC RMS currents, and the non-linear IEC 60909 factors kappa / m are applied once to the aggregate. i_p and I_th are never superposed directly.

Parameters:

Name Type Description Default
network Network

A network that already has results for fault_name (call :func:run_fault first).

 required
fault_name str

Name of the solved fault whose bus currents are checked.

 required
t_k float

Fault duration T_k in seconds. Defaults to Fault.t_k_s.

 None
kappa float

Peak factor kappa in (1, 2]. Mutually exclusive with r_to_x. Resolved from the sources when neither is given.

 None
r_to_x float

R/X ratio at the fault, converted to kappa internally. Mutually exclusive with kappa.

 None
n float

AC-decay heat factor n (IEC 60909-0). Defaults to Fault.n_factor, i.e. 1.0.

 None
f float

System frequency in Hz used in the m factor. Defaults to the lowest positive frequency of network, or 50.0.

 None
aggregation (weighted, max)

Only relevant when kappa is resolved from the sources; see :func:check_conductor_limits.

 'weighted'
elements sequence of str

Which of the two elements to report. Pass a single-element sequence to restrict the frame.

 ``('earthing_conductor', 'electrode')``

Returns:

Type Description
DataFrame

Long format, one row per bus and requested element, with the columns bus_name, element, I_rms_A (the bus current before the split), current_split, I_conductor_A, i_p_A, kappa, m, n, t_k_s, I_th_factor, I_th_A, material, cross_section_mm2, k, I_admissible_A, utilization and within_limit. within_limit is None where the BusType declares no material / cross-section for that element — the current columns are still filled, so an undimensioned bus shows what it would have to carry.

Raises:

Type Description
ValueError

For an unknown entry in elements, and for the same conditions as :func:check_conductor_limits (no results, missing t_k, unresolvable kappa, both kappa and r_to_x given, or a stored result that does not cover every active bus).
Notes

A bus that is both the fault bus and the only source bus gets i_inj = 0: the current never enters the grounding system, it circulates in the phase conductor. That is physically correct and not a missing value.

Results computed before ResultBus.i_inj existed default to 0 for the injection. Re-run :func:run_fault if an earthing-conductor row reads I_rms_A = 0 at a bus that should be carrying current.

Examples:

>>> gi.run_fault(net, "F1")
>>> gi.check_node_limits(net, "F1").filter(
...     pl.col("element") == "electrode"
... )
Source code in src/groundinsight/analysis/thermal.py 
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
def check_node_limits(
    network: Network,
    fault_name: str,
    t_k: Optional[float] = None,
    *,
    kappa: Optional[float] = None,
    r_to_x: Optional[float] = None,
    n: Optional[float] = None,
    f: Optional[float] = None,
    aggregation: str = "weighted",
    elements: Sequence[str] = ("earthing_conductor", "electrode"),
) -> pl.DataFrame:
    """
    Check every bus against the adiabatic thermal limits of its grounding.

    The node-side counterpart of :func:`check_conductor_limits`. Where the
    branch check assesses the shield / earth wire *between* buses, this one
    assesses the two elements *at* a bus, which EN 50522 / IEC 61936-1 keep
    strictly apart because they carry different currents:

    ``earthing_conductor``
        The **earthing conductor** (*Erdungsleiter*): the lumped connection
        that brings the earth-fault current from the installation into the
        grounding system. It carries the full injected current
        ``ResultBus.i_inj`` — at a source bus the source infeed, at the
        fault bus the total fault current — and is therefore the more
        heavily stressed of the two.

    ``electrode``
        The **earth electrode** (*Erder*): the buried part, which only
        carries the share actually dissipated into the soil at this bus,
        ``ResultBus.ia = u_EPR / Z_B``. In a well-meshed system that is a
        small fraction of the injection, so sizing the electrode for the
        full fault current is wasteful — and sizing the earthing conductor
        for the electrode current is dangerous.

    Both are reported separately, one row each, in a long-format frame.

    An element is assessed only if its :class:`BusType` carries both the
    ``*_material`` and the ``*_cross_section_mm2`` field of that prefix;
    otherwise the row is reported with ``within_limit = None``.

    .. note::
       The thermal data is **opt-in**. A network that declares none of it
       still solves and is still reported here — every active bus appears
       with its currents and only the judgement columns are ``None``.
       Declaring *half* of it (a material without a cross-section or the
       reverse) is the one case that is announced with a
       ``logging.WARNING``, because such a row is indistinguishable from a
       passing one at a glance.

    Current split
    -------------
    Each element carries a free factor ``current_split`` in ``(0, 1]``
    (``BusType.earthing_conductor_current_split`` /
    ``BusType.electrode_current_split``), applied as
    ``I_conductor = I_rms * current_split``. It expresses how the bus
    current divides between physically parallel paths that the nodal model
    lumps into one:

    * ``1.0`` (default) — one conductor / one electrode carries everything.
    * ``1/N`` — ``N`` equal parallel paths, e.g. ``0.25`` for four
      down-conductors from a substation steelwork to the grid.
    * ``0.5`` — a ring electrode fed at a single point: the current splits
      into the two halves of the ring.
    * any IEEE Std 80 style split / division factor, entered directly.

    The factor is deliberately *not* derived automatically: the split
    depends on the physical arrangement inside the substation, which the
    nodal grounding model does not resolve. Values above ``1`` are rejected
    at the model level — more current than the bus carries is an error, not
    a split.

    Superposition
    -------------
    Identical to the branch check: the frequency-domain solve superposes
    the *linear* AC RMS currents, and the *non-linear* IEC 60909 factors
    ``kappa`` / ``m`` are applied once to the aggregate. ``i_p`` and
    ``I_th`` are never superposed directly.

    Parameters
    ----------
    network : Network
        A network that already has results for ``fault_name`` (call
        :func:`run_fault` first).
    fault_name : str
        Name of the solved fault whose bus currents are checked.
    t_k : float, optional
        Fault duration ``T_k`` in seconds. Defaults to ``Fault.t_k_s``.
    kappa : float, optional
        Peak factor ``kappa`` in ``(1, 2]``. Mutually exclusive with
        ``r_to_x``. Resolved from the sources when neither is given.
    r_to_x : float, optional
        ``R/X`` ratio at the fault, converted to ``kappa`` internally.
        Mutually exclusive with ``kappa``.
    n : float, optional
        AC-decay heat factor ``n`` (IEC 60909-0). Defaults to
        ``Fault.n_factor``, i.e. ``1.0``.
    f : float, optional
        System frequency in Hz used in the ``m`` factor. Defaults to the
        lowest positive frequency of ``network``, or ``50.0``.
    aggregation : {'weighted', 'max'}, default 'weighted'
        Only relevant when ``kappa`` is resolved from the sources; see
        :func:`check_conductor_limits`.
    elements : sequence of str, default ``('earthing_conductor', 'electrode')``
        Which of the two elements to report. Pass a single-element
        sequence to restrict the frame.

    Returns
    -------
    polars.DataFrame
        Long format, one row per bus and requested element, with the
        columns ``bus_name``, ``element``, ``I_rms_A`` (the bus current
        before the split), ``current_split``, ``I_conductor_A``,
        ``i_p_A``, ``kappa``, ``m``, ``n``, ``t_k_s``, ``I_th_factor``,
        ``I_th_A``, ``material``, ``cross_section_mm2``, ``k``,
        ``I_admissible_A``, ``utilization`` and ``within_limit``.
        ``within_limit`` is ``None`` where the ``BusType`` declares no
        material / cross-section for that element — the current columns
        are still filled, so an undimensioned bus shows what it would have
        to carry.

    Raises
    ------
    ValueError
        For an unknown entry in ``elements``, and for the same conditions
        as :func:`check_conductor_limits` (no results, missing ``t_k``,
        unresolvable ``kappa``, both ``kappa`` and ``r_to_x`` given, or a
        stored result that does not cover every *active* bus).

    Notes
    -----
    A bus that is *both* the fault bus and the only source bus gets
    ``i_inj = 0``: the current never enters the grounding system, it
    circulates in the phase conductor. That is physically correct and not
    a missing value.

    Results computed before ``ResultBus.i_inj`` existed default to ``0``
    for the injection. Re-run :func:`run_fault` if an earthing-conductor
    row reads ``I_rms_A = 0`` at a bus that should be carrying current.

    Examples
    --------
    >>> gi.run_fault(net, "F1")                       # doctest: +SKIP
    >>> gi.check_node_limits(net, "F1").filter(       # doctest: +SKIP
    ...     pl.col("element") == "electrode"
    ... )
    """
    unknown = [e for e in elements if e not in _NODE_ELEMENTS]
    if unknown:
        raise ValueError(
            f"Unknown node element(s) {unknown}. Known elements: "
            f"{sorted(_NODE_ELEMENTS)}."
        )

    sc = _resolve_sc_inputs(
        network, fault_name, t_k, kappa, r_to_x, n, f, aggregation
    )

    result = network.results[fault_name]
    # ``solve_network`` emits a row for every *active* bus (inactive ones are
    # removed from the nodal system and carry no EPR), so the expectation is
    # the set of active buses. See _require_complete_results.
    _require_complete_results(
        fault_name,
        [name for name, bus in network.buses.items() if bus.active],
        [rb.name for rb in result.buses],
        "bus",
    )
    rows: List[Dict[str, Any]] = []
    incomplete: List[Tuple[str, str]] = []
    for result_bus in result.buses:
        name = result_bus.name
        bus = network.buses.get(name)
        btype = getattr(bus, "type", None)

        for element in elements:
            spec = _NODE_ELEMENTS[element]
            prefix = spec["prefix"]

            i_rms = float(getattr(result_bus, spec["current_attr"], 0.0) or 0.0)

            split = 1.0
            material = None
            cross_section = None
            k_val = None
            i_adm = None
            utilization = None
            within = None
            if btype is not None:
                split = float(
                    getattr(btype, f"{prefix}current_split", 1.0) or 1.0
                )
                material = getattr(btype, f"{prefix}material", None)
                cross_section = getattr(btype, f"{prefix}cross_section_mm2", None)

            i_conductor = i_rms * split
            i_th = i_conductor * sc.i_th_factor

            if material is not None and cross_section is not None:
                k_val = iec60949_k(
                    material,
                    theta_initial_C=getattr(btype, f"{prefix}theta_initial_C", 20.0),
                    theta_final_C=getattr(btype, f"{prefix}theta_final_C", None),
                )
                i_adm = admissible_short_circuit_current(
                    k_val, cross_section, sc.t_k
                )
                utilization = i_th / i_adm if i_adm > 0 else None
                within = bool(i_th <= i_adm)
            elif material is not None:
                incomplete.append((f"{name}/{element}", f"{prefix}cross_section_mm2"))
            elif cross_section is not None:
                incomplete.append((f"{name}/{element}", f"{prefix}material"))

            rows.append(
                {
                    "bus_name": name,
                    "element": element,
                    "I_rms_A": i_rms,
                    "current_split": split,
                    "I_conductor_A": i_conductor,
                    "i_p_A": peak_short_circuit_current(i_conductor, sc.kappa),
                    "kappa": sc.kappa,
                    "m": sc.m,
                    "n": sc.n,
                    "t_k_s": sc.t_k,
                    "I_th_factor": sc.i_th_factor,
                    "I_th_A": float(i_th),
                    "material": material,
                    "cross_section_mm2": cross_section,
                    "k": k_val,
                    "I_admissible_A": i_adm,
                    "utilization": utilization,
                    "within_limit": within,
                }
            )

    _warn_half_declared(incomplete, fault_name, "node element(s)")

    violations = [
        f"{r['bus_name']}/{r['element']}" for r in rows if r["within_limit"] is False
    ]
    if violations:
        logger.warning(
            "Thermal limit exceeded on %d node element(s) for fault '%s' "
            "(t_k=%.3g s, kappa=%.3f): %s",
            len(violations), fault_name, sc.t_k, sc.kappa, violations,
        )

    return pl.DataFrame(rows, schema=_NODE_SCHEMA)

classify 

classify(
    frame: pl.DataFrame,
    value: str,
    edges: Sequence[float],
    *,
    labels: Optional[Sequence[str]] = None,
    name: str = "class"
) -> pl.DataFrame

Add a class column by binning a numeric column at the given edges.

edges are the interior boundaries: n edges make n + 1 classes. Bins are closed on the right, so a value exactly on an edge falls into the lower class -- the conservative reading when the edge is a limit.

Parameters:

Name Type Description Default
frame DataFrame

Frame to extend.

 required
value str

Numeric column to classify.

 required
edges sequence of float

Strictly increasing interior boundaries, e.g. [80, 150] for the three bands "below 80", "80 to 150" and "above 150".

 required
labels sequence of str

Names for the classes, len(edges) + 1 of them. Defaults to readable interval labels built from the edges.

 None
name str

Name of the added column. Defaults to "class".

 'class'

Returns:

Type Description
DataFrame

The input with one column added. Rows whose value is null get a null class rather than being forced into the lowest band.

Raises:

Type Description
ValueError

If the column is missing or not numeric, if edges is empty or not strictly increasing, or if labels has the wrong length.

Examples:

>>> classify(study.buses(), "EPR_V", [80.0, 150.0],
...          labels=["ok", "check", "exceeded"])
Source code in src/groundinsight/analysis/statistics.py 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
def classify(
    frame: pl.DataFrame,
    value: str,
    edges: Sequence[float],
    *,
    labels: Optional[Sequence[str]] = None,
    name: str = "class",
) -> pl.DataFrame:
    """
    Add a class column by binning a numeric column at the given edges.

    ``edges`` are the interior boundaries: ``n`` edges make ``n + 1`` classes.
    Bins are closed on the right, so a value exactly on an edge falls into the
    lower class -- the conservative reading when the edge is a limit.

    Parameters
    ----------
    frame : pl.DataFrame
        Frame to extend.
    value : str
        Numeric column to classify.
    edges : sequence of float
        Strictly increasing interior boundaries, e.g. ``[80, 150]`` for the
        three bands "below 80", "80 to 150" and "above 150".
    labels : sequence of str, optional
        Names for the classes, ``len(edges) + 1`` of them. Defaults to readable
        interval labels built from the edges.
    name : str, optional
        Name of the added column. Defaults to ``"class"``.

    Returns
    -------
    pl.DataFrame
        The input with one column added. Rows whose value is null get a null
        class rather than being forced into the lowest band.

    Raises
    ------
    ValueError
        If the column is missing or not numeric, if ``edges`` is empty or not
        strictly increasing, or if ``labels`` has the wrong length.

    Examples
    --------
    >>> classify(study.buses(), "EPR_V", [80.0, 150.0],  # doctest: +SKIP
    ...          labels=["ok", "check", "exceeded"])
    """
    if value not in frame.columns:
        raise ValueError(
            f"Column '{value}' is not in the frame. Available: {frame.columns}."
        )
    if not frame.schema[value].is_numeric():
        raise ValueError(
            f"Column '{value}' has dtype {frame.schema[value]} and cannot be "
            f"binned numerically."
        )
    edges = [float(e) for e in edges]
    if not edges:
        raise ValueError(
            "classify needs at least one edge; with none there is only one "
            "class and nothing to decide."
        )
    if any(b <= a for a, b in zip(edges, edges[1:])):
        raise ValueError(
            f"Edges must be strictly increasing, got {edges}. Overlapping or "
            f"repeated edges would make the class assignment ambiguous."
        )

    if labels is None:
        labels = [f"<= {edges[0]:g}"]
        labels += [f"{a:g} - {b:g}" for a, b in zip(edges, edges[1:])]
        labels.append(f"> {edges[-1]:g}")
    labels = list(labels)
    if len(labels) != len(edges) + 1:
        raise ValueError(
            f"{len(edges)} edge(s) define {len(edges) + 1} classes, but "
            f"{len(labels)} label(s) were given."
        )

    expression = pl.when(pl.col(value).is_null()).then(None)
    for edge, label in zip(edges, labels):
        expression = expression.when(pl.col(value) <= edge).then(pl.lit(label))
    expression = expression.otherwise(pl.lit(labels[-1]))
    return frame.with_columns(expression.alias(name))

evaluate_max_epr_under_k 

evaluate_max_epr_under_k(
    network: Network,
    bus_names: List[str],
    k: KVector,
    *,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None
) -> Dict[str, float]

Evaluate the RMS EPR at each selected bus for a given k vector.

For every name in bus_names the bus grounding impedance is overwritten with the rho-f linear form Z(rho, f) = k1*rho + (k2+jk3)*f + (k4+jk5)*rho*f — using the bus' own :attr:Bus.specific_earth_resistance for rho — and the bus is swept as the active fault. An existing fault at that bus is reused; otherwise a temporary one is created with fault_scalings (or 1.0 at every simulation frequency by default) and removed at the end. The original bus impedances, the temporary faults and the previous active_fault are restored in a finally block, so the network is left exactly as it was on entry.

Parameters:

Name Type Description Default
network Network

The simulation network. Sources, branches, paths and the buses to be swept must already be configured.

 required
bus_names list of str

Buses whose impedance is rewritten and which are swept as fault locations one by one. Must be non-empty and refer to buses in network.

 required
k tuple of (float, float, float, float, float)

Real model parameters (k1, k2, k3, k4, k5).

 required
fault_scalings dict of float to float

Frequency-resolved scalings used for any fault that has to be created on the fly (no pre-existing fault at the swept bus). Defaults to {f: 1.0} for every simulation frequency. Existing faults are reused unmodified.

 None
run_fault_kwargs dict

Extra keyword arguments forwarded to :func:run_fault at every step.

 None

Returns:

Type Description
dict of str to float

Mapping of each bus_name to its RMS EPR in volts at that bus when it is the active fault.

Raises:

Type Description
ValueError

If bus_names is empty, k does not have length 5, or any name is not in the network.
Source code in src/groundinsight/analysis/inverse_rho_f.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
def evaluate_max_epr_under_k(
    network: Network,
    bus_names: List[str],
    k: KVector,
    *,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None,
) -> Dict[str, float]:
    """Evaluate the RMS EPR at each selected bus for a given ``k`` vector.

    For every name in ``bus_names`` the bus grounding impedance is
    overwritten with the rho-f linear form
    ``Z(rho, f) = k1*rho + (k2+jk3)*f + (k4+jk5)*rho*f`` — using the
    bus' own :attr:`Bus.specific_earth_resistance` for ``rho`` — and
    the bus is swept as the active fault. An existing fault at that bus
    is reused; otherwise a temporary one is created with
    ``fault_scalings`` (or ``1.0`` at every simulation frequency by
    default) and removed at the end. The original bus impedances, the
    temporary faults and the previous ``active_fault`` are restored in a
    ``finally`` block, so the network is left exactly as it was on entry.

    Parameters
    ----------
    network : Network
        The simulation network. Sources, branches, paths and the buses to
        be swept must already be configured.
    bus_names : list of str
        Buses whose impedance is rewritten and which are swept as fault
        locations one by one. Must be non-empty and refer to buses in
        ``network``.
    k : tuple of (float, float, float, float, float)
        Real model parameters ``(k1, k2, k3, k4, k5)``.
    fault_scalings : dict of float to float, optional
        Frequency-resolved scalings used for any fault that has to be
        created on the fly (no pre-existing fault at the swept bus).
        Defaults to ``{f: 1.0}`` for every simulation frequency. Existing
        faults are reused unmodified.
    run_fault_kwargs : dict, optional
        Extra keyword arguments forwarded to :func:`run_fault` at every
        step.

    Returns
    -------
    dict of str to float
        Mapping of each ``bus_name`` to its RMS EPR in volts at that bus
        when it is the active fault.

    Raises
    ------
    ValueError
        If ``bus_names`` is empty, ``k`` does not have length 5, or any
        name is not in the network.
    """
    if not bus_names:
        raise ValueError("bus_names must not be empty.")
    if len(k) != 5:
        raise ValueError(f"k must be a 5-tuple, got length {len(k)}.")
    missing = [b for b in bus_names if b not in network.buses]
    if missing:
        raise ValueError(f"Unknown bus(es) in network: {missing!r}.")

    rfk: Dict[str, Any] = dict(run_fault_kwargs) if run_fault_kwargs else {}
    if fault_scalings is None:
        fault_scalings = {float(f): 1.0 for f in network.frequencies}

    # Snapshot current bus impedances so the network can be restored.
    impedance_backup: Dict[str, Dict[float, ComplexNumber]] = {
        b: dict(network.buses[b].impedance) for b in bus_names
    }

    # Look up an existing fault per swept bus (if any).
    existing_faults_by_bus: Dict[str, str] = {}
    for fname, fault in network.faults.items():
        if fault.bus in bus_names and fault.bus not in existing_faults_by_bus:
            existing_faults_by_bus[fault.bus] = fname

    temp_faults_created: List[str] = []
    active_fault_backup = network.active_fault
    # Snapshot results of any *pre-existing* fault we are about to reuse so we
    # can restore them (run_fault overwrites results[fault] with the k-form EPR).
    reused_fault_names = set(existing_faults_by_bus.values())
    results_backup = {fn: network.results.get(fn) for fn in reused_fault_names}
    had_result_backup = {fn: fn in network.results for fn in reused_fault_names}

    # Snapshot existing paths and clear them: ``run_fault`` only triggers
    # ``define_paths`` when ``network.paths`` is empty. We need fresh paths
    # for every (source, fault) combination -- including the temporary
    # faults we are about to create -- so we drop the cache for the
    # duration of the sweep and restore it at the end.
    paths_backup = dict(network.paths)
    network.paths.clear()

    def _temp_fault_name(b: str) -> str:
        base = f"_inv_rhof_{b}"
        candidate = base
        n = 0
        while candidate in network.faults:
            n += 1
            candidate = f"{base}_{n}"
        return candidate

    epr_per_bus: Dict[str, float] = {}

    try:
        # Overwrite every selected bus impedance with the k-form.
        for b in bus_names:
            bus = network.buses[b]
            rho = float(bus.specific_earth_resistance)
            new_imp: Dict[float, ComplexNumber] = {}
            for f in network.frequencies:
                z = _z_rho_f(k, rho, float(f))
                new_imp[float(f)] = ComplexNumber(real=z.real, imag=z.imag)
            bus.impedance = new_imp

        # Pre-create any temporary faults the sweep needs *before* the
        # first ``run_fault`` so ``define_paths`` enumerates paths for
        # all swept (source, fault) combinations in a single pass.
        sweep_fault_names: List[str] = []
        for b in bus_names:
            if b in existing_faults_by_bus:
                sweep_fault_names.append(existing_faults_by_bus[b])
            else:
                fname = _temp_fault_name(b)
                network.add_fault(
                    Fault(name=fname, bus=b, scalings=dict(fault_scalings))
                )
                temp_faults_created.append(fname)
                sweep_fault_names.append(fname)

        # Sweep: each fault once. The first call populates ``network.paths``
        # for every (source, fault) pair; subsequent calls reuse the cache.
        for b, fname in zip(bus_names, sweep_fault_names):
            run_fault(network, fault_name=fname, **rfk)
            try:
                result_bus = next(
                    rb for rb in network.results[fname].buses if rb.name == b
                )
            except StopIteration as exc:
                # ``run_fault_kwargs={"buses": [...]}`` may filter the
                # swept bus out of the result frame. Surface the situation
                # with a clear lookup error rather than a bare
                # ``StopIteration`` from the generator expression.
                available = [
                    rb.name for rb in network.results[fname].buses
                ]
                raise LookupError(
                    f"Bus {b!r} is not present in the result frame of "
                    f"fault {fname!r}. Available buses: {available!r}. "
                    "If you are passing run_fault_kwargs={'buses': [...]}, "
                    "make sure every swept bus appears in that list."
                ) from exc
            epr_per_bus[b] = float(result_bus.uepr)
    finally:
        # Restore bus impedances regardless of success/failure.
        for b in bus_names:
            network.buses[b].impedance = impedance_backup[b]
        # Drop temporary faults and any results they produced.
        for fname in temp_faults_created:
            network.faults.pop(fname, None)
            network.results.pop(fname, None)
        # Restore results for reused (pre-existing) faults; drop any result
        # run_fault created for a reused fault that had none before.
        for fname in reused_fault_names:
            if had_result_backup.get(fname):
                network.results[fname] = results_backup[fname]
            else:
                network.results.pop(fname, None)
        # Restore the path cache (cleared at the start of the sweep).
        # Use an atomic dict swap rather than ``clear()`` + ``update()`` so
        # that concurrent readers never see an empty ``network.paths``.
        # Pydantic's ``BaseModel`` exposes the field as a mutable dict; we
        # mutate it in place but in two steps that look atomic from a
        # snapshot perspective (rebuild a new dict, then reassign).
        rebuilt = dict(paths_backup)
        try:
            network.paths = rebuilt  # atomic rebind when allowed
        except Exception:
            network.paths.clear()
            network.paths.update(rebuilt)
        # Restore active_fault and the per-fault _active flags exactly. When
        # the network had no active fault on entry (the common fresh-network
        # case), clear the flag rather than leaving it on a deleted temp fault.
        if active_fault_backup is None:
            network.active_fault = None
            for _flt in network.faults.values():
                _flt._set_active(False)
        elif active_fault_backup in network.faults:
            network.set_active_fault(active_fault_backup, keep_results=True)
        else:  # pragma: no cover -- backup fault vanished
            network.active_fault = None
            for _flt in network.faults.values():
                _flt._set_active(False)

    return epr_per_bus

final_temperature 

final_temperature(material: str, covering: str) -> float

Look up a tabulated maximum final temperature theta_f.

Convenience accessor for :data:FINAL_TEMPERATURES, meant to be used when building a BranchType or BusType so that the source of the number stays visible in the model::

gi.BusType(..., electrode_theta_final_C=gi.final_temperature("Steel", "bare"))

Which standard the answer comes from depends on the covering: EN 50522 Table 2 for an uninsulated conductor, IEC 60364-5-54 Table 54.2 for an insulated one. The insulated caps are much lower, so the two must not be confused — hence one accessor over one catalogue rather than a per-source lookup the caller has to choose between.

Parameters:

Name Type Description Default
material str

One of the keys of :data:IEC60949_MATERIALS ("Cu", "Al", "Steel").

 required
covering str

Covering / surface treatment key. Uninsulated: "bare", "tinned" (copper only) or "galvanized" (steel only). Insulated: "PVC", "XLPE" or "EPR". Only combinations that occur in practice are tabulated — there is no tinned aluminium and no galvanised copper — and "PE" is not tabulated at all.

 required

Returns:

Type Description
float

Maximum permissible final temperature in °C.

Raises:

Type Description
ValueError

If the material or the material / covering combination is not tabulated. The message lists what is available — the table is intentionally incomplete rather than populated with unsourced values, since an over-estimated theta_f permits too much current.
See Also

CABLE_INSULATION_LIMITS : the insulated subset, on its own, when the conductor material is not known or not relevant.

Examples:

>>> final_temperature("Cu", "bare")
300.0
>>> final_temperature("Cu", "tinned")
150.0
>>> final_temperature("Cu", "PVC")
160.0
Source code in src/groundinsight/analysis/thermal.py 
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
def final_temperature(material: str, covering: str) -> float:
    """
    Look up a tabulated maximum final temperature ``theta_f``.

    Convenience accessor for :data:`FINAL_TEMPERATURES`, meant to be used when
    *building* a ``BranchType`` or ``BusType`` so that the source of the
    number stays visible in the model::

        gi.BusType(..., electrode_theta_final_C=gi.final_temperature("Steel", "bare"))

    Which standard the answer comes from depends on the covering: EN 50522
    Table 2 for an uninsulated conductor, IEC 60364-5-54 Table 54.2 for an
    insulated one. The insulated caps are much lower, so the two must not be
    confused — hence one accessor over one catalogue rather than a per-source
    lookup the caller has to choose between.

    Parameters
    ----------
    material : str
        One of the keys of :data:`IEC60949_MATERIALS` (``"Cu"``, ``"Al"``,
        ``"Steel"``).
    covering : str
        Covering / surface treatment key. Uninsulated: ``"bare"``,
        ``"tinned"`` (copper only) or ``"galvanized"`` (steel only).
        Insulated: ``"PVC"``, ``"XLPE"`` or ``"EPR"``. Only combinations that
        occur in practice are tabulated — there is no tinned aluminium and no
        galvanised copper — and ``"PE"`` is not tabulated at all.

    Returns
    -------
    float
        Maximum permissible final temperature in °C.

    Raises
    ------
    ValueError
        If the material or the material / covering combination is not
        tabulated. The message lists what *is* available — the table is
        intentionally incomplete rather than populated with unsourced
        values, since an over-estimated ``theta_f`` permits too much
        current.

    See Also
    --------
    CABLE_INSULATION_LIMITS : the insulated subset, on its own, when the
        conductor material is not known or not relevant.

    Examples
    --------
    >>> final_temperature("Cu", "bare")
    300.0
    >>> final_temperature("Cu", "tinned")
    150.0
    >>> final_temperature("Cu", "PVC")
    160.0
    """
    if material not in FINAL_TEMPERATURES:
        raise ValueError(
            f"No tabulated final temperatures for material {material!r}. "
            f"Known materials: {sorted(FINAL_TEMPERATURES)}."
        )
    coverings = FINAL_TEMPERATURES[material]
    if covering not in coverings:
        raise ValueError(
            f"No tabulated final temperature for covering {covering!r} on "
            f"{material!r}. Tabulated coverings: {sorted(coverings)}. Add the "
            "value from EN 50522 Table 2 (uninsulated) or IEC 60364-5-54 "
            "Table 54.2 (insulated) to "
            "groundinsight.analysis.thermal.FINAL_TEMPERATURES, or pass "
            "theta_final_C explicitly."
        )
    return float(coverings[covering])

find_max_rho_f_scaling 

find_max_rho_f_scaling(
    network: Network,
    bus_names: List[str],
    u_limit: float,
    k_ref: KVector,
    *,
    c_bounds: Tuple[float, float] = (0.001, 1000.0),
    tol_rel: float = 0.001,
    max_iter: int = 60,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None
) -> Dict[str, Any]

Find the largest scaling factor c of a reference k-vector.

Performs a log-scale bisection on c so that k = c * k_ref yields a maximum RMS EPR (across all buses swept as faults) that is at or below u_limit. k_ref is typically obtained by fitting the rho-f model above to measured impedance points; this function answers how much head-room the network has around that characteristic.

Parameters:

Name Type Description Default
network Network

The simulation network. Buses, branches, sources and paths must already be configured.

 required
bus_names list of str

Buses whose impedance is parameterised by k and which are swept as fault locations.

 required
u_limit float

Upper bound on the RMS EPR (in volts) at any swept bus. Must be finite and strictly positive.

 required
k_ref tuple of (float, float, float, float, float)

Reference 5-tuple (k1_ref, ..., k5_ref). Must have length 5 and not be the zero vector.

 required
c_bounds tuple of (float, float)

Search interval for c. Finite, strictly positive, c_bounds[0] < c_bounds[1]. Defaults to (1e-3, 1e3).

 (0.001, 1000.0)
tol_rel float

Bisection tolerance on the relative bracket width (c_hi - c_lo) / c_lo. Must be finite and strictly positive. Defaults to 1e-3.

 0.001
max_iter int

Hard cap on bisection steps. Must be an int of at least 1. Defaults to 60.

 60
fault_scalings dict of float to float

See :func:evaluate_max_epr_under_k.

 None
run_fault_kwargs dict

Forwarded to :func:run_fault.

 None

Returns:

Type Description
dict

Mapping with keys

  • "c_max" (float): a scaling factor whose swept maximum EPR was evaluated and found compatible with u_limit. It is the largest such factor only when "converged" is True.
  • "k_max" (tuple of float): c_max * k_ref.
  • "max_epr_rms_at_c_max" (float): maximum RMS EPR across all swept buses at c_max, in volts.
  • "epr_rms_per_bus_at_c_max" (dict of str to float): RMS EPR per swept bus at c_max.
  • "iterations" (int): number of bisection steps taken.
  • "converged" (bool): True iff the bracket was closed to within tol_rel. Check this before using c_max or k_max as a design value.
  • "status" (str): "converged", "bracket_within_tol_on_entry", "max_iter_reached" or "bracket_fully_admissible".
  • "c_bracket" (tuple of float): the interval provably containing the true maximum; (c_hi, inf) when the whole bracket was admissible.
  • "bracket_rel_width" (float): (c_hi - c_lo) / c_lo of that interval, comparable against tol_rel.

As in :func:~groundinsight.analysis.find_max_rho_scaling, iterations == 0 does not identify a case on its own; read "status".

Raises:

Type Description
ValueError

For invalid input -- u_limit, tol_rel, max_iter and c_bounds are checked for finiteness as well as range -- or if c = c_bounds[0] already violates the EPR limit, or if the model returns a non-finite EPR at some trial factor.

Examples:

>>> import groundinsight as gi
>>> res = gi.find_max_rho_f_scaling(
...     network=net, bus_names=["b0", "b1"], u_limit=200.0,
...     k_ref=(0.01, 0.0, 0.0, 0.0, 0.0),
... )
Source code in src/groundinsight/analysis/inverse_rho_f.py 
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
def find_max_rho_f_scaling(
    network: Network,
    bus_names: List[str],
    u_limit: float,
    k_ref: KVector,
    *,
    c_bounds: Tuple[float, float] = (1e-3, 1e3),
    tol_rel: float = 1e-3,
    max_iter: int = 60,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None,
) -> Dict[str, Any]:
    """Find the largest scaling factor ``c`` of a reference k-vector.

    Performs a log-scale bisection on ``c`` so that ``k = c * k_ref``
    yields a maximum RMS EPR (across all buses swept as faults) that is
    at or below ``u_limit``. ``k_ref`` is typically obtained by fitting
    the rho-f model above to measured impedance points; this function
    answers *how much head-room the network has around that
    characteristic*.

    Parameters
    ----------
    network : Network
        The simulation network. Buses, branches, sources and paths must
        already be configured.
    bus_names : list of str
        Buses whose impedance is parameterised by ``k`` *and* which are
        swept as fault locations.
    u_limit : float
        Upper bound on the RMS EPR (in volts) at any swept bus. Must be
        finite and strictly positive.
    k_ref : tuple of (float, float, float, float, float)
        Reference 5-tuple ``(k1_ref, ..., k5_ref)``. Must have length 5
        and not be the zero vector.
    c_bounds : tuple of (float, float), optional
        Search interval for ``c``. Finite, strictly positive,
        ``c_bounds[0] < c_bounds[1]``. Defaults to ``(1e-3, 1e3)``.
    tol_rel : float, optional
        Bisection tolerance on the relative bracket width
        ``(c_hi - c_lo) / c_lo``. Must be finite and strictly positive.
        Defaults to ``1e-3``.
    max_iter : int, optional
        Hard cap on bisection steps. Must be an ``int`` of at least 1.
        Defaults to ``60``.
    fault_scalings : dict of float to float, optional
        See :func:`evaluate_max_epr_under_k`.
    run_fault_kwargs : dict, optional
        Forwarded to :func:`run_fault`.

    Returns
    -------
    dict
        Mapping with keys

        - ``"c_max"`` (float): a scaling factor whose swept maximum EPR
          was evaluated and found compatible with ``u_limit``. It is the
          *largest* such factor only when ``"converged"`` is ``True``.
        - ``"k_max"`` (tuple of float): ``c_max * k_ref``.
        - ``"max_epr_rms_at_c_max"`` (float): maximum RMS EPR across all
          swept buses at ``c_max``, in volts.
        - ``"epr_rms_per_bus_at_c_max"`` (dict of str to float): RMS EPR
          per swept bus at ``c_max``.
        - ``"iterations"`` (int): number of bisection steps taken.
        - ``"converged"`` (bool): ``True`` iff the bracket was closed to
          within ``tol_rel``. Check this before using ``c_max`` or
          ``k_max`` as a design value.
        - ``"status"`` (str): ``"converged"``,
          ``"bracket_within_tol_on_entry"``, ``"max_iter_reached"`` or
          ``"bracket_fully_admissible"``.
        - ``"c_bracket"`` (tuple of float): the interval provably
          containing the true maximum; ``(c_hi, inf)`` when the whole
          bracket was admissible.
        - ``"bracket_rel_width"`` (float): ``(c_hi - c_lo) / c_lo`` of
          that interval, comparable against ``tol_rel``.

        As in :func:`~groundinsight.analysis.find_max_rho_scaling`,
        ``iterations == 0`` does not identify a case on its own; read
        ``"status"``.

    Raises
    ------
    ValueError
        For invalid input -- ``u_limit``, ``tol_rel``, ``max_iter`` and
        ``c_bounds`` are checked for finiteness as well as range -- or if
        ``c = c_bounds[0]`` already violates the EPR limit, or if the
        model returns a non-finite EPR at some trial factor.

    Examples
    --------
    >>> import groundinsight as gi  # doctest: +SKIP
    >>> res = gi.find_max_rho_f_scaling(  # doctest: +SKIP
    ...     network=net, bus_names=["b0", "b1"], u_limit=200.0,
    ...     k_ref=(0.01, 0.0, 0.0, 0.0, 0.0),
    ... )
    """
    validate_limit(u_limit, "u_limit")
    validate_tol_rel(tol_rel)
    validate_max_iter(max_iter)
    if not bus_names:
        raise ValueError("bus_names must not be empty.")
    if len(k_ref) != 5:
        raise ValueError(f"k_ref must be a 5-tuple, got length {len(k_ref)}.")
    if all(x == 0 for x in k_ref):
        raise ValueError("k_ref must not be the zero vector.")
    c_lo_init, c_hi_init = validate_c_bounds(c_bounds)

    last_epr_per_bus: Dict[str, Dict[str, float]] = {"value": {}}

    def _max_epr_at(c: float) -> float:
        eprs = evaluate_max_epr_under_k(
            network, bus_names,
            k=tuple(c * x for x in k_ref),
            fault_scalings=fault_scalings,
            run_fault_kwargs=run_fault_kwargs,
        )
        last_epr_per_bus["value"] = eprs
        epr = max(eprs.values()) if eprs else 0.0
        if not math.isfinite(epr):
            # Same reasoning as in ``inverse_rho``: a non-finite EPR makes
            # every comparison against u_limit False, so the bisection would
            # take the same turn at every step and reach the lower bracket
            # bound without ever raising.
            raise ValueError(
                f"The swept maximum EPR evaluated to {epr!r} for the scaling "
                f"factor c={c:g} (k = {tuple(c * x for x in k_ref)!r}). A "
                "non-finite EPR cannot be compared against u_limit -- every "
                "comparison would be False and the bisection would walk to "
                "the lower bracket bound without raising. Check the k-vector "
                "at that scale against the buses' specific_earth_resistance."
            )
        return epr

    epr_lo = _max_epr_at(c_lo_init)
    eprs_lo_snapshot: Dict[str, float] = dict(last_epr_per_bus["value"])
    if epr_lo > u_limit:
        raise ValueError(
            f"u_limit={u_limit:g} V is below the maximum EPR at "
            f"c={c_lo_init:g}: max EPR_RMS={epr_lo:g} V — no scaling "
            f"factor in the bracket {c_bounds!r} satisfies the constraint."
        )

    epr_hi = _max_epr_at(c_hi_init)
    iterations = 0
    if epr_hi <= u_limit:
        logger.info(
            "Bracket fully admissible: max EPR_RMS at c=%g is %g V <= "
            "u_limit=%g V. Consider widening c_bounds.",
            c_hi_init, epr_hi, u_limit,
        )
        c_max, epr_at = c_hi_init, epr_hi
        epr_per_bus_at = dict(last_epr_per_bus["value"])
        status = STATUS_BRACKET_FULLY_ADMISSIBLE
        c_lo, c_hi = c_lo_init, c_hi_init
    else:
        c_lo, c_hi = c_lo_init, c_hi_init
        while iterations < max_iter and (c_hi - c_lo) / c_lo > tol_rel:
            c_mid = math.sqrt(c_lo * c_hi)  # geometric mean -> log bisection
            epr_mid = _max_epr_at(c_mid)
            if epr_mid <= u_limit:
                c_lo, epr_lo = c_mid, epr_mid
                eprs_lo_snapshot = dict(last_epr_per_bus["value"])
            else:
                c_hi = c_mid
            iterations += 1
        c_max, epr_at = c_lo, epr_lo
        epr_per_bus_at = eprs_lo_snapshot
        status = classify(iterations, c_lo, c_hi, tol_rel)
        if status == STATUS_MAX_ITER_REACHED:
            logger.warning(
                "Bisection stopped at the step cap max_iter=%d without "
                "closing the bracket: c in [%g, %g], relative width %g > "
                "tol_rel=%g. c_max=%g is admissible but may be well below "
                "the true maximum.",
                max_iter, c_lo, c_hi, (c_hi - c_lo) / c_lo, tol_rel, c_lo,
            )

    # Defensive consistency check: the maximum of the per-bus EPRs at
    # c_max must equal the reported ``max_epr_rms_at_c_max``. A drift
    # here would mean the snapshot and the headline figure refer to
    # different iterations — historically a bug surface (the
    # ``eprs_lo_snapshot`` shadow variable was added in 0.4.0 to plug
    # exactly this gap). Re-evaluate once at ``c_max`` to lock the
    # invariant in place; the cost is one extra ``run_fault`` per sweep.
    epr_at = _max_epr_at(c_max)
    epr_per_bus_at = dict(last_epr_per_bus["value"])

    return {
        "c_max": c_max,
        "k_max": tuple(c_max * x for x in k_ref),
        "max_epr_rms_at_c_max": epr_at,
        "epr_rms_per_bus_at_c_max": epr_per_bus_at,
        "iterations": iterations,
        **report(status, c_lo, c_hi),
    }

find_max_rho_scaling 

find_max_rho_scaling(
    network: Network,
    fault_name: str,
    bus_names: List[str],
    u_max: float,
    *,
    c_bounds: Tuple[float, float] = (0.001, 1000.0),
    tol_rel: float = 0.001,
    max_iter: int = 60,
    run_fault_kwargs: Dict[str, Any] = None
) -> Dict[str, Any]

Find the largest uniform rho-scaling factor compatible with an EPR limit.

Bisects the scalar c such that scaling every selected bus' specific earth resistance to c * rho_0 yields an RMS earth potential rise at the fault bus that just satisfies |u_EPR|_rms <= u_max. The bus impedance formula is re-evaluated through the existing :meth:Bus.calculate_impedance machinery, so any user-defined rho-f characteristic is honoured.

Parameters:

Name Type Description Default
network Network

The simulation network. Must already contain the named fault, the sources, the buses listed in bus_names and consistent paths from sources to fault.

 required
fault_name str

Name of the fault to evaluate. Used as active_fault during every bisection step.

 required
bus_names list of str

Names of the buses whose specific earth resistance is uniformly scaled by the same factor c. Must be non-empty and refer to buses in network.

 required
u_max float

Upper bound on the RMS earth potential rise at the fault bus, in volts. Must be finite and strictly positive. The RMS is taken over all simulation frequencies, matching :attr:ResultBus.uepr.

 required
c_bounds tuple of (float, float)

Search interval for the scaling factor c. Both bounds must be finite and strictly positive and c_bounds[0] < c_bounds[1]. Defaults to (1e-3, 1e3), i.e. six decades.

 (0.001, 1000.0)
tol_rel float

Relative tolerance on the bracket width (c_hi - c_lo) / c_lo at which the bisection terminates. Must be finite and strictly positive. Defaults to 1e-3.

 0.001
max_iter int

Hard cap on the number of bisection steps. Must be an int of at least 1. Defaults to 60, which is roughly four times what the default bracket and tolerance need.

 60
run_fault_kwargs dict

Extra keyword arguments forwarded to :func:run_fault at every step (e.g. {"auto_parallel_coefficients": True}). Defaults to None.

 None

Returns:

Type Description
dict

Mapping with keys

  • "c_max" (float): a scaling factor whose EPR was evaluated and found to satisfy the constraint. This is a guarantee in one direction only: c_max is always admissible, but it is the largest admissible factor only when "converged" is True.
  • "u_epr_rms_at_c_max" (float): the RMS EPR at the fault bus evaluated at c_max, in volts.
  • "rho_max" (dict of str to float): c_max * rho_0[bus] for every selected bus.
  • "iterations" (int): number of bisection steps taken.
  • "converged" (bool): True iff the search closed the bracket to within tol_rel around the threshold. Check this before using c_max as a design value.
  • "status" (str): which stopping condition applied -- "converged", "bracket_within_tol_on_entry", "max_iter_reached" or "bracket_fully_admissible". See :mod:groundinsight.analysis._bisection.
  • "c_bracket" (tuple of float): the interval that provably contains the true maximum admissible factor. For "bracket_fully_admissible" this is (c_hi, inf): nothing above the bracket was ever evaluated, so widen c_bounds.
  • "bracket_rel_width" (float): (c_hi - c_lo) / c_lo of that interval, directly comparable against tol_rel, and inf for "bracket_fully_admissible".

iterations == 0 on its own does not identify a case: it is produced by a fully admissible bracket, by a bracket that was already narrower than tol_rel, and (before the guards below existed) by a cap of zero steps -- outcomes whose c_max came from opposite ends of the bracket. Use "status".

Raises:

Type Description
ValueError

If u_max, tol_rel, max_iter or c_bounds is invalid (see the parameter descriptions -- non-finite values are rejected, not just out-of-range ones), bus_names is empty, any name is not in the network, the EPR at the lower bound c_bounds[0] already exceeds u_max, or the model returns a non-finite EPR at some trial factor.

Examples:

>>> import groundinsight as gi
>>> res = gi.find_max_rho_scaling(
...     network=net, fault_name="flt",
...     bus_names=["b0", "b1"], u_max=200.0,
... )
>>> if not res["converged"]:
...     print(res["status"], res["c_bracket"])
Source code in src/groundinsight/analysis/inverse_rho.py 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
def find_max_rho_scaling(
    network: Network,
    fault_name: str,
    bus_names: List[str],
    u_max: float,
    *,
    c_bounds: Tuple[float, float] = (1e-3, 1e3),
    tol_rel: float = 1e-3,
    max_iter: int = 60,
    run_fault_kwargs: Dict[str, Any] = None,
) -> Dict[str, Any]:
    """
    Find the largest uniform rho-scaling factor compatible with an EPR limit.

    Bisects the scalar ``c`` such that scaling every selected bus'
    specific earth resistance to ``c * rho_0`` yields an RMS earth
    potential rise at the fault bus that just satisfies
    ``|u_EPR|_rms <= u_max``. The bus impedance formula is re-evaluated
    through the existing :meth:`Bus.calculate_impedance` machinery, so
    any user-defined rho-f characteristic is honoured.

    Parameters
    ----------
    network : Network
        The simulation network. Must already contain the named fault, the
        sources, the buses listed in ``bus_names`` and consistent paths
        from sources to fault.
    fault_name : str
        Name of the fault to evaluate. Used as ``active_fault`` during
        every bisection step.
    bus_names : list of str
        Names of the buses whose specific earth resistance is uniformly
        scaled by the same factor ``c``. Must be non-empty and refer to
        buses in ``network``.
    u_max : float
        Upper bound on the RMS earth potential rise at the fault bus, in
        volts. Must be finite and strictly positive. The RMS is taken
        over all simulation frequencies, matching :attr:`ResultBus.uepr`.
    c_bounds : tuple of (float, float), optional
        Search interval for the scaling factor ``c``. Both bounds must be
        finite and strictly positive and ``c_bounds[0] < c_bounds[1]``.
        Defaults to ``(1e-3, 1e3)``, i.e. six decades.
    tol_rel : float, optional
        Relative tolerance on the bracket width
        ``(c_hi - c_lo) / c_lo`` at which the bisection terminates. Must
        be finite and strictly positive. Defaults to ``1e-3``.
    max_iter : int, optional
        Hard cap on the number of bisection steps. Must be an ``int``
        of at least 1. Defaults to ``60``, which is roughly four times
        what the default bracket and tolerance need.
    run_fault_kwargs : dict, optional
        Extra keyword arguments forwarded to :func:`run_fault` at every
        step (e.g. ``{"auto_parallel_coefficients": True}``). Defaults to
        ``None``.

    Returns
    -------
    dict
        Mapping with keys

        - ``"c_max"`` (float): a scaling factor whose EPR was evaluated
          and found to satisfy the constraint. **This is a guarantee in
          one direction only**: ``c_max`` is always admissible, but it is
          the *largest* admissible factor only when ``"converged"`` is
          ``True``.
        - ``"u_epr_rms_at_c_max"`` (float): the RMS EPR at the fault bus
          evaluated at ``c_max``, in volts.
        - ``"rho_max"`` (dict of str to float): ``c_max * rho_0[bus]`` for
          every selected bus.
        - ``"iterations"`` (int): number of bisection steps taken.
        - ``"converged"`` (bool): ``True`` iff the search closed the
          bracket to within ``tol_rel`` around the threshold. Check this
          before using ``c_max`` as a design value.
        - ``"status"`` (str): which stopping condition applied --
          ``"converged"``, ``"bracket_within_tol_on_entry"``,
          ``"max_iter_reached"`` or ``"bracket_fully_admissible"``. See
          :mod:`groundinsight.analysis._bisection`.
        - ``"c_bracket"`` (tuple of float): the interval that provably
          contains the true maximum admissible factor. For
          ``"bracket_fully_admissible"`` this is ``(c_hi, inf)``: nothing
          above the bracket was ever evaluated, so widen ``c_bounds``.
        - ``"bracket_rel_width"`` (float): ``(c_hi - c_lo) / c_lo`` of
          that interval, directly comparable against ``tol_rel``, and
          ``inf`` for ``"bracket_fully_admissible"``.

        ``iterations == 0`` on its own does **not** identify a case: it is
        produced by a fully admissible bracket, by a bracket that was
        already narrower than ``tol_rel``, and (before the guards below
        existed) by a cap of zero steps -- outcomes whose ``c_max`` came
        from opposite ends of the bracket. Use ``"status"``.

    Raises
    ------
    ValueError
        If ``u_max``, ``tol_rel``, ``max_iter`` or ``c_bounds`` is
        invalid (see the parameter descriptions -- non-finite values are
        rejected, not just out-of-range ones), ``bus_names`` is empty,
        any name is not in the network, the EPR at the lower bound
        ``c_bounds[0]`` already exceeds ``u_max``, or the model returns a
        non-finite EPR at some trial factor.

    Examples
    --------
    >>> import groundinsight as gi  # doctest: +SKIP
    >>> res = gi.find_max_rho_scaling(  # doctest: +SKIP
    ...     network=net, fault_name="flt",
    ...     bus_names=["b0", "b1"], u_max=200.0,
    ... )
    >>> if not res["converged"]:  # doctest: +SKIP
    ...     print(res["status"], res["c_bracket"])
    """
    validate_limit(u_max, "u_max")
    validate_tol_rel(tol_rel)
    validate_max_iter(max_iter)
    if not bus_names:
        raise ValueError("bus_names must not be empty.")
    c_lo_init, c_hi_init = validate_c_bounds(c_bounds)
    missing = [b for b in bus_names if b not in network.buses]
    if missing:
        raise ValueError(f"Unknown bus(es) in network: {missing!r}.")
    if fault_name not in network.faults:
        raise ValueError(f"Unknown fault {fault_name!r} in network.")

    rfk: Dict[str, Any] = dict(run_fault_kwargs) if run_fault_kwargs else {}

    fault_bus_name = network.faults[fault_name].bus

    # Snapshot original rhos so the network can be restored at the end.
    rho_0: Dict[str, float] = {
        b: float(network.buses[b].specific_earth_resistance) for b in bus_names
    }
    # Snapshot the state run_fault will mutate so the search leaves the network
    # exactly as it found it (the returned figures are already final).
    active_fault_backup = network.active_fault
    result_backup = network.results.get(fault_name)
    had_result = fault_name in network.results

    def _epr_rms_at(c: float) -> float:
        """Evaluate the RMS EPR at the fault bus for a given scaling factor."""
        for b in bus_names:
            bus = network.buses[b]
            bus.specific_earth_resistance = c * rho_0[b]
            bus.calculate_impedance(network.frequencies)
        run_fault(network, fault_name=fault_name, **rfk)
        result_bus = next(
            rb
            for rb in network.results[fault_name].buses
            if rb.name == fault_bus_name
        )
        epr = float(result_bus.uepr)
        if not math.isfinite(epr):
            # A non-finite EPR cannot be compared against u_max: every
            # comparison is False, so the bisection would take the same turn
            # at every step and walk silently to the lower bracket bound.
            # The impedance pipeline raises on NaN before a formula can get
            # this far today, so this is a second lock on the same door --
            # it defends the search below it, not the model above it.
            raise ValueError(
                f"The EPR at bus {fault_bus_name!r} evaluated to {epr!r} for "
                f"the scaling factor c={c:g}. A non-finite EPR cannot be "
                "compared against u_max -- every comparison would be False "
                "and the bisection would walk to the lower bracket bound "
                "without ever raising. This is a model problem, not a search "
                "problem: check the bus impedance formula at the scaled rho."
            )
        return epr

    iterations = 0
    try:
        c_lo, c_hi = c_lo_init, c_hi_init
        epr_lo = _epr_rms_at(c_lo)
        if epr_lo > u_max:
            raise ValueError(
                f"u_max={u_max:g} V is below the EPR at c_min={c_lo:g}: "
                f"|u_EPR|_rms(c_min)={epr_lo:g} V — no scaling factor in "
                f"the bracket {c_bounds!r} satisfies the constraint."
            )
        epr_hi = _epr_rms_at(c_hi)
        if epr_hi <= u_max:
            # Whole bracket is admissible. c_hi is a real, measured answer,
            # but it is a lower bound on the true maximum, not the maximum:
            # nothing above c_hi was evaluated. The status says so, and
            # ``c_bracket`` comes back as (c_hi, inf) to make that
            # machine-readable.
            logger.info(
                "Bracket fully admissible: |u_EPR|_rms(c_hi=%g)=%g V <= "
                "u_max=%g V. Consider widening c_bounds.",
                c_hi, epr_hi, u_max,
            )
            c_max, epr_at = c_hi, epr_hi
            status = STATUS_BRACKET_FULLY_ADMISSIBLE
        else:
            # From here on epr(c_lo) <= u_max < epr(c_hi) holds and is
            # preserved by every step, so c_lo is always a *verified*
            # admissible factor and the threshold stays bracketed.
            while iterations < max_iter and (c_hi - c_lo) / c_lo > tol_rel:
                c_mid = math.sqrt(c_lo * c_hi)  # geometric mean -> log bisection
                epr_mid = _epr_rms_at(c_mid)
                if epr_mid <= u_max:
                    c_lo, epr_lo = c_mid, epr_mid
                else:
                    c_hi, epr_hi = c_mid, epr_mid
                iterations += 1
            c_max, epr_at = c_lo, epr_lo
            status = classify(iterations, c_lo, c_hi, tol_rel)
            # Ask the classifier rather than re-deriving the same condition:
            # two copies of "did it close?" is exactly how a status and its
            # log message drift apart.
            if status == STATUS_MAX_ITER_REACHED:
                logger.warning(
                    "Bisection stopped at the step cap max_iter=%d without "
                    "closing the bracket: c in [%g, %g], relative width %g > "
                    "tol_rel=%g. c_max=%g is admissible but may be well below "
                    "the true maximum.",
                    max_iter, c_lo, c_hi, (c_hi - c_lo) / c_lo, tol_rel, c_lo,
                )
    finally:
        # Restore original rhos and recompute their impedances no matter what.
        for b in bus_names:
            bus = network.buses[b]
            bus.specific_earth_resistance = rho_0[b]
            bus.calculate_impedance(network.frequencies)
        # Restore the result cache and active fault so no trace of the search
        # remains on the network.
        if had_result:
            network.results[fault_name] = result_backup
        else:
            network.results.pop(fault_name, None)
        if active_fault_backup is None:
            network.active_fault = None
            for _flt in network.faults.values():
                _flt._set_active(False)
        elif active_fault_backup in network.faults:
            network.set_active_fault(active_fault_backup, keep_results=True)

    return {
        "c_max": c_max,
        "u_epr_rms_at_c_max": epr_at,
        "rho_max": {b: c_max * rho_0[b] for b in bus_names},
        "iterations": iterations,
        **report(status, c_lo, c_hi),
    }

iec60909_m 

iec60909_m(kappa: float, f: float, t_k: float) -> float

Heat-effect factor m of the aperiodic (DC) short-circuit component.

m = (exp(4*f*Tk*ln(kappa-1)) - 1) / (2*f*Tk*ln(kappa-1)) per IEC 60909-0. The limits are handled explicitly: kappa -> 1 (no DC offset) gives m -> 0; kappa -> 2 (non-decaying DC, X/R -> inf) gives m -> 2.

Notes

The kappa -> 2 limit is worth spelling out because it is easy to get backwards. Substituting a = ln(kappa - 1) -> 0 and expanding, (exp(4 f Tk a) - 1) / (2 f Tk a) -> 4 f Tk a / (2 f Tk a) = 2. A vanishing-resistance fault therefore carries the maximum DC heat, not none. pandapower's _calc_ith sets m = 0 for kappa > 1.99, which is why :func:groundinsight.io.pandapower_sc.read_shortcircuit_results recomputes I_th with this function instead of copying pandapower's ith_ka.

Parameters:

Name Type Description Default
kappa float

Peak factor kappa in (1, 2] (see :func:kappa_from_r_to_x).

 required
f float

System frequency in Hz (50 or 60).

 required
t_k float

Fault duration T_k in seconds. Must be strictly positive.

 required

Returns:

Type Description
float

The dimensionless factor m (>= 0).

Raises:

Type Description
ValueError

If t_k or f is not strictly positive.

Examples:

>>> iec60909_m(1.0, 50.0, 0.5)
0.0
>>> iec60909_m(2.0, 50.0, 0.5)
2.0
Source code in src/groundinsight/analysis/shortcircuit.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
def iec60909_m(kappa: float, f: float, t_k: float) -> float:
    """
    Heat-effect factor ``m`` of the aperiodic (DC) short-circuit component.

    ``m = (exp(4*f*Tk*ln(kappa-1)) - 1) / (2*f*Tk*ln(kappa-1))`` per
    IEC 60909-0. The limits are handled explicitly: ``kappa -> 1`` (no DC
    offset) gives ``m -> 0``; ``kappa -> 2`` (non-decaying DC, ``X/R -> inf``)
    gives ``m -> 2``.

    Notes
    -----
    The ``kappa -> 2`` limit is worth spelling out because it is easy to get
    backwards. Substituting ``a = ln(kappa - 1) -> 0`` and expanding,
    ``(exp(4 f Tk a) - 1) / (2 f Tk a) -> 4 f Tk a / (2 f Tk a) = 2``. A
    vanishing-resistance fault therefore carries the *maximum* DC heat, not
    none. pandapower's ``_calc_ith`` sets ``m = 0`` for ``kappa > 1.99``,
    which is why :func:`groundinsight.io.pandapower_sc.read_shortcircuit_results`
    recomputes ``I_th`` with this function instead of copying pandapower's
    ``ith_ka``.

    Parameters
    ----------
    kappa : float
        Peak factor ``kappa`` in ``(1, 2]`` (see :func:`kappa_from_r_to_x`).
    f : float
        System frequency in Hz (50 or 60).
    t_k : float
        Fault duration ``T_k`` in seconds. Must be strictly positive.

    Returns
    -------
    float
        The dimensionless factor ``m`` (``>= 0``).

    Raises
    ------
    ValueError
        If ``t_k`` or ``f`` is not strictly positive.

    Examples
    --------
    >>> iec60909_m(1.0, 50.0, 0.5)
    0.0
    >>> iec60909_m(2.0, 50.0, 0.5)
    2.0
    """
    if t_k <= 0:
        raise ValueError(f"t_k must be strictly positive, got {t_k!r}.")
    if f <= 0:
        raise ValueError(f"f must be strictly positive, got {f!r}.")
    if kappa <= 1.0:
        return 0.0
    if kappa >= 2.0:
        return 2.0
    a = math.log(kappa - 1.0)  # negative for 1 < kappa < 2
    denom = 2.0 * f * t_k * a
    if abs(denom) < 1e-12:  # kappa extremely close to 2 -> non-decaying DC
        return 2.0
    return (math.exp(4.0 * f * t_k * a) - 1.0) / denom

iec60949_k 

iec60949_k(
    material: str,
    theta_initial_C: float = 20.0,
    theta_final_C: Optional[float] = None,
) -> float

Material constant k of the adiabatic short-circuit equation.

k = K * sqrt(ln((theta_f + beta) / (theta_i + beta))) per IEC 60949 / IEC 60364-5-54, with the base constant K and beta taken from :data:IEC60949_MATERIALS.

Parameters:

Name Type Description Default
material str

One of the keys of :data:IEC60949_MATERIALS ("Cu", "Al", "Steel").

 required
theta_initial_C float

Initial conductor temperature theta_i in °C. Defaults to 20.0 (ambient). For a conductor that also carries load current use its maximum continuous operating temperature.

 20.0
theta_final_C float

Final (maximum permissible) conductor temperature theta_f in °C. Defaults to the material's bare-earthing-conductor value from :data:IEC60949_MATERIALS.

 None

Returns:

Type Description
float

The material constant k in A·s^0.5/mm².

Raises:

Type Description
ValueError

If material is unknown or theta_final_C <= theta_initial_C.

Examples:

>>> round(iec60949_k("Cu", theta_initial_C=90.0, theta_final_C=250.0))
143
Source code in src/groundinsight/analysis/thermal.py 
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
def iec60949_k(
    material: str,
    theta_initial_C: float = 20.0,
    theta_final_C: Optional[float] = None,
) -> float:
    """
    Material constant ``k`` of the adiabatic short-circuit equation.

    ``k = K * sqrt(ln((theta_f + beta) / (theta_i + beta)))`` per IEC 60949 /
    IEC 60364-5-54, with the base constant ``K`` and ``beta`` taken from
    :data:`IEC60949_MATERIALS`.

    Parameters
    ----------
    material : str
        One of the keys of :data:`IEC60949_MATERIALS` (``"Cu"``, ``"Al"``,
        ``"Steel"``).
    theta_initial_C : float
        Initial conductor temperature ``theta_i`` in °C. Defaults to
        ``20.0`` (ambient). For a conductor that also carries load current
        use its maximum continuous operating temperature.
    theta_final_C : float, optional
        Final (maximum permissible) conductor temperature ``theta_f`` in °C.
        Defaults to the material's bare-earthing-conductor value from
        :data:`IEC60949_MATERIALS`.

    Returns
    -------
    float
        The material constant ``k`` in A·s^0.5/mm².

    Raises
    ------
    ValueError
        If ``material`` is unknown or ``theta_final_C <= theta_initial_C``.

    Examples
    --------
    >>> round(iec60949_k("Cu", theta_initial_C=90.0, theta_final_C=250.0))
    143
    """
    if material not in IEC60949_MATERIALS:
        raise ValueError(
            f"Unknown conductor material {material!r}. Known materials: "
            f"{sorted(IEC60949_MATERIALS)}."
        )
    data = IEC60949_MATERIALS[material]
    theta_f = (
        data["theta_final_default_C"] if theta_final_C is None else float(theta_final_C)
    )
    beta = data["beta"]
    if theta_f <= theta_initial_C:
        raise ValueError(
            f"theta_final_C ({theta_f}) must exceed theta_initial_C "
            f"({theta_initial_C}) for material {material!r}."
        )
    return data["K"] * math.sqrt(
        math.log((theta_f + beta) / (theta_initial_C + beta))
    )

kappa_from_r_to_x 

kappa_from_r_to_x(r_to_x: float) -> float

Peak factor kappa from the R/X ratio (IEC 60909-0).

kappa = 1.02 + 0.98 * exp(-3 * R/X), bounded to (1, 2].

The R/X to insert is the ratio of the fault loop, not of the positive-sequence system alone. For a three-phase fault the two coincide; for the line-to-earth fault that dominates grounding studies the loop impedance is 2*Z1 + Z0, hence R/X = (2*R1 + R0) / (2*X1 + X0) (see :func:groundinsight.io.pandapower_sc.read_shortcircuit_results).

Parameters:

Name Type Description Default
r_to_x float

The R/X ratio of the fault loop. Must be non-negative.

 required

Returns:

Type Description
float

The peak factor kappa.

Raises:

Type Description
ValueError

If r_to_x is negative.

Examples:

>>> round(kappa_from_r_to_x(0.0), 3)
2.0
>>> round(kappa_from_r_to_x(0.1), 4)
1.746
Source code in src/groundinsight/analysis/shortcircuit.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def kappa_from_r_to_x(r_to_x: float) -> float:
    """
    Peak factor ``kappa`` from the ``R/X`` ratio (IEC 60909-0).

    ``kappa = 1.02 + 0.98 * exp(-3 * R/X)``, bounded to ``(1, 2]``.

    The ``R/X`` to insert is the ratio of the **fault loop**, not of the
    positive-sequence system alone. For a three-phase fault the two
    coincide; for the line-to-earth fault that dominates grounding studies
    the loop impedance is ``2*Z1 + Z0``, hence
    ``R/X = (2*R1 + R0) / (2*X1 + X0)``
    (see :func:`groundinsight.io.pandapower_sc.read_shortcircuit_results`).

    Parameters
    ----------
    r_to_x : float
        The ``R/X`` ratio of the fault loop. Must be non-negative.

    Returns
    -------
    float
        The peak factor ``kappa``.

    Raises
    ------
    ValueError
        If ``r_to_x`` is negative.

    Examples
    --------
    >>> round(kappa_from_r_to_x(0.0), 3)
    2.0
    >>> round(kappa_from_r_to_x(0.1), 4)
    1.746
    """
    if r_to_x < 0:
        raise ValueError(f"r_to_x must be non-negative, got {r_to_x!r}.")
    return 1.02 + 0.98 * math.exp(-3.0 * r_to_x)

peak_short_circuit_current 

peak_short_circuit_current(
    i_k: float, kappa: float
) -> float

Peak short-circuit current i_p = kappa * sqrt(2) * I_k'' (IEC 60909-0).

i_p drives the electrodynamic (mechanical) stress on conductors and supports; the thermal counterpart is :func:thermal_equivalent_current.

Parameters:

Name Type Description Default
i_k float

RMS short-circuit current I_k'' in amperes. Must be non-negative.

 required
kappa float

Peak factor from :func:kappa_from_r_to_x.

 required

Returns:

Type Description
float

The peak current in amperes.

Raises:

Type Description
ValueError

If i_k is negative.

Examples:

>>> round(peak_short_circuit_current(1000.0, 1.8), 2)
2545.58
Source code in src/groundinsight/analysis/shortcircuit.py 
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
def peak_short_circuit_current(i_k: float, kappa: float) -> float:
    """
    Peak short-circuit current ``i_p = kappa * sqrt(2) * I_k''`` (IEC 60909-0).

    ``i_p`` drives the *electrodynamic* (mechanical) stress on conductors and
    supports; the thermal counterpart is
    :func:`thermal_equivalent_current`.

    Parameters
    ----------
    i_k : float
        RMS short-circuit current ``I_k''`` in amperes. Must be
        non-negative.
    kappa : float
        Peak factor from :func:`kappa_from_r_to_x`.

    Returns
    -------
    float
        The peak current in amperes.

    Raises
    ------
    ValueError
        If ``i_k`` is negative.

    Examples
    --------
    >>> round(peak_short_circuit_current(1000.0, 1.8), 2)
    2545.58
    """
    if i_k < 0:
        raise ValueError(f"i_k must be non-negative, got {i_k!r}.")
    return kappa * math.sqrt(2.0) * i_k

resolve_fault_sc_characteristics 

resolve_fault_sc_characteristics(
    network: Network,
    fault_name: str,
    *,
    frequency: Optional[float] = None,
    aggregation: str = "weighted"
) -> FaultShortCircuitData

Aggregate the IEC 60909 data of a fault and its sources into one set of characteristics.

This is the executable form of the project's superposition rule (see the module docstring): the linear RMS currents stay superposed by the solve, and one effective kappa is resolved here so that the non-linear 60909 factors are applied exactly once, to the aggregate.

Parameters:

Name Type Description Default
network Network

The network holding the fault and its sources.

 required
fault_name str

Name of the fault to resolve.

 required
frequency float

Frequency in Hz at which the source injections are weighted. Defaults to the lowest positive frequency of network (the power frequency in a harmonic study), or 50.0.

 None
aggregation (weighted, max)

"weighted" (default) takes the current-weighted mean of the source kappa values, which reproduces the sum of the individual peak currents exactly. "max" takes the largest kappa, a strictly conservative bound.

 'weighted'

Returns:

Type Description
FaultShortCircuitData

The resolved characteristics. kappa is None when no contributing source carries 60909 data -- callers must then fall back to an explicit argument.

Raises:

Type Description
ValueError

If fault_name is unknown or aggregation is not one of the supported modes.

Examples:

>>> import groundinsight as gi
>>> net = gi.create_network(name="demo", frequencies=[50.0])
>>> bt = gi.BusType(name="b", system_type="s", voltage_level=20.0,
...                 impedance_formula="rho*0 + 1.0")
>>> _ = gi.create_bus(name="B1", type=bt, network=net)
>>> _ = gi.create_source(name="S1", bus="B1", values={50.0: 1000.0},
...                      network=net, r_to_x=0.1)
>>> _ = gi.create_fault(name="F1", bus="B1", scalings={50.0: 1.0},
...                     network=net, t_k_s=0.5)
>>> data = gi.resolve_fault_sc_characteristics(net, "F1")
>>> round(data.kappa, 4)
1.746
Source code in src/groundinsight/analysis/shortcircuit.py 
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
def resolve_fault_sc_characteristics(
    network: Network,
    fault_name: str,
    *,
    frequency: Optional[float] = None,
    aggregation: str = "weighted",
) -> FaultShortCircuitData:
    """
    Aggregate the IEC 60909 data of a fault and its sources into one set of
    characteristics.

    This is the executable form of the project's superposition rule (see the
    module docstring): the *linear* RMS currents stay superposed by the
    solve, and one effective ``kappa`` is resolved here so that the
    *non-linear* 60909 factors are applied exactly once, to the aggregate.

    Parameters
    ----------
    network : Network
        The network holding the fault and its sources.
    fault_name : str
        Name of the fault to resolve.
    frequency : float, optional
        Frequency in Hz at which the source injections are weighted.
        Defaults to the lowest positive frequency of ``network`` (the power
        frequency in a harmonic study), or ``50.0``.
    aggregation : {'weighted', 'max'}
        ``"weighted"`` (default) takes the current-weighted mean of the
        source ``kappa`` values, which reproduces the sum of the individual
        peak currents exactly. ``"max"`` takes the largest ``kappa``, a
        strictly conservative bound.

    Returns
    -------
    FaultShortCircuitData
        The resolved characteristics. ``kappa`` is ``None`` when no
        contributing source carries 60909 data -- callers must then fall
        back to an explicit argument.

    Raises
    ------
    ValueError
        If ``fault_name`` is unknown or ``aggregation`` is not one of the
        supported modes.

    Examples
    --------
    >>> import groundinsight as gi
    >>> net = gi.create_network(name="demo", frequencies=[50.0])
    >>> bt = gi.BusType(name="b", system_type="s", voltage_level=20.0,
    ...                 impedance_formula="rho*0 + 1.0")
    >>> _ = gi.create_bus(name="B1", type=bt, network=net)
    >>> _ = gi.create_source(name="S1", bus="B1", values={50.0: 1000.0},
    ...                      network=net, r_to_x=0.1)
    >>> _ = gi.create_fault(name="F1", bus="B1", scalings={50.0: 1.0},
    ...                     network=net, t_k_s=0.5)
    >>> data = gi.resolve_fault_sc_characteristics(net, "F1")
    >>> round(data.kappa, 4)
    1.746
    """
    if aggregation not in ("weighted", "max"):
        raise ValueError(
            f"aggregation must be 'weighted' or 'max', got {aggregation!r}."
        )
    fault = network.faults.get(fault_name)
    if fault is None:
        raise ValueError(
            f"Unknown fault {fault_name!r}. Known faults: "
            f"{sorted(network.faults)}."
        )

    if frequency is None:
        frequency = next((float(x) for x in network.frequencies if x > 0), 50.0)
    frequency = float(frequency)
    scaling = fault.scalings.get(frequency, 1)

    weights: Dict[str, float] = {}
    kappas: Dict[str, float] = {}
    r_values: List[float] = []
    without_kappa: List[str] = []
    i_k_total = 0.0

    for source_name in _sources_feeding(network, fault_name):
        source = network.sources[source_name]
        weight = _source_injection_magnitude(source, frequency, scaling)
        kappa_i = _source_kappa(source)
        if kappa_i is None:
            if weight > 0:
                without_kappa.append(source_name)
            i_k_total += weight
            continue
        # A source that carries 60909 data but injects nothing at this
        # frequency would silently vanish from a purely current-weighted
        # mean; give it a nominal weight so its kappa still counts.
        if weight <= 0:
            weight = float(getattr(source, "i_k_a", None) or 0.0)
        weights[source_name] = weight
        kappas[source_name] = kappa_i
        r_to_x_i = getattr(source, "r_to_x", None)
        if r_to_x_i is not None:
            r_values.append(float(r_to_x_i))
        i_k_total += weight

    n_factor = float(getattr(fault, "n_factor", 1.0) or 1.0)
    t_k_s = getattr(fault, "t_k_s", None)
    t_k_s = None if t_k_s is None else float(t_k_s)

    if not kappas:
        if without_kappa:
            logger.debug(
                "Fault %r: none of the feeding sources (%s) carries IEC 60909 "
                "data; kappa stays unresolved.",
                fault_name,
                ", ".join(without_kappa),
            )
        return FaultShortCircuitData(
            fault_name=fault_name,
            frequency=frequency,
            t_k_s=t_k_s,
            n_factor=n_factor,
            i_k_a=i_k_total,
            aggregation=aggregation,
            sources_without_kappa=without_kappa,
        )

    kappa_list = list(kappas.values())
    homogeneous = (max(kappa_list) - min(kappa_list)) <= 1e-9

    if aggregation == "max" or homogeneous:
        kappa_eff = max(kappa_list)
    else:
        weight_sum = sum(weights.values())
        if weight_sum > 0:
            kappa_eff = (
                sum(kappas[name] * weights[name] for name in kappas) / weight_sum
            )
        else:  # no current information at all -> plain mean
            kappa_eff = sum(kappa_list) / len(kappa_list)

    if not homogeneous:
        logger.warning(
            "Fault %r is fed by sources with different R/X (kappa between "
            "%.4f and %.4f). The non-linear IEC 60909 factors are applied "
            "once, with the %s kappa = %.4f. Use the transient solver for an "
            "exact mixed-R/X result.",
            fault_name,
            min(kappa_list),
            max(kappa_list),
            aggregation,
            kappa_eff,
        )
    if without_kappa:
        logger.warning(
            "Fault %r: source(s) %s inject current at %.1f Hz but carry no "
            "IEC 60909 data; they are excluded from the kappa aggregation.",
            fault_name,
            ", ".join(without_kappa),
            frequency,
        )

    r_to_x_eff = None
    if r_values and (max(r_values) - min(r_values)) <= 1e-12:
        r_to_x_eff = r_values[0]

    m = None
    if t_k_s is not None and t_k_s > 0:
        m = iec60909_m(kappa_eff, frequency, t_k_s)

    return FaultShortCircuitData(
        fault_name=fault_name,
        frequency=frequency,
        kappa=kappa_eff,
        r_to_x=r_to_x_eff,
        t_k_s=t_k_s,
        n_factor=n_factor,
        m=m,
        i_k_a=i_k_total,
        i_p_a=peak_short_circuit_current(i_k_total, kappa_eff),
        aggregation=aggregation,
        homogeneous=homogeneous,
        sources=sorted(kappas),
        sources_without_kappa=without_kappa,
    )

run_reference_cases 

run_reference_cases(
    cases: Optional[List[ReferenceCase]] = None,
) -> pl.DataFrame

Run the closed-form reference cases and report the comparison.

Parameters:

Name Type Description Default
cases list of ReferenceCase

Defaults to :data:REFERENCE_CASES.

 None

Returns:

Type Description
DataFrame

One row per case: case, quantity, conditions, closed_form, model, rel_deviation, tolerance, agrees.

Examples:

>>> import groundinsight as gi
>>> gi.run_reference_cases()
Source code in src/groundinsight/analysis/reference.py 
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
def run_reference_cases(
    cases: Optional[List[ReferenceCase]] = None,
) -> pl.DataFrame:
    """
    Run the closed-form reference cases and report the comparison.

    Parameters
    ----------
    cases : list of ReferenceCase, optional
        Defaults to :data:`REFERENCE_CASES`.

    Returns
    -------
    pl.DataFrame
        One row per case: ``case``, ``quantity``, ``conditions``,
        ``closed_form``, ``model``, ``rel_deviation``, ``tolerance``,
        ``agrees``.

    Examples
    --------
    >>> import groundinsight as gi  # doctest: +SKIP
    >>> gi.run_reference_cases()  # doctest: +SKIP
    """
    rows = [case.evaluate() for case in (cases or REFERENCE_CASES)]
    frame = pl.DataFrame(rows)
    failed = frame.filter(~pl.col("agrees"))
    if failed.height:
        logger.warning(
            "%d reference case(s) missed their closed form: %s. Either the model "
            "is wrong or a stated boundary condition was not met -- the second "
            "is the more common finding, so read the 'conditions' column first.",
            failed.height,
            ", ".join(failed["case"].to_list()),
        )
    return frame

select_rho_f_from_catalog 

select_rho_f_from_catalog(
    network: Network,
    bus_names: List[str],
    u_limit: float,
    candidates: Dict[str, KVector],
    *,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None,
    sort_by: Literal[
        "max_epr_asc", "max_epr_desc", "name"
    ] = "max_epr_asc"
) -> pl.DataFrame

Pick admissible rho-f characteristics from a user-provided catalog.

Evaluates every candidate k in the catalog with :func:evaluate_max_epr_under_k and reports, per candidate, the maximum RMS EPR observed across the bus sweep, the per-bus EPRs and whether the candidate satisfies the limit u_limit. The catalog is typically a hand-curated list of soil scenarios (e.g. dry sand, wet clay, permafrost) or rho-f fits from previous measurements.

Parameters:

Name Type Description Default
network Network

The simulation network. Sources, branches and the buses to be swept must already be configured.

 required
bus_names List[str]

Buses whose impedance is parameterised by every candidate k and which are swept as fault locations.

 required
u_limit float

Upper bound on the RMS EPR (in volts) at any swept bus. Must be strictly positive.

 required
candidates Dict[str, KVector]

Mapping {name: (k1, k2, k3, k4, k5)} of the candidate rho-f characteristics. Names must be unique (Python dict guarantees that), tuples must each have length 5. May be empty -- the result is then an empty DataFrame with the documented schema.

 required
fault_scalings Optional[Dict[float, float]]

Frequency-resolved scalings for any fault that has to be created on the fly. See :func:evaluate_max_epr_under_k.

 None
run_fault_kwargs Optional[Dict[str, Any]]

Forwarded to :func:run_fault.

 None
sort_by Literal['max_epr_asc', 'max_epr_desc', 'name']

How to sort the result rows.

  • "max_epr_asc" (default): admissible candidates first, tightest-EPR-margin candidates at the top.
  • "max_epr_desc": largest EPR first; useful to inspect the worst-case candidates.
  • "name": lexicographic by candidate name.
 'max_epr_asc'

Returns:

Type Description
A Polars DataFrame with one row per candidate and the columns
- ``"name"`` (str)
- ``"k1", "k2", "k3", "k4", "k5"`` (float)
- ``"max_epr_rms_V"`` (float)

maximum RMS EPR across the bus sweep, in volts.
- ``"admissible"`` (bool)

True iff max_epr_rms_V <= u_limit.
- one ``"epr_<bus>_V"`` column per swept bus, holding the RMS

EPR at that bus when it is the active fault, in volts.

Raises:

Type Description
ValueError

If u_limit is not finite and strictly positive, bus_names is empty, or any candidate has a wrong-length k or refers to an unknown bus (the underlying helper validates this). A NaN u_limit is rejected explicitly: it would otherwise pass a plain positivity check and mark the whole catalog inadmissible.

Examples:

>>> import groundinsight as gi
>>> from groundinsight.models.core_models import BusType, BranchType
>>> from groundinsight.analysis import select_rho_f_from_catalog
>>> bt = BusType(name="BT", system_type="Grounded",
...              voltage_level=20.0,
...              impedance_formula="rho * 0.01 + 0*f")
>>> brt = BranchType(name="BRT", grounding_conductor=True,
...                  self_impedance_formula="(0.25 + I*0.6)*l",
...                  mutual_impedance_formula="(0.0 + I*0.6)*l")
>>> net = gi.create_network(name="N", frequencies=[50])
>>> _ = gi.create_bus(name="b0", type=bt, network=net)
>>> _ = gi.create_bus(name="b1", type=bt, network=net)
>>> _ = gi.create_branch(name="br", type=brt, from_bus="b0",
...                      to_bus="b1", length=1.0, network=net)
>>> _ = gi.create_source(name="src", bus="b0",
...                      values={50: 100.0}, network=net)
>>> catalog = {
...     "low":   (0.005, 0.0, 0.0, 0.0, 0.0),
...     "med":   (0.01,  0.0, 0.0, 0.0, 0.0),
...     "high":  (0.05,  0.0, 0.0, 0.0, 0.0),
... }
>>> df = select_rho_f_from_catalog(
...     net, ["b0", "b1"], u_limit=20.0, candidates=catalog,
... )
>>> set(df.columns) >= {"name", "k1", "max_epr_rms_V", "admissible"}
True
Source code in src/groundinsight/analysis/inverse_rho_f.py 
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
def select_rho_f_from_catalog(
    network: Network,
    bus_names: List[str],
    u_limit: float,
    candidates: Dict[str, KVector],
    *,
    fault_scalings: Optional[Dict[float, float]] = None,
    run_fault_kwargs: Optional[Dict[str, Any]] = None,
    sort_by: Literal["max_epr_asc", "max_epr_desc", "name"] = "max_epr_asc",
) -> pl.DataFrame:
    """Pick admissible rho-f characteristics from a user-provided catalog.

    Evaluates every candidate ``k`` in the catalog with
    :func:`evaluate_max_epr_under_k` and reports, per candidate, the
    maximum RMS EPR observed across the bus sweep, the per-bus EPRs and
    whether the candidate satisfies the limit ``u_limit``. The catalog is
    typically a hand-curated list of soil scenarios (e.g. dry sand, wet
    clay, permafrost) or rho-f fits from previous measurements.

    Parameters
    ----------
    network
        The simulation network. Sources, branches and the buses
        to be swept must already be configured.
    bus_names
        Buses whose impedance is parameterised by every
        candidate ``k`` and which are swept as fault locations.
    u_limit
        Upper bound on the RMS EPR (in volts) at any swept bus.
        Must be strictly positive.
    candidates
        Mapping ``{name: (k1, k2, k3, k4, k5)}`` of the
        candidate rho-f characteristics. Names must be unique
        (Python dict guarantees that), tuples must each have length
        5. May be empty -- the result is then an empty DataFrame
        with the documented schema.
    fault_scalings
        Frequency-resolved scalings for any fault that
        has to be created on the fly. See
        :func:`evaluate_max_epr_under_k`.
    run_fault_kwargs
        Forwarded to :func:`run_fault`.
    sort_by
        How to sort the result rows.

        - ``"max_epr_asc"`` (default): admissible candidates first,
        tightest-EPR-margin candidates at the top.
        - ``"max_epr_desc"``: largest EPR first; useful to inspect
        the worst-case candidates.
        - ``"name"``: lexicographic by candidate name.

    Returns
    -------
    A Polars DataFrame with one row per candidate and the columns

    - ``"name"`` (str)
    - ``"k1", "k2", "k3", "k4", "k5"`` (float)
    - ``"max_epr_rms_V"`` (float)
        maximum RMS EPR across the bus
        sweep, in volts.
    - ``"admissible"`` (bool)
        ``True`` iff
        ``max_epr_rms_V <= u_limit``.
    - one ``"epr_<bus>_V"`` column per swept bus, holding the RMS
        EPR at that bus when it is the active fault, in volts.

    Raises
    ------
    ValueError
        If ``u_limit`` is not finite and strictly positive,
        ``bus_names`` is empty, or any candidate has a wrong-length
        ``k`` or refers to an unknown bus (the underlying helper
        validates this). A NaN ``u_limit`` is rejected explicitly:
        it would otherwise pass a plain positivity check and mark
        the whole catalog inadmissible.

    Examples
    --------
        >>> import groundinsight as gi
        >>> from groundinsight.models.core_models import BusType, BranchType
        >>> from groundinsight.analysis import select_rho_f_from_catalog
        >>> bt = BusType(name="BT", system_type="Grounded",
        ...              voltage_level=20.0,
        ...              impedance_formula="rho * 0.01 + 0*f")
        >>> brt = BranchType(name="BRT", grounding_conductor=True,
        ...                  self_impedance_formula="(0.25 + I*0.6)*l",
        ...                  mutual_impedance_formula="(0.0 + I*0.6)*l")
        >>> net = gi.create_network(name="N", frequencies=[50])
        >>> _ = gi.create_bus(name="b0", type=bt, network=net)
        >>> _ = gi.create_bus(name="b1", type=bt, network=net)
        >>> _ = gi.create_branch(name="br", type=brt, from_bus="b0",
        ...                      to_bus="b1", length=1.0, network=net)
        >>> _ = gi.create_source(name="src", bus="b0",
        ...                      values={50: 100.0}, network=net)
        >>> catalog = {
        ...     "low":   (0.005, 0.0, 0.0, 0.0, 0.0),
        ...     "med":   (0.01,  0.0, 0.0, 0.0, 0.0),
        ...     "high":  (0.05,  0.0, 0.0, 0.0, 0.0),
        ... }
        >>> df = select_rho_f_from_catalog(
        ...     net, ["b0", "b1"], u_limit=20.0, candidates=catalog,
        ... )
        >>> set(df.columns) >= {"name", "k1", "max_epr_rms_V", "admissible"}
        True
    """
    # Same check as in the two searches, and for the same reason -- but the
    # consequence here is worse. This function does not return a scalar a
    # reader might sanity-check, it returns a table with an ``admissible``
    # column, and every entry in it is ``max_epr <= u_limit``. Against NaN
    # every one of those comparisons is False, so the table would report
    # that *no* soil model in the catalog is usable, next to an EPR column
    # that is correct and finite. Nothing in the output would point at the
    # limit as the broken part.
    validate_limit(u_limit, "u_limit")
    if not bus_names:
        raise ValueError("bus_names must not be empty.")

    # Empty catalog -> return an empty DataFrame with the documented schema.
    if not candidates:
        schema = {
            "name": pl.Utf8,
            "k1": pl.Float64, "k2": pl.Float64, "k3": pl.Float64,
            "k4": pl.Float64, "k5": pl.Float64,
            "max_epr_rms_V": pl.Float64,
            "admissible": pl.Boolean,
        }
        for b in bus_names:
            schema[f"epr_{b}_V"] = pl.Float64
        return pl.DataFrame(schema=schema)

    rows: List[Dict[str, Any]] = []
    for name, k in candidates.items():
        if len(k) != 5:
            raise ValueError(
                f"Candidate {name!r}: k must be a 5-tuple, got length {len(k)}."
            )
        eprs = evaluate_max_epr_under_k(
            network, bus_names, k=tuple(k),
            fault_scalings=fault_scalings,
            run_fault_kwargs=run_fault_kwargs,
        )
        max_epr = max(eprs.values()) if eprs else 0.0
        row: Dict[str, Any] = {
            "name": name,
            "k1": float(k[0]),
            "k2": float(k[1]),
            "k3": float(k[2]),
            "k4": float(k[3]),
            "k5": float(k[4]),
            "max_epr_rms_V": float(max_epr),
            "admissible": bool(max_epr <= u_limit),
        }
        for b in bus_names:
            row[f"epr_{b}_V"] = float(eprs[b])
        rows.append(row)

    df = pl.DataFrame(rows)

    if sort_by == "max_epr_asc":
        # Admissible first (True -> 1 sorts after False -> 0 by default,
        # so descending on ``admissible`` puts True on top), then
        # ascending EPR within each block.
        df = df.sort(
            by=["admissible", "max_epr_rms_V"], descending=[True, False]
        )
    elif sort_by == "max_epr_desc":
        df = df.sort(by="max_epr_rms_V", descending=True)
    elif sort_by == "name":
        df = df.sort(by="name")
    else:  # pragma: no cover -- guarded by Literal type, defensive
        raise ValueError(
            f"sort_by must be one of "
            f"'max_epr_asc' | 'max_epr_desc' | 'name', got {sort_by!r}."
        )

    return df

summarize 

summarize(
    frame: pl.DataFrame,
    value: str,
    *,
    by: Optional[Sequence[str]] = None,
    quantiles: Sequence[float] = DEFAULT_QUANTILES
) -> pl.DataFrame

Reduce one column to count, spread, quantiles and extremes.

Parameters:

Name Type Description Default
frame DataFrame

Long-format results, e.g. from :meth:~groundinsight.simulation.sweep.SweepResult.buses.

 required
value str

Numeric column to summarise.

 required
by sequence of str

Grouping columns. Without them the whole frame is one group.

 None
quantiles sequence of float

Quantiles in [0, 1]. Each becomes a column p05, p50, ...

 DEFAULT_QUANTILES

Returns:

Type Description
DataFrame

One row per group: the grouping columns, then n, n_null, mean, std, min, the quantiles, and max. Sorted by the grouping columns so repeated runs produce identical output.

Raises:

Type Description
ValueError

If a named column is missing, if the value column is not numeric, or if a quantile lies outside [0, 1].

Examples:

>>> summarize(study.buses(), "EPR_V", by=["bus_name"])
Source code in src/groundinsight/analysis/statistics.py 
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def summarize(
    frame: pl.DataFrame,
    value: str,
    *,
    by: Optional[Sequence[str]] = None,
    quantiles: Sequence[float] = DEFAULT_QUANTILES,
) -> pl.DataFrame:
    """
    Reduce one column to count, spread, quantiles and extremes.

    Parameters
    ----------
    frame : pl.DataFrame
        Long-format results, e.g. from
        :meth:`~groundinsight.simulation.sweep.SweepResult.buses`.
    value : str
        Numeric column to summarise.
    by : sequence of str, optional
        Grouping columns. Without them the whole frame is one group.
    quantiles : sequence of float, optional
        Quantiles in ``[0, 1]``. Each becomes a column ``p05``, ``p50``, ...

    Returns
    -------
    pl.DataFrame
        One row per group: the grouping columns, then ``n``, ``n_null``,
        ``mean``, ``std``, ``min``, the quantiles, and ``max``. Sorted by the
        grouping columns so repeated runs produce identical output.

    Raises
    ------
    ValueError
        If a named column is missing, if the value column is not numeric, or if
        a quantile lies outside ``[0, 1]``.

    Examples
    --------
    >>> summarize(study.buses(), "EPR_V", by=["bus_name"])  # doctest: +SKIP
    """
    missing = [c for c in [value, *(by or [])] if c not in frame.columns]
    if missing:
        raise ValueError(
            f"Column(s) {missing} are not in the frame. Available: "
            f"{frame.columns}."
        )
    if not frame.schema[value].is_numeric():
        raise ValueError(
            f"Column '{value}' has dtype {frame.schema[value]}, which cannot be "
            f"summarised numerically. Note that res_buses() reports "
            f"'frequency_Hz' as a string because it carries the 'RMS' marker "
            f"row -- filter or cast it before grouping on it."
        )
    bad = [q for q in quantiles if not 0.0 <= q <= 1.0]
    if bad:
        raise ValueError(f"Quantile(s) {bad} lie outside [0, 1].")

    aggregations = [
        pl.col(value).count().alias("n"),
        pl.col(value).null_count().alias("n_null"),
        pl.col(value).mean().alias("mean"),
        pl.col(value).std().alias("std"),
        pl.col(value).min().alias("min"),
    ]
    aggregations += [
        pl.col(value).quantile(q, interpolation="linear").alias(_quantile_name(q))
        for q in quantiles
    ]
    aggregations.append(pl.col(value).max().alias("max"))

    if by:
        return frame.group_by(list(by)).agg(aggregations).sort(list(by))
    return frame.select(aggregations)

thermal_equivalent_current 

thermal_equivalent_current(
    i_k: float, m: float, n: float = 1.0
) -> float

Thermally equivalent short-time current I_th = I_k'' * sqrt(m + n).

Parameters:

Name Type Description Default
i_k float

RMS short-circuit current I_k'' in amperes. Must be non-negative.

 required
m float

DC heat-effect factor from :func:iec60909_m. Must be non-negative.

 required
n float

AC-decay heat factor in (0, 1]. 1.0 for the far-from-generator faults typical of distribution grounding studies.

 1.0

Returns:

Type Description
float

The thermally equivalent current in amperes.

Raises:

Type Description
ValueError

If i_k or m is negative, or n is not in (0, 1].
Notes

The bounds are enforced rather than clipped because violating them fails silently downwards: a negative m, or an n above one supplied in place of a correct smaller value, changes sqrt(m + n) without any visible symptom, and a too-small result under-estimates the thermal stress -- the unsafe direction for a limit check. A caller that hands in such a value has a bug upstream and should hear about it.

Examples:

>>> round(thermal_equivalent_current(1000.0, 0.0, 1.0), 6)
1000.0
Source code in src/groundinsight/analysis/shortcircuit.py 
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
def thermal_equivalent_current(i_k: float, m: float, n: float = 1.0) -> float:
    """
    Thermally equivalent short-time current ``I_th = I_k'' * sqrt(m + n)``.

    Parameters
    ----------
    i_k : float
        RMS short-circuit current ``I_k''`` in amperes. Must be
        non-negative.
    m : float
        DC heat-effect factor from :func:`iec60909_m`. Must be
        non-negative.
    n : float
        AC-decay heat factor in ``(0, 1]``. ``1.0`` for the
        far-from-generator faults typical of distribution grounding
        studies.

    Returns
    -------
    float
        The thermally equivalent current in amperes.

    Raises
    ------
    ValueError
        If ``i_k`` or ``m`` is negative, or ``n`` is not in ``(0, 1]``.

    Notes
    -----
    The bounds are enforced rather than clipped because violating them
    fails *silently downwards*: a negative ``m``, or an ``n`` above one
    supplied in place of a correct smaller value, changes ``sqrt(m + n)``
    without any visible symptom, and a too-small result under-estimates the
    thermal stress -- the unsafe direction for a limit check. A caller that
    hands in such a value has a bug upstream and should hear about it.

    Examples
    --------
    >>> round(thermal_equivalent_current(1000.0, 0.0, 1.0), 6)
    1000.0
    """
    if i_k < 0:
        raise ValueError(f"i_k must be non-negative, got {i_k!r}.")
    if m < 0:
        raise ValueError(f"m must be non-negative, got {m!r}.")
    if not 0.0 < n <= 1.0:
        raise ValueError(f"n must lie in (0, 1], got {n!r}.")
    return i_k * math.sqrt(m + n)