"""Распределённая МДС: проверка коэффициентов 1/2 и 1/3 в FEMM 5.

О. Б. Буль, с. 18–21: самостоятельный точный предельный пример.
Откройте программу во встроенном Python FEMM 5. Результаты записываются
в новый calculated-NN. Размеры в мм, глубина 20 мм. Это локальная модель
равномерного слоя токов, а не полный электромагнит: на верхней границе
A=0 подразумевается обратный поверхностный ток вне расчётной области.
На трёх остальных сторонах естественное условие dA/dn=0. Путь замыкания
идеализированного тока указан, интегральные величины относятся к слою.
"""
import os
import sys
import math
import json

source = globals().get('__file__', '') or (sys.argv[0] if sys.argv else '')
BASE = globals().get('BUL_OUTPUT_ROOT', os.path.dirname(os.path.abspath(source)) if source else os.getcwd())
number = 1
while os.path.exists(BASE+'/calculated-%02d' % number):
    number += 1
OUT = BASE+'/calculated-%02d' % number
os.mkdir(OUT)
MU0 = 4e-7*math.pi
N = 300
GAP = 50.0
HEIGHT = 40.0
DEPTH = 20.0
COIL = 'Равномерно распределённая обмотка'
MAT = 'Немагнитная область обмотки'
BOUND = 'Обратный ток на границе A=0'
RECORDS = []
CHECKS = []


def real(x):
    return float(complex(x).real)


def check(key, value, limit):
    CHECKS.append({'id': key, 'value': value, 'limit': limit,
                   'passed': bool(math.isfinite(value) and abs(value) <= limit)})


def measure():
    mo_groupselectblock(1)
    energy = real(mo_blockintegral(2))
    coenergy = real(mo_blockintegral(17))
    volume = real(mo_blockintegral(10))
    mo_clearblock()
    cp = mo_getcircuitproperties(COIL)
    a0 = real(mo_getpointvalues(GAP/2, 0.001)[0])
    ah = real(mo_getpointvalues(GAP/2, HEIGHT-0.001)[0])
    bmid = real(mo_getpointvalues(GAP/2, HEIGHT/2)[1])
    return [energy, coenergy, volume, real(cp[2]), a0, ah, bmid]


for k, (mesh, current) in enumerate([(1.0, 1.0), (0.5, 1.0), (0.25, 1.0), (0.5, 2.0)]):
    newdocument(0)
    mi_probdef(0, 'millimeters', 'planar', 1e-10, DEPTH, 25, 0)
    mi_setcomment('Распределённая МДС. Группа 1 — равномерная обмотка в немагнитной области. '
                  'Обратный поверхностный ток замыкается на верхней A=0. '
                  'Три естественные границы dA/dn=0. Локальный аналитический эталон; не полный аппарат. '
                  'N=300, ширина 50 мм, высота 40 мм, глубина 20 мм.')
    mi_addmaterial(MAT, 1, 1, 0, 0, 0)
    mi_addcircprop(COIL, current, 1)
    mi_addboundprop(BOUND, 0, 0, 0, 0, 0, 0, 0, 0, 0)
    corners = [(0, 0), (GAP, 0), (GAP, HEIGHT), (0, HEIGHT)]
    for x, y in corners:
        mi_addnode(x, y)
    for j in range(4):
        a, b = corners[j], corners[(j+1) % 4]
        mi_addsegment(a[0], a[1], b[0], b[1])
        mi_selectsegment((a[0]+b[0])/2, (a[1]+b[1])/2)
        mi_setsegmentprop(BOUND if j == 2 else '<None>', mesh, 0, 0, 1)
        mi_clearselected()
    mi_addblocklabel(GAP/2, HEIGHT/2)
    mi_selectlabel(GAP/2, HEIGHT/2)
    mi_setblockprop(MAT, 0, mesh, COIL, 0, 1, N)
    mi_clearselected()
    filename = '%02d-Распределённая-МДС-сетка-%gмм-ток-%gА.fem' % (k+1, mesh, current)
    path = OUT+'/'+filename
    mi_saveas(path)
    mi_close()
    opendocument(path)
    mi_analyze(1)
    mi_loadsolution()
    mo_smooth('off')
    values = measure()
    elements = int(mo_numelements())
    mo_close()
    mi_close()
    opendocument(path)
    mi_loadsolution()
    mo_smooth('off')
    repeated = measure()
    mo_close()
    mi_close()
    g, h, b = GAP*.001, HEIGHT*.001, DEPTH*.001
    geometric = MU0*b*h/g
    flux_exact = geometric*N*current/2
    psi_exact = geometric*N*N*current/3
    energy_exact = .5*current*psi_exact
    flux = b*(values[4]-values[5])
    record = {'id': 'distributed-%d' % k, 'family': 'distributed-mmf', 'model': filename,
              'solution': filename[:-4]+'.ans', 'mesh_mm': mesh, 'current_A': current,
              'turns': N, 'elements': elements, 'energy_J': values[0], 'coenergy_J': values[1],
              'linkage_Wb': values[3], 'leakage_flux_Wb': flux,
              'geometric_permeance_H': geometric,
              'flux_reduction_factor': flux/(geometric*N*current),
              'linkage_reduction_factor': values[3]/(geometric*N*N*current),
              'exact_flux_Wb': flux_exact, 'exact_linkage_Wb': psi_exact,
              'exact_energy_J': energy_exact, 'volume_m3': values[2],
              'saved_reopened_solved': True,
              'solution_reopened_matches': max(abs(a-c) for a,c in zip(values,repeated)) < 1e-12}
    RECORDS.append(record)
    prefix = record['id']+'-'
    check(prefix+'half-flux', flux/flux_exact-1, .002)
    check(prefix+'third-linkage', values[3]/psi_exact-1, .002)
    check(prefix+'field-energy', values[0]/energy_exact-1, .002)
    check(prefix+'energy-linkage', 2*values[0]/(current*values[3])-1, .002)
    check(prefix+'volume', values[2]/(g*h*b)-1, 1e-8)
    check(prefix+'field-ramp', values[6]/(-MU0*N*current/(2*g))-1, .04)
    check(prefix+'roundtrip', 0 if record['solution_reopened_matches'] else 1, 0)
    with open(OUT+'/progress.json', 'w') as f:
        json.dump({'completed_solves':len(RECORDS),'last':filename}, f)

check('mesh-linkage', RECORDS[2]['linkage_Wb']/RECORDS[1]['linkage_Wb']-1, .001)
check('current-scaling-linkage', RECORDS[3]['linkage_Wb']/(2*RECORDS[1]['linkage_Wb'])-1, 1e-8)
check('current-scaling-energy', RECORDS[3]['energy_J']/(4*RECORDS[1]['energy_J'])-1, 1e-8)
result = {'mode':'full', 'schema':'bul-distributed-mmf-v1', 'passed':all(c['passed'] for c in CHECKS),
          'implementation': str(sys.implementation), 'records':RECORDS, 'checks':CHECKS,
          'source_printed_pages':'Буль, 2005, с. 18–21',
          'scope':'Точный локальный плоский эталон равномерно распределённой МДС. '
                  'Обратный поверхностный ток на A=0. Скалярные коэффициенты 1/2 и 1/3 '
                  'относятся к равномерной намотке и принятому полю; не универсальны.'}
with open(OUT+'/result.json', 'w') as f:
    json.dump(result, f)
print('Распределённая МДС:', len(RECORDS), 'полей;', len(CHECKS), 'проверок;', result['passed'])
if not result['passed']:
    raise ValueError('Не пройдена проверка распределённой МДС')