"""Дополнение по О. Б. Булю: П6 и частично затенённый полюс.

Основной маршрут: Run Python Script во встроенном Python FEMM 5.
Каждый запуск создаёт новый calculated-NN, прежние поля не меняются.
П6: собственная осесимметричная симметричная геометрия; проводимости
трубок получаются из отдельного почти идеального линейного поля.
AC: собственный плоский длинный средний участок. Две стороны витка
и внешние торцевые соединения объявлены явно. Это двумерный аналог,
не трёхмерный расчёт контактного аппарата и не осесимметричное кольцо.
Поля линейного AC-аналога рассчитаны при DC как базис L/M. Комплексные
токи определяются по внешней цепи, а поле — строгой суперпозицией.
Вся численная часть, включая цепи, выполняется встроенным Python.
"""
import os, sys, math, json

MODE=globals().get('BUL_MODE','full')
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())
k=1
while os.path.exists(BASE+'/calculated-%02d'%k): k+=1
OUT=BASE+'/calculated-%02d'%k; os.mkdir(OUT)
MU0=4e-7*math.pi
AIR='Воздух'; COIL='Обмотка — внешний расчёт сопротивления'
IDEAL='Сталь — почти идеальная, мю 100000'
STEEL='Сталь — условная нелинейная'
ACSTEEL='Сталь — линейная AC, мю 1500'
BC='Внешняя граница A=0'
PRIMARY='Основная обмотка'; SECONDARY='Частичный затеняющий виток'
BH=[[0,0],[.4,.4/(MU0*2000)],[.8,.8/(MU0*2000)],
    [1.2,550],[1.45,1000],[1.6,2200],[1.75,6000],[1.9,18000],
    [2.1,65000],[2.3,180000],[2.6,400000],[3,750000],[4,1545775]]
RECORDS=[]; CHECKS=[]; NODES={}; EDGES={}
def save(name,value):
    with open(OUT+'/'+name,'w') as f: json.dump(value,f)
def check(name,passed,observed,limit,scope):
    CHECKS.append({'id':name,'passed':bool(passed),'observed':observed,'value':observed,'limit':limit,'scope':scope})
def real(x): return float(complex(x).real)
def pair(x): x=complex(x); return [x.real,x.imag]
def rel(a,b): return abs(a-b)/max(abs(b),1e-25)
def node(p):
    p=(round(p[0],9),round(p[1],9))
    if p not in NODES: mi_addnode(p[0],p[1]); NODES[p]=True
    return p
def edge(a,b,h=1,group=0,boundary='<None>'):
    a=node(a); b=node(b); 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,h,0,0,group); mi_clearselected()
def poly(points,h=1,group=0):
    for i in range(len(points)): edge(points[i],points[(i+1)%len(points)],h,group)
def rect(x0,y0,x1,y1,h=1,group=0): poly([(x0,y0),(x1,y0),(x1,y1),(x0,y1)],h,group)
def label(x,y,mat,h=1,group=0,circuit='<None>',turns=0):
    mi_addblocklabel(x,y); mi_selectlabel(x,y); mi_setblockprop(mat,0,h,circuit,0,group,turns); mi_clearselected()
def start(axi=False):
    global NODES,EDGES
    NODES={}; EDGES={}; newdocument(0)
    mi_probdef(0,'millimeters','axi' if axi else 'planar',1e-9,30,25,0)
    for name,mu in [(AIR,1),(COIL,1),(IDEAL,100000),(STEEL,2000),(ACSTEEL,1500)]: mi_addmaterial(name,mu,mu,0,0,0)
    for b,h in BH: mi_addbhpoint(STEEL,b,h)
    mi_addboundprop(BC,0,0,0,0,0,0,0,0,0)
def line(a,b,kind):
    mo_clearcontour(); mo_addcontour(a[0],a[1]); mo_addcontour(b[0],b[1]); x=mo_lineintegral(kind); mo_clearcontour(); return real(x[0])
def phi(z): return abs(line((0,z),(8,z),0))
def hpath(points):
    return abs(sum(line(points[i],points[i+1],1) for i in range(len(points)-1)))
