"""Учебное исследование магнитного поля во встроенном Python FEMM 5.

Геометрия — мм, результаты — СИ, гармонические токи — амплитуды.
Каждый запуск создаёт calculated-NN; прежние модели и решения сохраняются.
Каждая модель сохраняется, закрывается, повторно открывается и решается.
Группы: 0 — воздух, 1 — сердечник, 2 — исследуемое тело, 3/4 — обмотки.
Условная B(H) взята из учебного комплекта курса, это не марка стали.
Источники методики и точные допущения приведены в Порядок-работы.md.
"""
import os
import sys
import math
import cmath
import json

MU0 = 4e-7 * math.pi
BH = [(0, 0), (.3, 100), (.6, 200), (1, 350), (1.3, 700),
      (1.5, 1800), (1.65, 5000), (1.8, 15000), (2, 60000), (2.2, 200000)]
source = globals().get('__file__', '') or (sys.argv[0] if sys.argv else '')
HERE = os.path.dirname(os.path.abspath(source)) if source else os.getcwd()
counter = 1
while os.path.exists(HERE + '/calculated-%02d' % counter):
    counter += 1
OUT = HERE + '/calculated-%02d' % counter
os.mkdir(OUT)
ROWS = []
CHECKS = []
NODES = {}
EDGES = {}


def real(value):
    return complex(value).real


def parts(value):
    return [complex(value).real, complex(value).imag]


def rel(a, b, floor=1e-12):
    return abs(a - b) / max(abs(a), abs(b), floor)


def check(name, error, limit):
    CHECKS.append({'name': name, 'error': error, 'limit': limit,
                   'passed': math.isfinite(error) and error <= limit})


def edge(a, b, size=0, boundary='<None>'):
    for point in (a, b):
        if point not in NODES:
            mi_addnode(*point)
            NODES[point] = True
    key = tuple(sorted((a, b)))
    if key not in EDGES:
        mi_addsegment(a[0], a[1], b[0], b[1])
        EDGES[key] = True
    mi_selectsegment((a[0] + b[0]) / 2, (a[1] + b[1]) / 2)
    mi_setsegmentprop(boundary, size, 0 if size else 1, 0, 0)
    mi_clearselected()


def polygon(points, size=0):
    for k in range(len(points)):
        edge(points[k], points[(k + 1) % len(points)], size)


def rect(x0, y0, x1, y1, size=0):
    polygon([(x0, y0), (x1, y0), (x1, y1), (x0, y1)], size)


def label(x, y, material, size, group=0, circuit='<None>', turns=0, angle=0):
    mi_addblocklabel(x, y)
    mi_selectlabel(x, y)
    mi_setblockprop(material, 0, size, circuit, angle, group, turns)
    mi_clearselected()


def begin(axis=True, frequency=0, depth=20):
    global NODES, EDGES
    NODES, EDGES = {}, {}
    newdocument(0)
    mi_probdef(frequency, 'millimeters', 'axi' if axis else 'planar', 1e-10, depth, 25)
    mi_addmaterial('Воздух', 1, 1, 0, 0, 0)
    mi_addmaterial('Обмотка без потерь', 1, 1, 0, 0, 0)
    mi_addmaterial('Сталь линейная учебная', 1000, 1000, 0, 0, 0)
    mi_addmaterial('Сталь с учебной кривой насыщения', 1, 1, 0, 0, 0)
    for b, h in BH:
        mi_addbhpoint('Сталь с учебной кривой насыщения', b, h)
    mi_addboundprop('Внешняя граница A=0', 0, 0, 0, 0, 0, 0, 0, 0, 0)


def exterior(axis=True, scale=1, mesh=4, cuts=()):
    bound = 'Внешняя граница A=0'
    radius, bottom, top = 120 * scale, -120 * scale, 120 * scale
    if axis:
        levels = sorted(set([bottom, top] + list(cuts)))
        for a, b in zip(levels[:-1], levels[1:]):
            edge((0, a), (0, b), 0, bound)
        for a, b in [((0, bottom), (radius, bottom)),
                     ((radius, bottom), (radius, top)), ((radius, top), (0, top))]:
            edge(a, b, 0, bound)
        label(radius * .8, 0, 'Воздух', mesh)
    else:
        for a, b in [((-radius, bottom), (radius, bottom)),
                     ((radius, bottom), (radius, top)), ((radius, top), (-radius, top)),
                     ((-radius, top), (-radius, bottom))]:
            edge(a, b, 0, bound)
        label(radius * .8, 0, 'Воздух', mesh)


