"""Собственный учебный опыт: поле рассеяния над поперечным пазом.

Запуск: FEMM 5 → Run Script → сохранённый reproduce.py.
Встроенный Python; команды mi_*/mo_* доступны глобально.
Новый запуск создаёт calculated-NN и не перезаписывает опубликованные файлы.
Плоская магнитостатика, мм; условное намагничивание удалённой границей.
Это не модель реального ярма, гистерезиса или локальной трёхмерной трещины.
"""
import os
import sys
import math
import json

saved = globals().get('__file__', '') or (sys.argv[0] if sys.argv else '')
base = os.path.dirname(os.path.abspath(saved)) if saved else os.getcwd()
index = 1
while os.path.exists(base + '/calculated-%02d' % index):
    index += 1
directory = base + '/calculated-%02d' % index
os.mkdir(directory)
records = []
BH = [(0, 0), (0.5, 100), (1.0, 200), (1.4, 400),
      (1.6, 1000), (1.8, 5000), (2.0, 20000), (2.2, 80000)]
CONFIG = {'plate_width_mm': 80, 'plate_height_mm': 10, 'depth_mm': 10,
          'notch_width_mm': 1, 'notch_depth_mm': 2,
          'linear_relative_permeability': 500, 'far_field_T': 0.01,
          'scan_x_mm': [-12 + k * 0.5 for k in range(49)],
          'lift_off_mm': [1, 2, 4], 'BH_T_A_per_m': BH,
          'BH_origin': 'Авторская условная характеристика; не паспорт марки стали.'}


def polygon(points, boundary=None, spacing=0):
    for x, y in points:
        mi_addnode(x, y)
    for i, (x, y) in enumerate(points):
        x2, y2 = points[(i + 1) % len(points)]
        mi_addsegment(x, y, x2, y2)
        if boundary:
            mi_selectsegment((x + x2) / 2, (y + y2) / 2)
            mi_setsegmentprop(boundary, spacing, 0, 0, 0)
            mi_clearselected()


def label(x, y, name, mesh, group):
    mi_addblocklabel(x, y)
    mi_selectlabel(x, y)
    mi_setblockprop(name, 0, mesh, '<None>', 0, group, 0)
    mi_clearselected()


def field(notch, level, outer, bfar, nonlinear=False, tag='base'):
    newdocument(0)
    mi_probdef(0, 'millimeters', 'planar', 1e-9, CONFIG['depth_mm'], 30, 0)
    mi_addmaterial('Воздух', 1, 1, 0, 0, 0)
    metal = 'Условная сталь: нелинейная' if nonlinear else 'Условная сталь: линейная'
    mi_addmaterial(metal, 500, 500, 0, 0, 0)
    if nonlinear:
        for b, h in BH:
            mi_addbhpoint(metal, b, h)
    # Внешняя функция A = Bfar*y. FEMM использует мм в координатном полиноме.
    # Bx = dA/dy; коэффициент в Вб/(м·мм) равен Bfar/1000.
    boundary = 'Удалённое поле вдоль x'
    mi_addboundprop(boundary, 0, 0, bfar / 1000, 0, 0, 0, 0, 0, 0)
    polygon([(-outer, -outer), (outer, -outer), (outer, outer), (-outer, outer)], boundary, outer/12)
    # Вспомогательная граница одного и того же воздуха задаёт локальную сетку.
    polygon([(-45, -10), (45, -10), (45, 15), (-45, 15)])
    half_width = CONFIG['plate_width_mm'] / 2
    half_height = CONFIG['plate_height_mm'] / 2
    half_notch = CONFIG['notch_width_mm'] / 2
    bottom_notch = half_height - CONFIG['notch_depth_mm']
    # Одинаковые подобласти и ограничения сетки в целой и дефектной модели.
    # Здоровый вариант заполняет подобласть паза той же сталью, дефектный — воздухом.
    polygon([(-half_width, -half_height), (half_width, -half_height),
             (half_width, half_height), (half_notch, half_height),
             (-half_notch, half_height), (-half_width, half_height)])
    for x, y in [(-half_notch, bottom_notch), (half_notch, bottom_notch)]:
        mi_addnode(x, y)
    mi_addsegment(-half_notch, half_height, -half_notch, bottom_notch)
    mi_addsegment(-half_notch, bottom_notch, half_notch, bottom_notch)
    mi_addsegment(half_notch, bottom_notch, half_notch, half_height)
    scale = 2 ** level
    label(0, 0, metal, 0.8 / scale, 1)
    label(0, half_height-CONFIG['notch_depth_mm']/2,
          'Воздух' if notch else metal, 0.25 / scale, 4)
    label(0, 10, 'Воздух', 0.8 / scale, 2)
    label(0, outer * 0.7, 'Воздух', outer / 10, 3)
    # Ограничиваем ребро паза отдельно; точечный максимум в углу не измеряем.
    for x, y in [(half_notch, (half_height+bottom_notch)/2),
                 (0, bottom_notch), (-half_notch, (half_height+bottom_notch)/2)]:
        mi_selectsegment(x, y)
        mi_setsegmentprop('<None>', 0.25 / scale, 0, 0, 1)
        mi_clearselected()
    stem = '%s-%s-m%d' % (tag, 'notch' if notch else 'whole', level + 1)
    path = directory + '/' + stem + '.fem'
    mi_saveas(path)
    mi_close()
    opendocument(path)
    mi_analyze(1)
    mi_loadsolution()
    mo_smooth('on')
    scans = []
    for height in CONFIG['lift_off_mm']:
        values = []
        for x in CONFIG['scan_x_mm']:
            p = mo_getpointvalues(x, half_height + height)
            values.append([x, float(p[1]), float(p[2])])
        scans.append({'lift_off_mm': height, 'points_x_Bx_By': values})
    p = mo_getpointvalues(half_width/2, 0)
    row = {'model': stem + '.fem', 'notch': bool(notch), 'mesh_level': level + 1,
           'outer_halfwidth_mm': outer, 'far_field_T': bfar, 'nonlinear': nonlinear,
           'saved_reopened_solved': True, 'elements': mo_numelements(),
           'nodes': mo_numnodes(), 'scans': scans,
           'B_sample_T': math.sqrt(float(p[1]) ** 2 + float(p[2]) ** 2),
           'scope': 'Плоский паз, неограниченный вдоль глубины; без гистерезиса.'}
    records.append(row)
    mo_close()
    mi_close()
    print(stem, row['elements'])
    return row