def getb(x,y):
    p=mo_getpointvalues(x,y); return [p[1],p[2]]
def p6build(current,ideal=False,mesh=1):
    start(True)
    mi_setcomment('Собственный симметричный двухполюсный магнит для П6 Буля. Половина z>=0; на z=0 естественная симметрия Br=0. Полный зазор 1 мм; полюс r=8 мм. Группы: 1 — непрерывные полюс и ярмо; 3 — обмотка, 150 витков в половине; 4 — воздух. Условная B(H), не марка стали. Вторая зеркальная обмотка даёт полную МДС 300I.')
    mi_addcircprop(PRIMARY,current,1)
    # Один непрерывный многоугольник: нет перекрывающихся рёбер.
    poly([(0,.5),(8,.5),(8,20),(15,20),(15,40),(45,40),
          (45,0),(55,0),(55,55),(0,55)],.45*mesh,1)
    rect(17,22,35,38,.9*mesh,3)
    mat=IDEAL if ideal else STEEL
    label(4,10,mat,.45*mesh,1)
    label(25,30,COIL,.9*mesh,3,PRIMARY,150)
    # z=0 остаётся естественной границей; r=0 — A=0.
    for z0,z1 in [(0,.5),(.5,20),(20,40),(40,55),(55,100)]: edge((0,z0),(0,z1),.3*mesh,4,BC)
    edge((0,100),(100,100),3*mesh,4,BC); edge((100,100),(100,0),3*mesh,4,BC)
    edge((100,0),(55,0),2*mesh,4); edge((45,0),(0,0),.35*mesh,4)
    label(70,35,AIR,2*mesh,4); label(30,10,AIR,1*mesh,4)
    mi_selectsegment(4,.5); mi_setsegmentprop('<None>',.08*mesh,0,0,1); mi_clearselected()
def begin_solve(tag):
    stem='%02d-%s'%(len(RECORDS)+1,tag); path=OUT+'/'+stem+'.fem'
    mi_saveas(path); mi_close(); opendocument(path); mi_analyze(1); mi_loadsolution()
    return stem,path
def finish(stem,path,row):
    row.update(model=stem+'.fem',solution=stem+'.ans',elements=mo_numelements(),nodes=mo_numnodes(),saved_reopened_solved=True)
    linkage=complex(mo_getcircuitproperties(PRIMARY)[2]); row['linkage_Wb']=linkage.real; row['native_linkage_complex_Wb']=pair(linkage)
    mo_close(); mi_close(); opendocument(path); mi_loadsolution()
    row['russian_circuit_reopened']=rel(complex(mo_getcircuitproperties(PRIMARY)[2]),linkage)<1e-10
    row['solution_reopened_matches']=row['russian_circuit_reopened']; row['id']=stem
    mo_close(); mi_close(); RECORDS.append(row); save(stem+'-measurements.json',row); save('progress.json',{'completed_solves':len(RECORDS),'last':stem}); print('SOLVED',stem,row['elements'])
    return row
def p6solve(tag,current,ideal=False,mesh=1):
    p6build(current,ideal,mesh); stem,path=begin_solve(tag)
    zs=[.501+(19.998-.501)*i/64 for i in range(65)]
    values=[phi(z) for z in zs]
    gap_b=real(getb(.5,.0001)[1])
    gap_u=hpath([(.5,.0001),(.5,.4999)])
    pole_u=hpath([(1,.501),(1,19.999)])
    tail_u=hpath([(1,20.001),(1,47),(50,47),(50,.0001)])
    row={'family':'П6-двухполюсный','tag':tag,'current_A':current,'half_mmf_At':150*current,'ideal':ideal,'mesh_factor':mesh,'gap_Bcenter_T':gap_b,'gap_U_At':gap_u,'pole_U_centerline_At':pole_u,'tail_U_centerline_At':tail_u,'sections_z_mm':zs,'section_flux_Wb':values,'section_average_B_T':[v/(math.pi*.008**2) for v in values],'yoke_B_T':abs(real(getb(50,25)[1]))}
    probes=[(r,z) for r in [1,4,8,12] for z in [25,30,35]]
    probes += [(r,z) for r in [1,5,15,30,45,50] for z in [43,47,52]]
    probes += [(r,z) for r in [47,50,53] for z in [5,15,25,35]]
    row['yoke_bulk_probes']=[{'point_mm':[r,z],'B_abs_T':math.sqrt(sum(real(v)**2 for v in getb(r,z)))} for r,z in probes]
    row['transition_shank_probe_max_B_T']=max(p['B_abs_T'] for p in row['yoke_bulk_probes'][:12])
    row['yoke_bulk_probe_max_B_T']=max(p['B_abs_T'] for p in row['yoke_bulk_probes'][12:])
    row['ampere_path_relative_error']=rel(gap_u+pole_u+tail_u,150*current)
    return finish(stem,path,row)
