"""Энергия за цикл синтетической динамической B-H-петли.

Среда: встроенный MicroPython FEMM 5.
Вход: амплитуды Bpeak [Т], Hsin и Hcos [А/м], f [Гц], rho [кг/м³];
Ns и Np — числа витков, area [м²] и length [м] — эффективные размеры.
Выход: synthetic-loop.csv и results.json с энергией [Дж/м³] и потерями [Вт/кг].
Это авторский математический сигнал, а не марка стали FEMM и не измерение.
Предполагаются однородные B и H, одинаковый поток витков и замкнутый период.
OUTPUT должен отсутствовать; для повторного запуска выберите новое имя.
"""

import json
import math
import os

OUTPUT = "measurement-bh-run"
# Существующий каталог не заменяется: сохраняем предыдущую запись опыта.
os.mkdir(OUTPUT)

Bpeak = 1.2
Hsin = 80.0
Hcos = 30.0
f = 50.0
rho = 7800.0
Ns = 100
area = 1e-4
Np = 100
length = 0.25
count = 4000  # Число интервалов; последняя точка замыкает цикл.
omega = 2 * math.pi * f
dt = 1 / (f * count)
emf_peak = Ns * area * omega * Bpeak

# Интеграл H*dB вычисляем потоково: нужны только предыдущие B и H.
# Положительное Hcos задаёт положительную энергию при выбранном обходе.
# Модуль интеграла намеренно не берём: он скрыл бы ошибку полярности.
area_loop = 0.0
previous_b = 0.0
previous_h = 0.0
with open(OUTPUT + "/synthetic-loop.csv", "w") as stream:
    stream.write(
        "time_s,synthetic_B_T,synthetic_H_A_per_m,synthetic_current_A,synthetic_voltage_V\n"
    )
    for k in range(count + 1):
        angle = 2 * math.pi * k / count
        sin_angle = math.sin(angle)
        cos_angle = math.cos(angle)
        b = Bpeak * sin_angle
        h = Hsin * sin_angle + Hcos * cos_angle

        # Закон полного тока H*l = Np*i и закон Фарадея e = -Ns*S*dB/dt.
        # Оба выражения используют явно указанные идеализации магнитной цепи.
        current = h * length / Np
        voltage = -emf_peak * cos_angle
        if k > 0:
            # Трапеция на плоскости B-H: Дж/м³ = (А/м)*Т.
            area_loop += 0.5 * (previous_h + h) * (b - previous_b)

        # CSV сохраняем построчно без массива из count + 1 строк.
        row = (k * dt, b, h, current, voltage)
        stream.write(",".join(str(value) for value in row) + "\n")
        previous_b = b
        previous_h = h

# Точный интеграл: Hsin даёт ноль, интеграл cos² за цикл равен pi.
# Умножение на частоту даёт Вт/м³, деление на плотность — Вт/кг.
exact = math.pi * Bpeak * Hcos
relative_error = abs(area_loop / exact - 1)
assert relative_error < 1e-6, "Квадратура не достигла заданной точности"
result = {
    "origin": "synthetic_analytical_ellipse_not_FEMM_or_measured_steel",
    "loop_energy_J_per_m3": area_loop,
    "exact_J_per_m3": exact,
    "power_density_W_per_m3": f * area_loop,
    "specific_loss_W_per_kg": f * area_loop / rho,
    "relative_quadrature_error": relative_error,
}
with open(OUTPUT + "/results.json", "w") as stream:
    json.dump(result, stream)
print(json.dumps(result))