def coil(r0, r1, z0, z1, name, turns, current, mesh, group=3):
    mi_addcircprop(name, current, 1)
    rect(r0, z0, r1, z1, mesh)
    label((r0 + r1) / 2, (z0 + z1) / 2, 'Обмотка без потерь', mesh, group, name, turns)


def solve(tag, circuits=(), force_group=0, probe=(4, 0), **parameters):
    path = OUT + '/' + tag + '.fem'
    mi_setcomment('Учебная адаптация магнитных полей. Группы: 0 воздух, 1 сердечник, '
                  '2 исследуемое тело, 3 и 4 обмотки. Размеры мм, результаты СИ. ' + tag)
    mi_saveas(path)
    mi_close()
    opendocument(path)
    mi_analyze(1)
    mi_loadsolution()
    mo_clearblock()
    mo_groupselectblock()
    row = {'tag': tag, 'model': tag + '.fem', 'solution': tag + '.ans',
           'elements': mo_numelements(), 'energy_J': real(mo_blockintegral(2)),
           'coenergy_J': real(mo_blockintegral(17)), 'loss_W': real(mo_blockintegral(4))}
    mo_clearblock()
    if force_group:
        mo_groupselectblock(force_group)
        row['force_N'] = real(mo_blockintegral(19))
        row['lorentz_N'] = real(mo_blockintegral(12))
        mo_clearblock()
    row['circuits'] = []
    for name in circuits:
        current, voltage, psi = mo_getcircuitproperties(name)
        row['circuits'].append({'name': name, 'current_A': parts(current),
                                'voltage_V': parts(voltage), 'psi_Wb': parts(psi)})
    if probe:
        values = mo_getpointvalues(*probe)
        row.update(Bx_T=real(values[1]), By_T=real(values[2]),
                   Hx_Am=real(values[5]), Hy_Am=real(values[6]))
        row['Mx_Am'] = row['Bx_T'] / MU0 - row['Hx_Am']
        row['My_Am'] = row['By_T'] / MU0 - row['Hy_Am']
    row.update(parameters)
    mo_close()
    mi_close()
    ROWS.append(row)
    with open(OUT + '/progress.json', 'w') as handle:
        json.dump({'solved': len(ROWS), 'last': tag}, handle)
    print('SOLVED', len(ROWS), tag, row['elements'])
    return row


def psi(row, index=0):
    return complex(*row['circuits'][index]['psi_Wb'])


def energy_check(row):
    work = sum(real(complex(*c['current_A']).conjugate() * complex(*c['psi_Wb']))
               for c in row['circuits'])
    check('Энергия и коэнергия ' + row['tag'],
          rel(row['energy_J'] + row['coenergy_J'], work, 1e-9), .01)


def finish(identifier, derived, settings):
    result = {'id': identifier, 'passed': all(c['passed'] for c in CHECKS),
              'settings': settings, 'records': ROWS, 'derived': derived, 'checks': CHECKS}
    with open(OUT + '/result.json', 'w') as handle:
        json.dump(result, handle)
    print('RESULT', identifier, result['passed'], len(ROWS), len(CHECKS))
    assert result['passed'], 'Не пройдены проверки; изучите result.json.'


# F07a: одна пластина, одинаковая средняя амплитуда B при всех t и f.
def sheet(thickness, frequency, cells=8):
    begin(False, frequency, 1000)
    mi_addmaterial('Проводящая пластина', 100, 100, 0, 0, 6)
    a = .2 * thickness * .001 / 2
    mi_addboundprop('Верхняя поверхность', a, 0, 0, 0, 0, 0, 0, 0, 0)
    mi_addboundprop('Нижняя поверхность', -a, 0, 0, 0, 0, 0, 0, 0, 0)
    rect(-10, -thickness / 2, 10, thickness / 2)
    for y, name in ((thickness / 2, 'Верхняя поверхность'),
                    (-thickness / 2, 'Нижняя поверхность')):
        mi_selectsegment(0, y)
        mi_setsegmentprop(name, 1, 0, 0, 0)
        mi_clearselected()
    label(0, 0, 'Проводящая пластина', thickness / cells, 1)
    row = solve('F07-t%g-f%g-n%d' % (thickness, frequency, cells),
                probe=(0, 0), thickness_mm=thickness, frequency_Hz=frequency,
                cells=cells, family='sheet')
    volume = 20 * thickness * 1e-6
    row['loss_density_W_m3'] = row['loss_W'] / volume
    row['thin_sheet_W_m3'] = math.pi ** 2 / 6 * 6e6 * frequency ** 2 * .2 ** 2 * (thickness * .001) ** 2
    skin = math.sqrt(2 / (2 * math.pi * frequency * MU0 * 100 * 6e6))
    z = (1 + 1j) * (thickness * .001 / 2) / skin
    tanh = (cmath.exp(2 * z) - 1) / (cmath.exp(2 * z) + 1)
    inverse_mu = z / (MU0 * 100 * tanh)
    row['exact_slab_W_m3'] = math.pi * frequency * .2 ** 2 * inverse_mu.imag
    row['skin_depth_mm'] = skin * 1000
    check('Аналитический проводящий слой ' + row['tag'],
          rel(row['loss_density_W_m3'], row['exact_slab_W_m3']), .03)
    if thickness * .001 / skin < .2:
        check('Тонколистовой предел ' + row['tag'],
              rel(row['loss_density_W_m3'], row['thin_sheet_W_m3']), .03)
    return row


