import math
# Постоянная C = μ₀·N²·S в согласованных единицах СИ; L = C/g.
# Здесь проверяется производная аналитической коэнергии, а не режим проницаемости FEMM.
constant = 4*math.pi*1e-7*200**2*4e-4
current_a, gap_m = 1.0, 0.001
# Коэнергия, Дж, как функция зазора, м; ток во всех состояниях сохраняется.
def coenergy_j(gap):
    return constant*current_a**2/(2*gap)
# Независимая производная даёт силу, Н. Отрицательный знак означает уменьшение зазора.
exact_force_n = -constant*current_a**2/(2*gap_m**2)
# Сравниваем два механических шага. Для переноса геометрии FEMM шаг переводят в мм.
for step_m in [0.00005, 0.000025]:
    force_n = (coenergy_j(gap_m+step_m)-coenergy_j(gap_m-step_m))/(2*step_m)
    # Порядок вывода: шаг, м; разностная сила, Н; точная аналитическая сила, Н.
    print(step_m, force_n, exact_force_n)