def h_of_b(b):
    b=abs(b)
    for j in range(1,len(BH)):
        if b<=BH[j][0]:
            b0,h0=BH[j-1]; b1,h1=BH[j]; return h0+(b-b0)*(h1-h0)/(b1-b0)
    return BH[-1][1]+(b-BH[-1][0])/MU0
def p6method(template,target_b,n,method):
    # Λ_i=(Φ_{i+1}-Φ_i)/Uδ: геометрические трубки отдельного поля.
    # П6.1–П6.9: H(B) и средняя H соседних сечений. П6.13–14:
    # предварительное полусечение уточняет H до середины участка.
    ug=target_b*.0005/MU0
    scale=ug/template['gap_U_At']
    flux0=template['section_flux_Wb'][0]*scale
    area=math.pi*.008**2; u=ug; p=flux0; rows=[]
    stride=64//n
    for i in range(n):
        j=i*stride; q=(i+1)*stride
        length=(template['sections_z_mm'][q]-template['sections_z_mm'][j])*.001
        lam=(template['section_flux_Wb'][q]-template['section_flux_Wb'][j])/template['gap_U_At']
        hi=h_of_b(p/area)
        if method==1:
            umid=u+hi*length/2
        else:
            pmid=p+lam*u/2
            hmid=h_of_b(pmid/area)
            umid=u+(hi+hmid)*length/4
        pnext=p+lam*umid; hn=h_of_b(pnext/area)
        du=(hi+hn)*length/2
        rows.append({'z0_mm':template['sections_z_mm'][j],'z1_mm':template['sections_z_mm'][q],'lambda_H':lam,'phi_in_Wb':p,'phi_out_Wb':pnext,'B_in_T':p/area,'B_out_T':pnext/area,'U_tube_At':umid,'delta_U_At':du})
        p=pnext; u+=du
    # Широкое ярмо в слабом поле: μr=2000, геометрический tail R
    # взят с почти идеального поля и приведён по отношению μ.
    tail_r=template['tail_U_centerline_At']/template['section_flux_Wb'][-1]*100000/2000
    tail_u=tail_r*p
    return {'method':method,'sections':n,'target_Bcenter_T':target_b,'gap_U_At':ug,'pole_U_At':u-ug,'tail_U_At':tail_u,'half_mmf_At':u+tail_u,'full_mmf_At':2*(u+tail_u),'rows':rows,'source':'Буль, П6.1–П6.14, печатные с.319–324; перенесено на собственные трубки и геометрию.'}