sheet_series = []
for thickness in (.1, .2, .4):
    for frequency in (10, 50, 200, 1000):
        sheet_series.append(sheet(thickness, frequency)['tag'])
fine = sheet(.4, 1000, 12)
base = [r for r in ROWS if r['tag'] == sheet_series[-1]][0]
check('Сетка потерь пластины 8/12 элементов по толщине', rel(base['loss_W'], fine['loss_W']), .03)


# F07b: магнитное заполнение при sigma=0, штатная модель пакета FEMM.
# Плоский пакет: листы лежат в плоскости xy, нормаль к листам z.
def packing(fill, analytical=False, mesh=.4, domain=1):
    begin(False)
    name = 'Шихтованный пакет в плоскости модели'
    if analytical:
        mu_parallel = fill * 1000 + 1 - fill
        mi_addmaterial(name, mu_parallel, mu_parallel, 0, 0, 0)
    else:
        mi_addmaterial(name, 1000, 1000, 0, 0, 0, .35, 0, fill, 0)
    rect(-30, -30, 30, 30, mesh)
    rect(-20, -20, 20, 20, mesh)
    label(0, 25, name, mesh, 1)
    mi_addcircprop('Обмотка возбуждения', 1, 1)
    for x0, x1, turns in ((-40, -31, 200), (-19, -10, -200)):
        rect(x0, -10, x1, 10, mesh)
        label((x0 + x1) / 2, 0, 'Обмотка без потерь', mesh, 3, 'Обмотка возбуждения', turns)
    label(0, 0, 'Воздух', mesh * 2)
    exterior(False, domain, 5)
    row = solve('F07-fill%g-a%d-h%g-d%g' % (fill, analytical, mesh, domain),
                ['Обмотка возбуждения'], probe=(0, 25), fill=fill, analytical=analytical,
                mesh_mm=mesh, domain=domain, family='packing')
    row['L_H'] = real(psi(row))
    energy_check(row)
    return row


fill_series = []
for fill in (.7, .85, 1):
    native = packing(fill)
    equivalent = packing(fill, True)
    check('Штатная модель пакета и аналитический тензор mu, f=%g' % fill,
          rel(native['L_H'], equivalent['L_H']), .01)
    fill_series.append(native['tag'])
fine = packing(.85, mesh=.25)
large = packing(.85, mesh=.25, domain=1.5)
base = [r for r in ROWS if r['tag'] == fill_series[1]][0]
check('Сетка коэффициента заполнения', rel(base['L_H'], fine['L_H']), .01)
check('Внешняя область коэффициента заполнения', rel(fine['L_H'], large['L_H']), .01)
finish('lamination', {'sheet_series': sheet_series, 'fill_series': fill_series},
       {'sheet_B_average_peak_T': .2, 'sheet_mu_r': 100, 'sheet_sigma_S_m': 6e6,
        'sheet_width_mm': 20, 'sheet_depth_mm': 1000,
        'packing_sigma_S_m': 0, 'packing_mu_steel': 1000, 'packing_LamType': 0,
        'packing_N_I': [200, 1], 'packing_orientation': 'Листы в плоскости xy, нормаль z',
        'packing_scope': 'Плоский квадратный сердечник 60/40 мм, Depth=20 мм; габариты и NI постоянны.',
        'effective_mu_formula': 'mu_parallel=f*mu+(1-f); 1/mu_normal=f/mu+(1-f)'})