# Учебная редакция: назначение, единицы и этапы алгоритма пояснены по-русски.
# Исходная программа остаётся самостоятельным свидетельством прежней проверки.
# Развиваем характеристику II02-2.py до взаимной системы двух обмоток.
# Функция расчёта FEMM должна решать новое поле при каждом вызове и неизменном угле.
def linkage(i):
    """Вернуть взаимные аналитические потокосцепления двух обмоток, Вб, по токам, А.
    
    Одинаковые перекрёстные коэффициенты 0,02 Гн обеспечивают взаимность.
    Кубические члены представляют учебную характеристику, не опытную сталь."""
    return [.1*i[0]+.02*i[1]+.002*i[0]**3, .02*i[0]+.08*i[1]+.001*i[1]**3]


def jacobian(solve, currents, step):
    """Вычислить матрицу dΨ/dI центральными разностями, Гн.
    
    Строки соответствуют потокосцеплениям Ψ1, Ψ2, столбцы — токам I1, I2.
    currents содержит два тока в А; step — положительное возмущение в А.
    При замене solve на FEMM каждый вызов должен решать новое поле
    при одной и той же геометрии и неизменном угле ротора."""
    if step <= 0:
        raise ValueError("Positive current step required")
    columns = []
    for k in range(2):
        # Создаём независимые копии токов: меняем одну координату и сохраняем вторую.
        plus, minus = list(currents), list(currents)
        plus[k] += step
        minus[k] -= step
        a, b = solve(plus), solve(minus)
        # Расстояние между двумя пробными токами равно 2·step.
        columns.append([(a[j]-b[j])/(2*step) for j in range(2)])
    return [[columns[k][j] for k in range(2)] for j in range(2)]


def run_checks():
    """Сопоставить два шага с точной производной, взаимностью и положительным определителем."""
    one, two = jacobian(linkage, [2., 1.], .1), jacobian(linkage, [2., 1.], .05)
    # Точная производная задана независимо от конечно-разностной реализации.
    exact = [[.124, .02], [.02, .083]]
    err = lambda m: max(abs(m[j][k]-exact[j][k]) for j in range(2) for k in range(2))
    assert abs(one[0][0]-.12402) < 1e-12 and abs(two[1][1]-.0830025) < 1e-12
    # Уменьшение шага вдвое должно уменьшить ошибку второго порядка в четыре раза.
    assert abs(err(one)/err(two)-4) < 1e-6 and abs(two[0][1]-two[1][0]) < 1e-12
    assert two[0][0]*two[1][1]-two[0][1]*two[1][0] > 0
    return {"units": "H", "order": "rows Psi1,Psi2; columns I1,I2", "step_0_1": one, "step_0_05": two, "exact": exact}