def p6campaign(target=1.25,stress=False,template=None):
    prefix='П6-%gТл-'%target
    if template is None:
        template=p6solve('П6-идеальные-трубки',1,True,.7)
        check('П6-positive-leakage',min(template['section_flux_Wb'][i+1]-template['section_flux_Wb'][i] for i in range(64))>0,min(template['section_flux_Wb'][i+1]-template['section_flux_Wb'][i] for i in range(64)),0,'В узком полюсе поток увеличивается к ярму; направления трубок проверены.')
    results=[p6method(template,target,n,m) for m in [1,2] for n in [4,8,16,32,64]]
    for m in [1,2]:
        r=[x for x in results if x['method']==m]
        e=rel(r[-2]['half_mmf_At'],r[-1]['half_mmf_At'])
        check(prefix+'method-%d-section-convergence'%m,e<.025,e,.025,'32 и 64 участка при неизменных геометрических трубках.')
    predicted=(results[-1]['half_mmf_At']+results[4]['half_mmf_At'])/2
    # Независимое нелинейное поле при предсказанной МДС; затем локальная
    # секущая по току подгоняет только измеренный Bδ к target.
    row=p6solve(prefix+'нелинейный-прогноз',predicted/150,False,.7)
    trial=p6solve(prefix+'нелинейный-секущая',predicted/150*1.035,False,.7)
    for i in range(3):
        b0=row['gap_Bcenter_T']; b1=trial['gap_Bcenter_T']; i0=row['current_A']; i1=trial['current_A']
        current=max(.1,i1+(target-b1)*(i1-i0)/(b1-b0))
        row,trial=trial,p6solve(prefix+'нелинейный-уточнение-%d'%(i+1),current,False,.7)
    fine=p6solve(prefix+'нелинейный-сетка',trial['current_A'],False,.5)
    check(prefix+'target-B',rel(trial['gap_Bcenter_T'],target)<.005,trial['gap_Bcenter_T'],target,'Независимое FEMM поле при найденном токе, центр воздушного промежутка.')
    check(prefix+'field-mesh',rel(fine['gap_Bcenter_T'],trial['gap_Bcenter_T'])<.012,rel(fine['gap_Bcenter_T'],trial['gap_Bcenter_T']),.012,'Сгущение сетки 0.7 → 0.5 при неизменном токе.')
    check(prefix+'ampere-path',fine['ampere_path_relative_error']<.025,fine['ampere_path_relative_error'],.025,'Интеграл H по пути зазор–полюс–ярмо охватывает 150I.')
    check(prefix+'yoke-unsaturated',fine['yoke_bulk_probe_max_B_T']<.8,fine['yoke_bulk_probe_max_B_T'],.8,'30 контрольных точек широкого фланца и наружного возврата. Ещё 12 точек переходной шейки опубликованы отдельно: она может быть умеренно нелинейной. Это проверка выборки, не доказательство глобального максимума B.')
    comparisons=[]
    for x in results:
        comparisons.append({'method':x['method'],'sections':x['sections'],'half_mmf_At':x['half_mmf_At'],'FEMM_half_mmf_At':fine['half_mmf_At'],'mmf_relative_error':rel(x['half_mmf_At'],fine['half_mmf_At']),'pole_relative_error':rel(x['pole_U_At'],fine['pole_U_centerline_At'])})
    # Не маскируем модельную погрешность цепной редукции сеточной.
    error=max(comparisons[4]['mmf_relative_error'],comparisons[-1]['mmf_relative_error'])
    if stress:
        check(prefix+'reduction-limit-detected',error>.20,error,.20,'Сетка/дробление сходятся, но геометрическая редукция теряет точность: ошибка не принимается за точный расчёт.')
    else:
        check(prefix+'reduction-bound',error<.08,error,.08,'Проверка общей МДС основного примера; отдельно публикуется ошибка МДС именно полюса.')
    save('p6-stress-results.json' if stress else 'p6-results.json',{'template_model':template['model'],'algorithms':results,'field':fine,'comparison':comparisons,'scope':'Новые размеры, условная B(H); Λ из почти идеального поля не учитывает перестройку трубок при насыщении. H(B) цепи линейно интерполируется между опубликованными узлами; FEMM использует собственную интерполяцию той же таблицы. Расхождение включает также это различие, не только трубки и неоднородность по сечению. Применимость оценивается сравнением с независимым FEMM, а не переносом утверждения книги о 5 Тл.'})
    return template

