"""Проводимость и поворотный якорь: самостоятельные модели FEMM 5.

Методический источник: О. Б. Буль, 2005, гл. 3 и приложение П4.
Геометрия, параметры и программы созданы для курса. Это не копии
библиотечных файлов и не воспроизведение численных таблиц книги.
Основной маршрут — Run Python Script во встроенном Python FEMM 5.
Каждый запуск сохраняет новые результаты calculated-NN. Старые не меняет.
"""
import os
import sys
import math
import json

MODE = globals().get('BUL_MODE', 'full')
source_path = globals().get('__file__', '') or (sys.argv[0] if sys.argv else '')
BASE = globals().get('BUL_OUTPUT_ROOT', os.path.dirname(os.path.abspath(source_path)) if source_path 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
AIR = 'Воздух — немагнитная область'
STEEL = 'Сталь — линейная учебная, мю 2000'
COIL = 'Обмотка — немагнитная область'
CIRCUIT = 'Катушка возбуждения — ток 1 А'
ZERO = 'Линия потока — A=0'
FLUX = 'Линия потока — заданное A'
OUTER = 'Внешняя граница — A=0'
GROUPS = {'1':'Неподвижный магнитопровод', '2':'Поворотный якорь',
          '3':'Катушка возбуждения', '4':'Воздух', '5':'Идеальная поверхность полюса'}
CONTRACT = {
 'revision':'bul-completion-fields-v1', 'length_unit':'mm', 'groups':GROUPS,
 'source':'О. Б. Буль, 2005, гл. 3, с. 86–168; П4, с. 287–296.',
 'adaptation':'Самостоятельные размеры, построение и код. Заимствована методика, не исходные файлы и не числа решения.',
 'tube-rectangle':{'depth_mm':20, 'width_mm':20, 'gap_mm':2,
    'A_difference_Wb_per_m':.002, 'boundaries':'A=ΔA слева, A=0 справа; сверху и снизу естественное ∂A/∂n=0, идеальное железо вне области.',
    'analytic':'P=mu0*b*w/g; Φ=b*ΔA; Um=ΔA*g/(mu0*w); B=ΔA/w.'},
 'tube-sector':{'depth_mm':20, 'radii_mm':[20,40], 'angle_deg':30,
    'A_difference_Wb_per_m':.002, 'boundaries':'A=ΔA на внутренней дуге, A=0 на внешней; радиальные стороны — идеальные поверхности полюсов.',
    'analytic':'P=mu0*b*ln(r2/r1)/alpha; Φ=b*ΔA; Um=ΔA*alpha/(mu0*ln(r2/r1)); Bθ=ΔA/(r*ln(r2/r1)).'},
 'pole-plane':{'depth_mm':20, 'outer_arc_radius_mm':50, 'gap_mm':25, 'overlap_width_mm':40,
    'A_difference_Wb_per_m':.002,
    'polygon_mm':[[-50,0],[0,50],[0,25],[40,25],[40,0]],
    'boundaries':'Четверть окружности: A=ΔA; правая вертикаль: A=0. Нижняя плоскость и ломаная верхняя поверхность — естественные идеальные полюсные границы.',
    'normalization':'λ=Φ/(mu0*b*Um)=P/(mu0*b); безразмерная. Удельная проводимость P/b имеет единицу Гн/м.',
    'scope':'Заданная смешанная краевая задача идеального железа. Дуга является границей линии потока; это не расчёт открытой 3D-конструкции.'},
 'pivot-u':{'depth_mm':20, 'turns':300, 'current_A':1, 'relative_permeability':2000,
    'core_polygon_mm':[[-30,0],[30,0],[30,50],[20,50],[20,10],[-20,10],[-20,50],[-30,50]],
    'initial_armature_mm':[[-30,52],[30,52],[30,62],[-30,62]],
    'hinge_mm':[-30,52], 'angle_sign':'Положительный угол против часовой стрелки открывает правый зазор. Момент притяжения отрицателен.',
    'kinematics':'Все четыре вершины, метка и вектор перемещения якоря вращаются вокруг одной неподвижной оси. Якорь не переносится вертикально.',
    'torque':'T_hinge=T_origin-x_hinge*Fy+y_hinge*Fx; метры и Н. T_hinge=dWprime/dtheta при неизменном токе, theta в радианах.',
    'scope':'Плоская модель конечной глубины; без торцевого 3D-рассеяния, механики шарнира и контактов. Линейная сталь выбрана для прозрачной энергетической проверки.'},
 'omissions':['торцевое рассеяние 3D','механика','нагрев','гистерезис','вихревые токи','аттестация промышленного материала']}
with open(OUT+'/geometry-contract.json','w') as handle:
    json.dump(CONTRACT,handle)
RECORDS = []
CHECKS = []
NODES = {}

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

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

def check(name,observed,limit,scope,condition=None):
    CHECKS.append({'id':name,'passed':bool(observed<limit if condition is None else condition),
        'observed':observed,'limit':limit,'scope':scope})

def node(p):
    key = (round(p[0],10),round(p[1],10))
    if key not in NODES:
        mi_addnode(key[0],key[1])
        NODES[key] = True

def segment(a,b,h,group=4,boundary='<None>'):
    node(a)
    node(b)
    mi_addsegment(a[0],a[1],b[0],b[1])
    mi_selectsegment((a[0]+b[0])/2,(a[1]+b[1])/2)
    mi_setsegmentprop(boundary,h,0,0,group)
    mi_clearselected()

def arc(a,b,angle,maxangle,group,boundary,probe):
    node(a)
    node(b)
    mi_addarc(a[0],a[1],b[0],b[1],angle,maxangle)
    mi_selectarcsegment(probe[0],probe[1])
    mi_setarcsegmentprop(maxangle,boundary,0,group)
    mi_clearselected()

def polygon(points,h,group):
    for k in range(len(points)):
        segment(points[k],points[(k+1)%len(points)],h,group)

def rectangle(x0,y0,x1,y1,h,group):
    polygon([[x0,y0],[x1,y0],[x1,y1],[x0,y1]],h,group)

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

def start(depth,da,comment):
    global NODES
    NODES = {}
    newdocument(0)
    mi_probdef(0,'millimeters','planar',1e-10,depth,28,0)
    mi_addmaterial(AIR,1,1,0,0,0)
    mi_addboundprop(ZERO,0,0,0,0,0,0,0,0,0)
    mi_addboundprop(FLUX,da,0,0,0,0,0,0,0,0)
    mi_setcomment(comment+' Все длины мм. Глубина %.9g мм. '%depth+
      'Имена свойств русские. Группы: 1 — неподвижный магнитопровод; 2 — поворотный якорь; 3 — обмотка; 4 — воздух; 5 — идеальная поверхность полюса. '+
      'Самостоятельная учебная адаптация методики Буля; не исходная модель книги.')

def line(points,kind):
    mo_clearcontour()
    for p in points:
        mo_addcontour(p[0],p[1])
    values = mo_lineintegral(kind)
    mo_clearcontour()
    return real(values[0])

def energy():
    mo_clearblock()
    mo_groupselectblock()
    w = real(mo_blockintegral(2))
    wp = real(mo_blockintegral(17))
    mo_clearblock()
    return w,wp

def prepare(stem):
    path = OUT+'/'+stem+'.fem'
    mi_saveas(path)
    mi_close()
    opendocument(path)
    mi_analyze(1)
    mi_loadsolution()
    return path

def finish(row,path,probe):
    row['elements'] = mo_numelements()
    row['nodes'] = mo_numnodes()
    a = real(mo_getpointvalues(probe[0],probe[1])[0])
    mo_close()
    mi_close()
    opendocument(path)
    mi_loadsolution()
    reopened = real(mo_getpointvalues(probe[0],probe[1])[0])
    row['reopened_potential_Wb_per_m'] = reopened
    row['roundtrip_relative_error'] = rel(reopened,a)
    row['solution_reopened_matches'] = row['roundtrip_relative_error']<1e-12
    row['saved_reopened_solved'] = True
    mo_close()
    mi_close()
    RECORDS.append(row)
    check(row['id']+'-roundtrip',row['roundtrip_relative_error'],1e-12,
        'Точная сохранённая русифицированная модель закрыта до solve и решение открыто повторно.')
    with open(OUT+'/progress.json','w') as handle:
        json.dump({'completed_solves':len(RECORDS),'last':row['id']},handle)
    print('BUL_COMPLETION_SOLVED',row['id'],row['elements'])
    return row

def tube(family,h=.3,da=.002,depth=20,arc_step=1,tag='основной'):
    start(depth,da,'Точная полевая трубка: '+family+'. Вне расчётной области железо с бесконечной проницаемостью; свободного тока внутри нет.')
    alpha = math.pi/6
    if family=='tube-rectangle':
        segment((0,0),(20,0),h,5)
        segment((20,0),(20,2),h,4,ZERO)
        segment((20,2),(0,2),h,5)
        segment((0,2),(0,0),h,4,FLUX)
        label(10,1,AIR,h,4)
        probe = [10,1]
        hpath = [[10,.00001],[10,1.99999]]
        cross = [[.00001,1],[19.99999,1]]
        exact_p = MU0*(depth*.001)*.020/.002
        exact_um = da*.002/(MU0*.020)
    else:
        r1,r2 = 20,40
        c,s = math.cos(alpha),math.sin(alpha)
        segment((r1,0),(r2,0),h,5)
        segment((r1*c,r1*s),(r2*c,r2*s),h,5)
        arc((r1,0),(r1*c,r1*s),30,arc_step,4,FLUX,(r1*math.cos(alpha/2),r1*math.sin(alpha/2)))
        arc((r2,0),(r2*c,r2*s),30,arc_step,4,ZERO,(r2*math.cos(alpha/2),r2*math.sin(alpha/2)))
        label(30*math.cos(alpha/2),30*math.sin(alpha/2),AIR,h,4)
        probe = [30*math.cos(alpha/2),30*math.sin(alpha/2)]
        hpath = [[30*math.cos(alpha*k/120),30*math.sin(alpha*k/120)] for k in range(121)]
        cross = [[(r1+.00001)*math.cos(alpha/2),(r1+.00001)*math.sin(alpha/2)],
            [(r2-.00001)*math.cos(alpha/2),(r2-.00001)*math.sin(alpha/2)]]
        exact_p = MU0*depth*.001*math.log(2)/alpha
        exact_um = da*alpha/(MU0*math.log(2))
    stem = '%02d-%s-%s'%(len(RECORDS)+1,family,tag)
    path = prepare(stem)
    phi = abs(line(cross,0))
    um = abs(line(hpath,1))
    w,wp = energy()
    p = phi/um
    row = {'id':family+'-'+tag,'family':family,'tag':tag,'model':stem+'.fem','solution':stem+'.ans',
      'parameters':{'mesh_mm':h,'arc_step_deg':arc_step,'depth_mm':depth,'A_difference_Wb_per_m':da},
      'metrics':{'flux_Wb':phi,'magnetic_voltage_A':um,'permeance_H':p,'energy_J':w,'coenergy_J':wp,
        'exact_flux_Wb':da*depth*.001,'exact_magnetic_voltage_A':exact_um,'exact_permeance_H':exact_p,
        'flux_relative_error':rel(phi,da*depth*.001),'permeance_relative_error':rel(p,exact_p),
        'energy_identity_relative_error':rel(w,.5*phi*um)}}
    check(row['id']+'-flux',row['metrics']['flux_relative_error'],.0002,'Φ=b ΔA, плоский векторный потенциал без лишнего μ0.')
    check(row['id']+'-analytic',row['metrics']['permeance_relative_error'],.005,'Независимая точная аналитическая проводимость геометрии в СИ.')
    check(row['id']+'-energy',row['metrics']['energy_identity_relative_error'],.004,'W=Φ Um/2 для линейного поля без свободного тока.')
    return finish(row,path,probe)

def pole(h=.5,arc_step=.5,da=.002,depth=20,tag='основной'):
    start(depth,da,'Полюс–плоскость: смешанная воздушная краевая задача. A=const на дуге и правой границе, естественное Ht=0 на идеальном железе.')
    segment((0,50),(0,25),h,5)
    segment((0,25),(40,25),h,5)
    segment((40,25),(40,0),h,4,ZERO)
    segment((40,0),(-50,0),h,5)
    arc((0,50),(-50,0),90,arc_step,4,FLUX,(-50/math.sqrt(2),50/math.sqrt(2)))
    label(20,12,AIR,h,4)
    stem = '%02d-pole-plane-%s'%(len(RECORDS)+1,tag)
    path = prepare(stem)
    um_paths = [abs(line([[x,.0001],[x,24.9999]],1)) for x in [10,20,30]]
    um = sum(um_paths)/len(um_paths)
    phi = abs(line([[-49.999,0.0001],[39.999,.0001]],0))
    w,wp = energy()
    p = phi/um
    row = {'id':'pole-plane-'+tag,'family':'pole-plane','tag':tag,'model':stem+'.fem','solution':stem+'.ans',
      'parameters':{'mesh_mm':h,'arc_step_deg':arc_step,'depth_mm':depth,'A_difference_Wb_per_m':da},
      'metrics':{'flux_Wb':phi,'exact_flux_Wb':depth*.001*da,'magnetic_voltage_A':um,
        'path_magnetic_voltages_A':um_paths,'path_relative_spread':(max(um_paths)-min(um_paths))/um,
        'permeance_H':p,'specific_permeance_H_per_m':p/(depth*.001),
        'normalized_permeance_lambda':p/(MU0*depth*.001),'energy_J':w,'coenergy_J':wp,
        'energy_identity_relative_error':rel(w,.5*phi*um),'flux_relative_error':rel(phi,depth*.001*da)}}
    check(row['id']+'-flux',row['metrics']['flux_relative_error'],.0005,'Φ=b ΔA, интеграл по контуру между заданными линиями потока.')
    check(row['id']+'-paths',row['metrics']['path_relative_spread'],.004,'Um между идеальными полюсами по трём путям, в воздухе rot H=0.')
    check(row['id']+'-energy',row['metrics']['energy_identity_relative_error'],.004,'Интегральная энергетическая проверка W=Φ Um/2.')
    return finish(row,path,[20,12])

def rotate(p,angle):
    c,s = math.cos(angle),math.sin(angle)
    x,y = p[0]+30,p[1]-52
    return [-30+c*x-s*y,52+s*x+c*y]

def circuit_estimate(angle_deg):
    theta = angle_deg*math.pi/180
    plegs = []
    for xl,xr in [(-30,-20),(20,30)]:
        if abs(theta)<1e-12:
            plegs.append(MU0*.020*.010/.002)
        else:
            plegs.append(MU0*.020/math.tan(theta)*math.log((2+(xr+30)*math.tan(theta))/(2+(xl+30)*math.tan(theta))))
    rsteel = .140/(MU0*2000*.010*.020)
    phi = 300/(1/plegs[0]+1/plegs[1]+rsteel)
    return plegs,phi,.5*300*phi

def pivot(angle_deg=3,mf=1,df=1,tag='основной'):
    theta = angle_deg*math.pi/180
    start(20,.002,'П-магнит с поворотным якорем. Ось (-30,52) мм; угол %.9g градус. Угол положителен против часовой стрелки.'%angle_deg)
    mi_addmaterial(STEEL,2000,2000,0,0,0)
    mi_addmaterial(COIL,1,1,0,0,0)
    mi_addcircprop(CIRCUIT,1,1)
    mi_addboundprop(OUTER,0,0,0,0,0,0,0,0,0)
    hs,hg,ha = .8*mf,.14*mf,2.2*mf
    polygon(CONTRACT['pivot-u']['core_polygon_mm'],hs,1)
    moved = [rotate(p,theta) for p in CONTRACT['pivot-u']['initial_armature_mm']]
    polygon(moved,hs,2)
    midpoint = rotate([0,52],theta)
    mi_selectsegment(midpoint[0],midpoint[1])
    mi_setsegmentprop('<None>',hg,0,0,2)
    mi_clearselected()
    armature_label = rotate([0,57],theta)
    label(armature_label[0],armature_label[1],STEEL,hs,2)
    label(-25,25,STEEL,hs,1)
    for x0,x1,turns in [(-40,-32,300),(-18,-10,-300)]:
        rectangle(x0,16,x1,44,hs,3)
        label((x0+x1)/2,30,COIL,hs,3,CIRCUIT,turns)
    for x in [-25,25]:
        mi_selectsegment(x,50)
        mi_setsegmentprop('<None>',hg,0,0,1)
        mi_clearselected()
    x0,x1,y0,y1 = -110*df,110*df,-60*df,130*df
    for a,b in [((x0,y0),(x1,y0)),((x1,y0),(x1,y1)),((x1,y1),(x0,y1)),((x0,y1),(x0,y0))]:
        segment(a,b,ha,4,OUTER)
    label(0,-20,AIR,ha,4)
    stem = '%02d-pivot-u-%s'%(len(RECORDS)+1,tag)
    path = prepare(stem)
    current,voltage,psi = mo_getcircuitproperties(CIRCUIT)
    mo_clearblock()
    mo_groupselectblock(2)
    fx,fy,t0,volume = [real(mo_blockintegral(k)) for k in [18,19,22,10]]
    mo_clearblock()
    th = t0+.030*fy+.052*fx
    w,wp = energy()
    # Вертикальная высота нижней поверхности над полюсами при данном x.
    gap_left = 2+5*math.tan(theta)
    gap_right = 2+55*math.tan(theta)
    # Трубки без выпучивания: интеграл dP=mu0*b*dx/g(x).
    # Точная для этих вертикальных трубок геометрическая проводимость, не точное поле.
    plegs,estimate_flux,estimate_energy = circuit_estimate(angle_deg)
    estimate_torque = (circuit_estimate(angle_deg+.01)[2]-circuit_estimate(angle_deg-.01)[2])/(.02*math.pi/180)
    gap_fluxes = [abs(line([[xl+.001,51],[xr-.001,51]],0)) for xl,xr in [(-30,-20),(20,30)]]
    row = {'id':'pivot-u-'+tag,'family':'pivot-u','tag':tag,'model':stem+'.fem','solution':stem+'.ans',
      'parameters':{'angle_deg':angle_deg,'hinge_mm':[-30,52],'mesh_factor':mf,'domain_factor':df,
        'mesh_mm':{'steel':hs,'pole_edges':hg,'air':ha},'armature_vertices_mm':moved},
      'metrics':{'Fx_N':fx,'Fy_N':fy,'torque_origin_Nm':t0,'torque_hinge_Nm':th,
        'energy_J':w,'coenergy_J':wp,'linkage_Wb':real(psi),'current_A':real(current),
        'energy_identity_relative_error':rel(w+wp,real(current)*real(psi)),
        'armature_volume_m3':volume,'armature_volume_exact_m3':.06*.01*.02,
        'armature_volume_relative_error':rel(volume,.06*.01*.02),
        'center_vertical_gaps_mm':[gap_left,gap_right],
        'tube_gap_permeances_H':plegs,'circuit_flux_estimate_Wb':estimate_flux,
        'gap_fluxes_Wb':gap_fluxes,'circuit_coenergy_estimate_J':estimate_energy,
        'circuit_torque_estimate_Nm':estimate_torque,
        'circuit_scope':'Вертикальные трубки, нет выпучивания/рассеяния. Равномерная сталь длиной 140 мм — только явно заданная инженерная оценка.'}}
    check(row['id']+'-volume',row['metrics']['armature_volume_relative_error'],1e-8,'Объём повернутого якоря не меняется; область момента группа 2.')
    check(row['id']+'-energy',row['metrics']['energy_identity_relative_error'],.003,'W+Wprime=I Psi, линейная одна цепь без постоянных магнитов.')
    check(row['id']+'-restoring',th,0,'При положительном раскрывающем угле электромагнитный момент закрывает зазор.')
    return finish(row,path,armature_label)

rect = tube('tube-rectangle')
sector = tube('tube-sector',h=.5,arc_step=1)
pole_base = pole(h=1,arc_step=1,tag='сетка-грубая')
pivot_base = pivot(tag='сетка-грубая')
derivatives = []
if MODE!='smoke':
    sector_fine = tube('tube-sector',h=.25,arc_step=.25,tag='сетка-мелкая')
    check('tube-sector-convergence',rel(sector['metrics']['permeance_H'],sector_fine['metrics']['permeance_H']),.003,'Сетка воздуха и аппроксимация дуг сгущены одновременно; точный эталон проверен отдельно.')
    pm = pole(h=.5,arc_step=.5,tag='сетка-средняя')
    pf = pole(h=.25,arc_step=.25,tag='сетка-мелкая')
    check('pole-plane-mesh',rel(pm['metrics']['normalized_permeance_lambda'],pf['metrics']['normalized_permeance_lambda']),.003,'Последняя пара сеток; углы дуг также уменьшены.')
    ps = pole(h=.25,arc_step=.25,da=.004,depth=40,tag='масштаб-потока-глубины')
    check('pole-plane-lambda-invariance',rel(ps['metrics']['normalized_permeance_lambda'],pf['metrics']['normalized_permeance_lambda']),1e-8,'A увеличено 2 раза, глубина 2 раза: Φ в 4 раза, Um в 2, λ неизменна.')
    check('pole-plane-flux-scaling',rel(ps['metrics']['flux_Wb'],4*pf['metrics']['flux_Wb']),1e-8,'Точное плоское масштабирование глубины и возбуждения.')
    pivot_medium = pivot(mf=.7,tag='сетка-средняя')
    pivot_fine = pivot(mf=.5,tag='сетка-мелкая')
    pivot_domain = pivot(mf=.5,df=1.4,tag='область-больше')
    check('pivot-u-mesh',rel(pivot_medium['metrics']['torque_hinge_Nm'],pivot_fine['metrics']['torque_hinge_Nm']),.02,'Момент вокруг шарнира; последняя пара сеток.')
    check('pivot-u-domain',rel(pivot_domain['metrics']['torque_hinge_Nm'],pivot_fine['metrics']['torque_hinge_Nm']),.015,'Внешний прямоугольник увеличен 1.4 раза, локальная сетка неизменна.')
    for step in [.1,.05]:
        minus = pivot(3-step,.5,tag='угол-минус-%g'%step)
        plus = pivot(3+step,.5,tag='угол-плюс-%g'%step)
        torque = (plus['metrics']['coenergy_J']-minus['metrics']['coenergy_J'])/(2*step*math.pi/180)
        derivatives.append({'step_deg':step,'torque_from_coenergy_Nm':torque,
           'torque_weighted_Maxwell_Nm':pivot_fine['metrics']['torque_hinge_Nm'],
           'relative_difference':rel(torque,pivot_fine['metrics']['torque_hinge_Nm']),
           'minus_model':minus['model'],'plus_model':plus['model']})
    check('pivot-u-torque-methods',derivatives[-1]['relative_difference'],.04,'Взвешенный Maxwell и dWprime/dtheta при постоянном I; theta в радианах, момент перенесён к шарниру.')
    check('pivot-u-derivative-step',rel(derivatives[0]['torque_from_coenergy_Nm'],derivatives[1]['torque_from_coenergy_Nm']),.02,'Два независимых центральных шага по углу.')
    for angle in [0,1,6]:
        pivot(angle,.7,tag='характеристика-%g-град'%angle)
with open(OUT+'/pivot-u-energy-derivative.json','w') as handle:
    json.dump(derivatives,handle)
result = {'schema':'bul-completion-fields-native-v1','implementation':repr(sys.implementation),'mode':MODE,
 'records':RECORDS,'checks':CHECKS,'passed':all(c['passed'] for c in CHECKS),'derivatives':derivatives,
 'lookup':{r['id']:r for r in RECORDS},'geometry_revision':CONTRACT['revision'],
 'scope':'Собственные 2D магнитостатические адаптации методик книги. Проверенные поля FEMM 5, не 3D/физические испытания аппарата.'}
with open(OUT+'/result.json','w') as handle:
    json.dump(result,handle)
print('BUL_COMPLETION_PASSED',result['passed'],len(RECORDS),OUT)
if not result['passed']:
    raise RuntimeError('Не пройдены checks: сохраните квитанцию и изучите result.json')