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


# F01: неизменная обмотка, меняется только материал центрального стержня.
MATERIALS = ['Воздух', 'Сталь линейная учебная', 'Сталь с учебной кривой насыщения']


def single(material, current, mesh=.5, domain=1, suffix=''):
    begin()
    rect(0, -30, 8, 30, mesh)
    label(4, 0, MATERIALS[material], mesh, 1)
    coil(10, 16, -10, 10, 'Обмотка возбуждения', 200, current, mesh)
    exterior(scale=domain, mesh=5, cuts=(-30, 30))
    row = solve('F01-%d-I%g-h%g-d%g%s' % (material, current, mesh, domain, suffix),
                ['Обмотка возбуждения'], material=material, current_A=current,
                mesh_mm=mesh, domain=domain, family='single')
    row['L_sec_H'] = real(psi(row)) / current
    energy_check(row)
    return row


single_results = []
for material in range(3):
    series = [single(material, current) for current in (1, 5, 10, 20, 40)]
    derivatives = []
    for step in (.4, .2):
        low = single(material, 40 - step, suffix='-minus')
        high = single(material, 40 + step, suffix='-plus')
        derivatives.append(real(psi(high) - psi(low)) / (2 * step))
    check('Шаг производной собственной L, материал %d' % material,
          rel(*derivatives), .01)
    coarse = single(material, 40, .8)
    finer = single(material, 40, .32)
    enlarged = single(material, 40, .32, 1.5)
    check('Сетка собственной L, материал %d' % material,
          rel(finer['L_sec_H'], series[-1]['L_sec_H']), .01)
    check('Граница собственной L, материал %d' % material,
          rel(finer['L_sec_H'], enlarged['L_sec_H']), .01)
    check('Рост числа элементов, материал %d' % material,
          0 if finer['elements'] > coarse['elements'] else 1, 0)
    single_results.append({'material': MATERIALS[material], 'L_diff_40A_H': derivatives[-1],
                           'L_sec_40A_H': series[-1]['L_sec_H'],
                           'series': [r['tag'] for r in series]})
check('Насыщение уменьшает касательную относительно секущей',
      0 if single_results[2]['L_diff_40A_H'] < single_results[2]['L_sec_40A_H'] else 1, 0)


# F02: расстояние задано между обращёнными торцами; стержень неизменен.
def pair(material, gap, currents, mesh=.5, domain=1, suffix=''):
    begin()
    rect(0, -55, 8, 55, mesh)
    label(4, 0, MATERIALS[material], mesh, 1)
    for j, sign in enumerate((-1, 1)):
        center = sign * (10 + gap / 2)
        coil(10, 16, center - 10, center + 10, 'Обмотка %d' % (j + 1),
             (200, 100)[j], currents[j], mesh, 3 + j)
    exterior(scale=3 * domain, mesh=5, cuts=(-55, 55))
    row = solve('F02-%d-g%g-i%g-%g-h%g-d%g%s' %
                (material, gap, currents[0], currents[1], mesh, domain, suffix),
                ['Обмотка 1', 'Обмотка 2'], material=material, gap_mm=gap,
                currents_A=list(currents), mesh_mm=mesh, domain=domain, family='pair')
    energy_check(row)
    return row


matrices = []
for material in (0, 1):
    for gap in (5, 15, 30):
        a = pair(material, gap, (1, 0))
        b = pair(material, gap, (0, 1))
        l11, l21 = real(psi(a)), real(psi(a, 1))
        l12, l22 = real(psi(b)), real(psi(b, 1))
        check('Взаимность матрицы %d/%d' % (material, gap), rel(l12, l21), .01)
        k = (l12 + l21) / (2 * math.sqrt(l11 * l22))
        check('Положительная энергия и допустимый k %d/%d' % (material, gap),
              0 if l11 > 0 and l22 > 0 and abs(k) < 1 else 1, 0)
        matrices.append({'material': MATERIALS[material], 'gap_mm': gap,
                         'L11_H': l11, 'L22_H': l22, 'M12_H': l12, 'M21_H': l21, 'k': k})
    fine = pair(material, 15, (1, 0), .32)
    large = pair(material, 15, (1, 0), .32, 1.5)
    check('Сетка взаимной L, материал %d' % material,
          rel(real(psi(fine, 1)), matrices[-2]['M21_H']), .01)
    check('Граница взаимной L, материал %d' % material,
          rel(psi(fine, 1), psi(large, 1)), .01)

# При насыщении вычисляется Якобиан в одной рабочей точке (40, 20) А.
local = []
for step in (.4, .2):
    columns = []
    for j in range(2):
        low, high = [40, 20], [40, 20]
        low[j] -= step
        high[j] += step
        a = pair(2, 15, low, suffix='-jac')
        b = pair(2, 15, 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)
local_refinements = []
for domain in (1, 1.5):
    columns = []
    for j in range(2):
        low, high = [40, 20], [40, 20]
        low[j] -= .2
        high[j] += .2
        a = pair(2, 15, low, .32, domain, '-jac-fine')
        b = pair(2, 15, high, .32, domain, '-jac-fine')
        columns.append([real(psi(b, k) - psi(a, k)) / .4 for k in range(2)])
    local_refinements.append(columns)
for j in range(2):
    for k in range(2):
        check('Сетка нелинейного Якобиана %d%d' % (j, k), rel(local[-1][j][k], local_refinements[0][j][k]), .01)
        check('Граница нелинейного Якобиана %d%d' % (j, k), rel(local_refinements[0][j][k], local_refinements[1][j][k]), .01)
finish('coils', {'single': single_results, 'linear_matrices': matrices,
                  'local_matrix_columns_H': local[-1]},
       {'N': [200, 100], 'coil_r_mm': [10, 16], 'coil_length_mm': 20,
        'single_core_r_z_mm': [8, -30, 30], 'pair_core_r_z_mm': [8, -55, 55],
        'gap_definition': 'Между обращёнными торцами катушек', 'pair_domain_r_z_mm': [360, -360, 360],
        'BH_T_Am': BH, 'local_currents_A': [40, 20], 'J_window_relation': 'J=NI/S'})