def acbuild(i1,i2,mesh=1):
    start(False)
    mi_setcomment('Плоский модельный аналог частичного затенения полюса. Глубина 30 мм. Основная обмотка N=200. КЗВ: N=1, две стороны с противоположными знаками; ток замыкается вне плоскости торцевыми перемычками, их полное R задаёт внешняя цепь. Открытый паз в левом полюсе x=-21..-19,y=46..50; затенена полоса x=-19..-10. Группы: 1 — ярмо и полюса; 2 — якорь; 3 — основная обмотка; 4 — КЗВ; 5 — воздух. Сталь линейная, гистерезис и skin/proximity не моделируются.')
    mi_addcircprop(PRIMARY,i1,1); mi_addcircprop(SECONDARY,i2,1)
    poly([(-30,0),(30,0),(30,50),(10,50),(10,10),(-10,10),
          (-10,50),(-19,50),(-19,46),(-21,46),(-21,50),(-30,50)],.7*mesh,1)
    rect(-30,50.5,30,60.5,.7*mesh,2)
    for x0,x1,turns in [(-40,-32,200),(-8,0,-200)]:
        rect(x0,15,x1,35,.8*mesh,3); label((x0+x1)/2,25,COIL,.8*mesh,3,PRIMARY,turns)
    # Проводник в пазу и проводник снаружи: внешний торцевой возврат
    # включён в R2. Паз — часть геометрии обеих сравниваемых схем.
    for x0,x1,turns in [(-21,-19,1),(-9,-7,-1)]:
        rect(x0,46,x1,50,.25*mesh,4); label((x0+x1)/2,48,COIL,.25*mesh,4,SECONDARY,turns)
    label(-25,25,ACSTEEL,.7*mesh,1); label(0,55,ACSTEEL,.7*mesh,2)
    for x in [-25,25]:
        mi_selectsegment(x,50); mi_setsegmentprop('<None>',.1*mesh,0,0,1); mi_clearselected()
    mi_selectsegment(0,50.5); mi_setsegmentprop('<None>',.1*mesh,0,0,2); mi_clearselected()
    for a,b in [((-100,-50),(100,-50)),((100,-50),(100,110)),((100,110),(-100,110)),((-100,110),(-100,-50))]: edge(a,b,2.5*mesh,5,BC)
    label(0,-20,AIR,2*mesh,5)
def acsolve(tag,i1,i2,mesh=1,extract=True):
    acbuild(i1,i2,mesh); stem,path=begin_solve(tag)
    props=[mo_getcircuitproperties(x) for x in [PRIMARY,SECONDARY]]
    row={'family':'AC-частичный-полюс','tag':tag,'currents_A':[i1,i2],'linkages_Wb':[real(p[2]) for p in props],'mesh_factor':mesh}
    mo_clearblock(); mo_groupselectblock(); row['energy_J']=real(mo_blockintegral(2)); mo_clearblock()
    mo_groupselectblock(2); row['armature_attraction_WST_N']=-real(mo_blockintegral(19)); mo_clearblock()
    if extract:
        # Поле в середине зазора; одинаковые абсциссы для обеих баз.
        xs=[-60+120*(i+.5)/1200 for i in range(1200)]
        row['gap_x_mm']=xs; row['gap_B_T']=[[real(x) for x in getb(t,50.25)] for t in xs]
        row['unshaded_flux_Wb']=abs(line((-29.999,50.25),(-21.001,50.25),0))
        row['shaded_flux_Wb']=abs(line((-18.999,50.25),(-10.001,50.25),0))
        # Для суперпозиции нужны ориентированные значения, не модули.
        row['signed_unshaded_flux_Wb']=line((-29.999,50.25),(-21.001,50.25),0)
        row['signed_shaded_flux_Wb']=line((-18.999,50.25),(-10.001,50.25),0)
    return finish(stem,path,row)
def conjugate(x): return complex(x).conjugate()
def real_force(row):
    return sum((p[1]**2-p[0]**2) for p in row['gap_B_T'])*.030*.0001/(2*MU0)
