"""Учебное исследование магнитного поля во встроенном 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.'


# F03: та же квадратная геометрия, что в занятиях 8.7–8.11 курса.
def core_state(currents, nonlinear=True, mesh=.55, domain=1, suffix=''):
    begin(False)
    # Разбиение на четыре непересекающихся прямоугольных стержня.
    polygon([(-30, -30), (30, -30), (30, 30), (-30, 30)], mesh)
    polygon([(-20, -20), (20, -20), (20, 20), (-20, 20)], mesh)
    steel = 'Сталь с учебной кривой насыщения' if nonlinear else 'Сталь линейная учебная'
    label(0, 25, steel, mesh, 1)
    # Стороны обмоток отделены от стали воздушным технологическим промежутком.
    # Размеры 60/40 мм, Depth=20 мм и витки 100/50 сохранены из основной модели.
    for j, entries in enumerate((((-40, -31, 100), (-19, -10, -100)),
                                 ((10, 19, -50), (31, 40, 50)))):
        name = 'Обмотка %d' % (j + 1)
        mi_addcircprop(name, currents[j], 1)
        for x0, x1, turns in entries:
            rect(x0, -10, x1, 10, mesh)
            label((x0 + x1) / 2, 0, 'Обмотка без потерь', mesh, 3 + j, name, turns)
    label(0, 0, 'Воздух', mesh * 2)
    exterior(False, domain, 6)
    row = solve('F03-n%d-i%g-%g-h%g-d%g%s' %
                (nonlinear, currents[0], currents[1], mesh, domain, suffix),
                ['Обмотка 1', 'Обмотка 2'], probe=(0, 25),
                currents_A=list(currents), nonlinear=nonlinear, mesh_mm=mesh, domain=domain)
    energy_check(row)
    return row


curves = []
for nonlinear in (False, True):
    for j in range(2):
        for current in (1, 5, 20, 50):
            currents = [0, 0]
            currents[j] = current
            row = core_state(currents, nonlinear)
            curves.append({'nonlinear': nonlinear, 'excited_winding': j + 1,
                           'current_A': current, 'psi1_Wb': real(psi(row)),
                           'psi2_Wb': real(psi(row, 1)),
                           'B_T': math.sqrt(row['Bx_T'] ** 2 + row['By_T'] ** 2),
                           'M_Am': math.sqrt(row['Mx_Am'] ** 2 + row['My_Am'] ** 2), 'tag': row['tag']})
local = []
for step in (.2, .1):
    columns = []
    for j in range(2):
        low, high = [20, 10], [20, 10]
        low[j] -= step
        high[j] += step
        a, b = core_state(low, suffix='-jac'), core_state(high, suffix='-jac')
        columns.append([real(psi(b, k) - psi(a, k)) / (2 * step) for k in range(2)])
    check('Симметрия локальной матрицы при шаге %g' % step,
          rel(columns[0][1], columns[1][0]), .01)
    local.append(columns)
for j in range(2):
    for k in range(2):
        check('Шаг Якобиана %d%d' % (j, k), rel(local[0][j][k], local[1][j][k]), .01)
base = core_state([20, 10])
fine = core_state([20, 10], mesh=.35)
large = core_state([20, 10], mesh=.35, domain=1.5)
for k in range(2):
    check('Сетка потокосцепления %d' % k, rel(psi(base, k), psi(fine, k)), .01)
    check('Внешняя область %d' % k, rel(psi(fine, k), psi(large, k)), .01)
for j in range(2):
    low, high = [20, 10], [20, 10]
    low[j] -= .1
    high[j] += .1
    a, b = core_state(low, mesh=.35, suffix='-jac-fine'), core_state(high, mesh=.35, suffix='-jac-fine')
    for k in range(2):
        value = real(psi(b, k) - psi(a, k)) / .2
        check('Сетка локального Якобиана %d%d' % (j, k), rel(local[-1][j][k], value), .01)
field_map = []
opendocument(OUT + '/' + fine['model'])
mi_loadsolution()
mo_smooth('off')
for x in range(-28, 29, 4):
    for y in range(-28, 29, 4):
        if abs(x) > 20 or abs(y) > 20:
            value = mo_getpointvalues(x, y)
            bx, by, hx, hy = [real(value[k]) for k in (1, 2, 5, 6)]
            field_map.append({'x_mm': x, 'y_mm': y, 'Bx_T': bx, 'By_T': by,
                              'Hx_Am': hx, 'Hy_Am': hy,
                              'Mx_Am': bx / MU0 - hx, 'My_Am': by / MU0 - hy})
mo_close()
mi_close()
finish('core-matrix', {'curves': curves, 'local_matrix_columns_H': local[-1],
                       'field_map': field_map, 'field_map_model': fine['model']},
       {'outer_square_mm': 60, 'inner_square_mm': 40, 'depth_mm': 20,
        'turns': [100, 50], 'local_currents_A': [20, 10], 'BH_T_Am': BH,
        'winding_clearance_mm': 1,
        'geometry_origin': 'Учебная квадратная модель основной части курса; добавлен явный зазор обмотка-сталь 1 мм.'})