def difference(whole, notch):
    scans = []
    for a, b in zip(whole['scans'], notch['scans']):
        points = [[p[0], q[1]-p[1], q[2]-p[2]]
                  for p, q in zip(a['points_x_Bx_By'], b['points_x_Bx_By'])]
        scans.append({'lift_off_mm': a['lift_off_mm'], 'points_x_deltaBx_deltaBy': points,
                      'deltaBy_peak_to_peak_T': max(p[2] for p in points)-min(p[2] for p in points),
                      'deltaBx_peak_abs_T': max(abs(p[1]) for p in points)})
    return {'whole_model': whole['model'], 'notch_model': notch['model'],
            'mesh_level': whole['mesh_level'], 'outer_halfwidth_mm': whole['outer_halfwidth_mm'],
            'far_field_T': whole['far_field_T'], 'nonlinear': whole['nonlinear'], 'scans': scans}


signals = []
for level in range(3):
    a = field(False, level, 160, 0.01)
    b = field(True, level, 160, 0.01)
    signals.append(difference(a, b))
a = field(False, 2, 240, 0.01, tag='outer')
b = field(True, 2, 240, 0.01, tag='outer')
signals.append(difference(a, b))
for bfar in [0.01, 0.1, 0.2, 0.25]:
    tag = 'nonlinear-%03d' % round(bfar * 1000)
    a = field(False, 2, 160, bfar, True, tag)
    b = field(True, 2, 160, bfar, True, tag)
    signals.append(difference(a, b))

def ratio(a, b):
    return abs(a-b) / max(abs(b), 1e-20)

checks = [
    {'name': 'Конечность всех отсчётов', 'passed': all(math.isfinite(v) for r in records for s in r['scans'] for p in s['points_x_Bx_By'] for v in p)},
    {'name': 'Рост зазора ослабляет размах deltaBy', 'passed': all(a['deltaBy_peak_to_peak_T'] > b['deltaBy_peak_to_peak_T'] for a,b in zip(signals[2]['scans'][:-1],signals[2]['scans'][1:]))},
]
for i, height in enumerate(CONFIG['lift_off_mm']):
    base = signals[2]['scans'][i]['deltaBy_peak_to_peak_T']
    mesh = signals[1]['scans'][i]['deltaBy_peak_to_peak_T']
    outer = signals[3]['scans'][i]['deltaBy_peak_to_peak_T']
    checks.append({'name': 'Положительный разностный сигнал', 'lift_off_mm': height,
                   'value_T': base, 'passed': base > 1e-6})
    checks.append({'name': 'Последние две сетки: размах deltaBy', 'lift_off_mm': height,
                   'relative_change': ratio(mesh, base), 'limit': 0.05, 'passed': ratio(mesh, base) < 0.05})
    checks.append({'name': 'Удаление границы 160→240 мм: размах deltaBy', 'lift_off_mm': height,
                   'relative_change': ratio(outer, base), 'limit': 0.05, 'passed': ratio(outer, base) < 0.05})
result = {'schema': 'femm5-mfl-notch-teaching-v1', 'runtime': 'FEMM 5 / embedded Python',
          'parameters': CONFIG, 'records': records, 'signals': signals, 'checks': checks,
          'passed': all(c['passed'] for c in checks),
          'scope': 'Самостоятельная учебная модель MFL. Сходимость размаха проверяется на всех трёх высотах линейной базовой серии. Нелинейные варианты служат отдельным исследованием на финальной сетке; для них индивидуальная сходимость не заявляется. Расчётная проверка не является физической аттестацией дефектоскопа.'}
with open(directory + '/result.json', 'w') as stream:
    json.dump(result, stream)
print('RESULT', directory, result['passed'])
if not result['passed']:
    raise RuntimeError('Проверки MFL не пройдены; изучите result.json без изменения допусков задним числом.')