def acharmonic(state):
    i1=complex(*state['I1_RMS_A']); i2=complex(*state['I2_RMS_A'])
    acbuild(math.sqrt(2)*i1,math.sqrt(2)*i2,.7)
    mi_probdef(50,'millimeters','planar',1e-9,30,25,0)
    stem,path=begin_solve('AC-50Гц-комплексная-проверка')
    psi=[complex(mo_getcircuitproperties(x)[2])/math.sqrt(2) for x in [PRIMARY,SECONDARY]]
    mo_clearblock(); mo_groupselectblock(2); force=-real(mo_blockintegral(19)); mo_clearblock()
    row={'family':'AC-частичный-полюс-гармоника','tag':'50Гц-комплексная-проверка','frequency_Hz':50,'I_peak_A':[pair(math.sqrt(2)*i1),pair(math.sqrt(2)*i2)],'linkages_RMS_Wb':[pair(p) for p in psi],'force_mean_WST_N':force,'scope':'Токи RMS внешней цепи перенесены в FEMM как peak с множителем sqrt(2). Линейная сталь σ=0; медные R учтены во внешней цепи.'}
    return finish(stem,path,row)
def acstate(bases,r2,closed=True,loss_resistance=0):
    # Все комплексные токи и U ниже — RMS; L получены из DC полей.
    # e^{jωt}, пассивное правило U=RI+jωΨ, I2=-jωM I1/(R2+jωL22).
    a,b=bases; l11=a['linkages_Wb'][0]; l22=b['linkages_Wb'][1]
    m=(a['linkages_Wb'][1]+b['linkages_Wb'][0])/2
    omega=2*math.pi*50; r1=8; voltage=24
    ratio=-1j*omega*m/(r2+1j*omega*l22) if closed else 0j
    z=r1+loss_resistance+1j*omega*l11+1j*omega*m*ratio
    i1=voltage/z; i2=ratio*i1
    psi1=l11*i1+m*i2; psi2=m*i1+l22*i2
    pinput=(voltage*conjugate(i1)).real
    pcu1=r1*abs(i1)**2; pcu2=r2*abs(i2)**2
    # Представление постоянных потерь стали последовательным Rfe
    # допустимо лишь в локальной линейной рабочей точке. Значение
    # назначено для обучения, не вычислено из DC B(H).
    pfe=loss_resistance*abs(i1)**2
    force0=0; force2=0j; dx=.120/1200; depth=.030
    for k in range(1200):
        bx=a['gap_B_T'][k][0]*i1+b['gap_B_T'][k][0]*i2
        by=a['gap_B_T'][k][1]*i1+b['gap_B_T'][k][1]*i2
        force0+=(abs(by)**2-abs(bx)**2)*depth*dx/(2*MU0)
        force2+=(by*by-bx*bx)*depth*dx/(2*MU0)
    fu=a['signed_unshaded_flux_Wb']*i1+b['signed_unshaded_flux_Wb']*i2
    fs=a['signed_shaded_flux_Wb']*i1+b['signed_shaded_flux_Wb']*i2
    angle=math.atan2((fs*conjugate(fu)).imag,(fs*conjugate(fu)).real)*180/math.pi
    # RMS → peak: sqrt(2) множитель; F0 и F2 уже учитывают его.
    primary_residual=voltage-(r1+loss_resistance)*i1-1j*omega*psi1
    secondary_residual=r2*i2+1j*omega*psi2 if closed else i2
    return {'closed_secondary':closed,'voltage_RMS_V':voltage,'frequency_Hz':50,'R1_ohm':r1,'R2_ohm':r2,'Rfe_local_ohm':loss_resistance,'I1_RMS_A':pair(i1),'I2_RMS_A':pair(i2),'I1_RMS_abs_A':abs(i1),'I2_RMS_abs_A':abs(i2),'impedance_ohm':pair(z),'Pinput_W':pinput,'Pcu1_W':pcu1,'Pcu2_W':pcu2,'Pfe_prescribed_W':pfe,'power_balance_relative_error':rel(pinput,pcu1+pcu2+pfe),'primary_KVL_relative_error':abs(primary_residual)/voltage,'secondary_KVL_V':abs(secondary_residual),'unshaded_flux_RMS_Wb':pair(fu),'shaded_flux_RMS_Wb':pair(fs),'shaded_minus_unshaded_phase_deg':angle,'force_mean_N':force0,'force_2f_amplitude_N':abs(force2),'force_min_N':force0-abs(force2),'force_max_N':force0+abs(force2),'force_scope':'Двумерный тензор Maxwell по воздушной линии 120 мм в середине зазора; длинный средний участок 30 мм, линейный материал. Торцевые поля, насыщение, механика, пружина и движение отсутствуют.','phasor_convention':'RMS, e^(jωt), физическая величина sqrt(2)Re(X exp(jωt)); FEMM базис DC, не AC peak.'}
def accampaign():
    a=acsolve('AC-базис-основная',1,0,.7); b=acsolve('AC-базис-виток',0,1,.7)
    fine=acsolve('AC-сетка-основная',1,0,.5)
    l11=a['linkages_Wb'][0]; l22=b['linkages_Wb'][1]; m=(a['linkages_Wb'][1]+b['linkages_Wb'][0])/2
    check('AC-reciprocity',rel(a['linkages_Wb'][1],b['linkages_Wb'][0])<.005,rel(a['linkages_Wb'][1],b['linkages_Wb'][0]),.005,'Две независимые единичные возбуждающие цепи.')
    check('AC-positive-energy',l11>0 and l22>0 and l11*l22-m*m>0,[l11,l22,m,l11*l22-m*m],0,'Положительная определённость L, пассивность без потерь.')
    check('AC-field-mesh',rel(l11,fine['linkages_Wb'][0])<.01,rel(l11,fine['linkages_Wb'][0]),.01,'Сетка 0.7 → 0.5 по L11.')
    check('AC-force-mesh',rel(real_force(a),real_force(fine))<.01,rel(real_force(a),real_force(fine)),.01,'Сгущение сетки по независимому интегралу тензора в зазоре.')
    check('AC-force-methods',rel(real_force(a),a['armature_attraction_WST_N'])<.01,rel(real_force(a),a['armature_attraction_WST_N']),.01,'Линейное поле: тензор по линии зазора против взвешенного Maxwell всего якоря.')
    # Прямоугольная медная петля: две активные стороны по30 мм и
    # две торцевые перемычки по12 мм; площадь сечения 2*4 мм².
    # σ=58 МСм/м при принятой температуре; skin/proximity опущены.
    r2=(2*.030+2*.012)/(58e6*8e-6)
    states=[acstate([a,b],r2,False,0),acstate([a,b],r2,True,0),acstate([a,b],r2,True,1.2)]
    for j,s in enumerate(states):
        check('AC-power-%d'%j,s['power_balance_relative_error']<1e-10,s['power_balance_relative_error'],1e-10,'U RMS, медь основной/вторичной цепи и отдельно заданные потери стали.')
        check('AC-KVL-%d'%j,s['primary_KVL_relative_error']<1e-10 and s['secondary_KVL_V']<1e-10,[s['primary_KVL_relative_error'],s['secondary_KVL_V']],1e-10,'Пассивная вторичная цепь без источника напряжения.')
    check('AC-partial-phase',abs(states[1]['shaded_minus_unshaded_phase_deg'])>8,states[1]['shaded_minus_unshaded_phase_deg'],8,'Фазовое разделение потоков именно двух частей одного полюса.')
    check('AC-force-minimum',states[1]['force_min_N']>states[0]['force_min_N']+1e-4,[states[0]['force_min_N'],states[1]['force_min_N']],1e-4,'Минимум мгновенной силы рассчитан по полю, а не только средняя сила.')
    # Третье независимое поле проверяет суперпозицию и энергию.
    mixed=acsolve('AC-проверка-суперпозиции',.7,-50,.7)
    expected=[.7*l11-50*m,.7*m-50*l22]
    check('AC-superposition',max(rel(mixed['linkages_Wb'][i],expected[i]) for i in [0,1])<.006,max(rel(mixed['linkages_Wb'][i],expected[i]) for i in [0,1]),.006,'Реальные линейные токи .7 и -50 А, отдельное FEMM поле.')
    energy=.5*(.7*.7*l11+2*.7*(-50)*m+2500*l22)
    check('AC-field-energy',rel(mixed['energy_J'],energy)<.006,rel(mixed['energy_J'],energy),.006,'Интеграл энергии FEMM против 1/2 i^T L i.')
    harmonic=acharmonic(states[1])
    check('AC-harmonic-force',rel(harmonic['force_mean_WST_N'],states[1]['force_mean_N'])<.015,rel(harmonic['force_mean_WST_N'],states[1]['force_mean_N']),.015,'Нативный комплексный solve 50 Гц с peak-токами против цепного RMS и интеграла поля.')
    i1=complex(*states[1]['I1_RMS_A']); i2=complex(*states[1]['I2_RMS_A'])
    expected_psi=[l11*i1+m*i2,m*i1+l22*i2]
    e=max(rel(complex(*harmonic['linkages_RMS_Wb'][i]),expected_psi[i]) for i in [0,1])
    check('AC-harmonic-linkage',e<.006,e,.006,'Нативное гармоническое FEMM поле и потокосцепления матрицы L/M; peak → RMS.')
    save('ac-results.json',{'basis_models':[a['model'],b['model']],'harmonic_model':harmonic['model'],'L_H':[[l11,m],[m,l22]],'states':states,'R2_origin':'R2=l/(σS): l=2*30+2*12=84 мм, S=2*4=8 мм², σ=58 МСм/м; обе отсутствующие в поле торцевые перемычки входят в R2. Температура фиксирована, skin/proximity опущены.','steel_loss_origin':'Rfe=1.2 Ω — заданный локальный эквивалент активных потерь; из линейного DC поля не выведен. Сравнение без потерь/с потерями явно раздельно.','scope':'Плоский средний участок; обе активные стороны витка и замыкание вне плоскости объявлены. Не физическая 3D аттестация контактора.'})

save('geometry-contract.json',{'source':'О.Б.Буль, 2005: П6 с.319–324; гл.2 с.44–85, частичный КЗВ с.78–85. Размеры и материалы адаптации самостоятельные.','P6':{'problem':'axi, половина z>=0','full_gap_mm':1,'pole_radius_mm':8,'pole_z_mm':[.5,20],'shank_radius_mm':15,'yoke_r_mm':[45,55],'coil_rz_mm':[17,22,35,38],'turns_half':150,'full_mmf':'300I','section_count_finest':64,'group_names':{'1':'Полюс и непрерывное ярмо','3':'Катушка','4':'Воздух'}},'AC':{'problem':'planar','depth_mm':30,'gap_each_mm':.5,'primary_turns':200,'secondary_turns':1,'secondary_cross_sections_mm':[[-21,46,-19,50],[-9,46,-7,50]],'closure':'Торцевые перемычки вне плоскости; их сопротивление входит в полное R2.','group_names':{'1':'Магнитопровод','2':'Якорь','3':'Основная обмотка','4':'Частичный виток','5':'Воздух'}},'BH_T_A_per_m':BH,'BH_origin':'Условная монотонная кривая; не марка стали и не измеренная характеристика.'})
template=p6campaign()
if MODE=='full': p6campaign(1.45,True,template)
if MODE!='p6': accampaign()
for r in RECORDS: check('reopen-'+r['model'],r['russian_circuit_reopened'],r['russian_circuit_reopened'],True,'Сохранение, повторное открытие, решение и повторное чтение русской цепи FEMM 5.')
save('result.json',{'schema':'bul-completion-circuits-v1','implementation':'FEMM 5 embedded MicroPython: all fields, P6 recurrences, external AC circuit, force and energy arithmetic.','mode':MODE,'records':RECORDS,'checks':CHECKS,'passed':all(x['passed'] for x in CHECKS)})
print('COMPLETED',len(RECORDS),'models;',len(CHECKS),'checks;',all(x['passed'] for x in CHECKS))