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

Запустите из отдельной рабочей папки. Создаются calculated-NN, модели,
решения и results.json. Координаты в мм, активная длина плоских задач 1 м.
Комплексные токи и поля — амплитудные фазоры; P — средняя мощность.
"""
import os, sys, math, json
MU0=4e-7*math.pi
def sinh(x):return (math.exp(x)-math.exp(-x))/2
def cosh(x):return (math.exp(x)+math.exp(-x))/2
SCRIPT=globals().get('__file__','')
BASE=os.path.dirname(os.path.abspath(SCRIPT)) if SCRIPT else os.getcwd()
if globals().get('SC_OUTPUT'):
    OUT=SC_OUTPUT
else:
    k=1
    while os.path.exists(BASE+'/calculated-%02d'%k):k+=1
    OUT=BASE+'/calculated-%02d'%k
if not os.path.exists(OUT):os.mkdir(OUT)
MODE=globals().get('SC_MODE','all')
ROWS=[]
CHECKS=[]
NODES={};SEGMENTS={}

def real(z):return float(complex(z).real)
def pair(z):return [complex(z).real,complex(z).imag]
def relative(a,b):return abs(a-b)/max(abs(b),1e-20)
def save():
    with open(OUT+'/results.json','w') as f:json.dump({'rows':ROWS,'checks':CHECKS,'passed':all(c['passed'] for c in CHECKS),'runtime':sys.version},f)
def check(name,value,limit):
    CHECKS.append({'name':name,'value':value,'limit':limit,'passed':value<=limit})
    save()
def node(x,y):
    key=(round(x,9),round(y,9))
    if key not in NODES:mi_addnode(x,y);NODES[key]=1
def segment(a,b,h,boundary='<None>',group=0):
    node(a[0],a[1]);node(b[0],b[1])
    key=tuple(sorted((a,b)))
    if key not in SEGMENTS:mi_addsegment(a[0],a[1],b[0],b[1]);SEGMENTS[key]=1
    mi_selectsegment((a[0]+b[0])/2,(a[1]+b[1])/2)
    mi_setsegmentprop(boundary,h,0,0,group);mi_clearselected()
def polygon(points,h,boundary='<None>',group=0):
    for i in range(len(points)):segment(points[i],points[(i+1)%len(points)],h,boundary,group)
def rectangle(x1,y1,x2,y2,h,boundary='<None>',group=0):polygon([(x1,y1),(x2,y1),(x2,y2),(x1,y2)],h,boundary,group)
def label(x,y,material,h,group=0,circuit='<None>',turns=0):
    mi_addblocklabel(x,y);mi_selectlabel(x,y)
    mi_setblockprop(material,0,h,circuit,0,group,turns);mi_clearselected()
def bound(name,value=0):mi_addboundprop(name,value,0,0,0,0,0,0,0,0)
def start(f=0,axi=False):
    global NODES,SEGMENTS
    NODES={};SEGMENTS={}
    newdocument(0);mi_probdef(f,'millimeters','axi' if axi else 'planar',1e-9,1000,25,0)
    mi_setcomment('Учебная самостоятельная геометрия. Русские материалы и цепи. Плоская глубина 1000 мм; осесимметричная задача — полный оборот. Группа 1: воздух; 2: сталь; 10–99: проводники.')
    mi_addmaterial('Воздух',1,1,0,0,0)
def solve(case,lab,parameters):
    path=OUT+'/'+case+'.fem'
    mi_saveas(path);mi_close();opendocument(path)
    mi_analyze(1);mi_loadsolution();mo_smooth('off')
    row={'id':case,'lab_id':lab,'parameters':parameters,'model':case+'.fem','solution':case+'.ans','roundtrip':True,'nodes':mo_numnodes(),'elements':mo_numelements()}
    return row
def finish(row):
    ROWS.append(row);save();mo_close();mi_close();print('SC_DONE',row['id'],row['elements'])
def integral(group,number):
    mo_clearblock();mo_groupselectblock(group);v=mo_blockintegral(number);mo_clearblock();return v
def optional_integral(group,number):
    try:return integral(group,number)
    except Exception as e:
        if 'No area selected' in str(e):return 0
        raise
def total(number):
    mo_clearblock();mo_groupselectblock();v=mo_blockintegral(number);mo_clearblock();return v
def circuit(name):return [pair(z) for z in mo_getcircuitproperties(name)]

def slot_domain(width=8,depth=20,gap=2,pitch=40,h=.3,B=.8):
    p=pitch/2;b=width/2
    points=[(-p,0),(-b,0),(-b,-depth),(b,-depth),(b,0),(p,0),(p,gap),(-p,gap)] if width else [(-p,0),(p,0),(p,gap),(-p,gap)]
    bound('Левый заданный поток',B*p*.001);bound('Правый заданный поток',-B*p*.001)
    polygon(points,h)
    segment((-p,0),(-p,gap),h,'Левый заданный поток')
    segment((p,0),(p,gap),h,'Правый заданный поток')

def gap_carter():
    lab='slots-gap-carter'
    for width in [0,4,8,12]:
        for h in ([.6,.3,.15] if width in [0,8] else [.3]):
            start();slot_domain(width=width,h=h);label(0,1,'Воздух',h,1)
            row=solve('sc-carter-b%02d-h%03d'%(width,int(h*100)),lab,{'opening_mm':width,'pitch_mm':40,'gap_mm':2,'slot_depth_mm':20,'mesh_mm':h,'Bmean_T':.8})
            W=real(total(2));delta=2*MU0*W/(.8**2*.04)
            row['metrics']={'energy_J':W,'effective_gap_mm':delta*1000,'carter':delta/.002,'tooth_mid_B_T':real(mo_getpointvalues(14,1)[2])}
            rho=width/4
            decrement=8/math.pi*(rho*math.atan(rho)-.5*math.log(1+rho*rho))
            row['metrics']['carter_conformal']=40/(40-decrement)
            row['profile']=[{'x_mm':-19.95+k*39.9/159,'By_T':real(mo_getpointvalues(-19.95+k*39.9/159,1)[2])} for k in range(160)]
            finish(row)
            check('Картер: конформная оценка '+row['id'],relative(row['metrics']['carter'],row['metrics']['carter_conformal']),.01)
    fine=[r for r in ROWS if r['id']=='sc-carter-b08-h015'][0]
    medium=[r for r in ROWS if r['id']=='sc-carter-b08-h030'][0]
    smooth=[r for r in ROWS if r['id']=='sc-carter-b00-h015'][0]
    check('Плоский гладкий зазор: энергия',relative(smooth['metrics']['carter'],1),.003)
    check('Коэффициент Картера: сгущение 0.30→0.15 мм',relative(medium['metrics']['carter'],fine['metrics']['carter']),.015)
    # Возбуждение одинаковым пространственным токовым слоем; собственные токи
    # интегрируются по ширине ячеек. Нулевая проводимость исключает вихревые токи.
    for width in [0,8]:
        for count in [24,48]:
            h=.22;pitch=40;gap=2;thick=.3
            start(50)
            slot_domain(width=width,h=h,B=0)
            # Dirichlet на сечениях синусоидальной симметрии, cos ток и sin B.
            for n in range(count):
                x1=-20+40*n/count;x2=-20+40*(n+1)/count;xm=(x1+x2)/2
                j=complex(math.cos(2*math.pi*xm/pitch),0)+complex(.3*math.cos(6*math.pi*xm/pitch)*math.cos(.6),.3*math.cos(6*math.pi*xm/pitch)*math.sin(.6))
                name='Токовый слой %02d'%n
                mi_addmaterial(name,1,1,0,10*j,0)
                rectangle(x1,gap-thick,x2,gap,h)
                label(xm,gap-thick/2,name,h,10)
            label(0,.7,'Воздух',h,1)
            row=solve('sc-harmonic-b%02d-n%02d'%(width,count),lab,{'opening_mm':width,'segments':count,'gap_mm':gap,'source_thickness_mm':thick,'J1_peak_A_per_mm2':10,'J3_ratio':.3,'J3_phase_rad':.6,'sample_y_mm':1,'frequency_Hz':50,'conductivity_S_per_m':0})
            samples=[];spec=[];ns=256
            for k in range(ns):
                x=-20+(k+.5)*40/ns;by=complex(mo_getpointvalues(x,1)[2]);samples.append([x,pair(by)])
            for nu in [1,3,5,7]:
                z=sum(complex(v[1][0],v[1][1])*complex(math.cos(-2*math.pi*nu*v[0]/40),math.sin(-2*math.pi*nu*v[0]/40)) for v in samples)*2/ns
                spec.append({'order':nu,'By_complex_T':pair(z),'amplitude_T':abs(z),'phase_rad':math.atan2(z.imag,z.real)})
            row['harmonics']=spec;row['profile']=samples
            k1=2*math.pi/.04
            row['metrics']={'B1_T':spec[0]['amplitude_T'],'B3_T':spec[1]['amplitude_T'],'analytic_B1_T':MU0*1e7*sinh(k1*.0003)*cosh(k1*.001)/(k1*sinh(k1*.002))}
            finish(row)

    sm=[r for r in ROWS if r['id']=='sc-harmonic-b00-n48'][0]
    check('Гармонический гладкий зазор: решение Лапласа',relative(sm['metrics']['B1_T'],sm['metrics']['analytic_B1_T']),.025)
    for w in [0,8]:
        a=[r for r in ROWS if r['id']=='sc-harmonic-b%02d-n24'%w][0];b=[r for r in ROWS if r['id']=='sc-harmonic-b%02d-n48'%w][0]
        check('Дробление токового слоя b=%d'%w,relative(a['metrics']['B1_T'],b['metrics']['B1_T']),.025)

def circle(radius,h,boundary='<None>',group=0):
    for p in [(radius,0),(0,radius),(-radius,0),(0,-radius)]:node(p[0],p[1])
    for a,b in [((radius,0),(0,radius)),((0,radius),(-radius,0)),((-radius,0),(0,-radius)),((0,-radius),(radius,0))]:
        mi_addarc(a[0],a[1],b[0],b[1],90,3)
        angle=(math.atan2(a[1],a[0])+math.pi/4)
        mi_selectarcsegment(radius*math.cos(angle),radius*math.sin(angle));mi_setarcsegmentprop(3,boundary,0,group);mi_clearselected()

def stirrer():
    # Двухполюсное вращающееся поле: фазоры J(theta)=Jm exp(-j theta).
    # Длинный проводящий цилиндр неподвижен; движение жидкости здесь не задаётся.
    configurations=[(f,.35,48,160,'') for f in [10,50,100,200,225,250,400,800]]+[(400,.175,48,160,'-fine'),(400,.35,96,160,'-source96'),(400,.35,48,320,'-domain')]
    for f,hc,count,rout,suffix in configurations:
        start(f);bound('Внешняя граница A=0');circle(rout,3,'Внешняя граница A=0')
        mi_addmaterial('Проводящий цилиндр 20.5 МСм на метр',1,1,0,0,20.5);circle(10,hc,group=10);label(0,0,'Проводящий цилиндр 20.5 МСм на метр',hc,10)
        J=1000/(2*19.5)
        for k in range(count):
            t1=2*math.pi*k/count;t2=2*math.pi*(k+1)/count;tm=(t1+t2)/2
            name='Фазированный токовый слой %02d'%(k+1);mi_addmaterial(name,1,1,0,J*complex(math.cos(tm),-math.sin(tm)),0)
            pts=[(r*math.cos(t),r*math.sin(t)) for r,t in [(19,t1),(20,t1),(20,t2),(19,t2)]]
            polygon(pts,.4);label(19.5*math.cos(tm),19.5*math.sin(tm),name,.4,20)
        label(15,0,'Воздух',.5,1);label(35,0,'Воздух',3 if rout==160 else 6,1)
        row=solve('sc-rotating-f%03d'%f+suffix,'conductors-ac-transposition',{'frequency_Hz':f,'pole_pairs':1,'cylinder_radius_mm':10,'conductivity_MS_per_m':20.5,'source_ring_radii_mm':[19,20],'source_segments':count,'source_J_peak_A_per_mm2':J,'mechanical_speed_rpm':0,'active_length_m':1,'outer_radius_mm':rout,'mesh_conductor_mm':hc})
        P=real(integral(10,4));T=real(integral(10,15));Tw=real(integral(10,22))
        row['metrics']={'loss_W_per_m':P,'Lorentz_torque_Nm_per_m':T,'WST_torque_Nm_per_m':Tw,'synchronous_speed_rpm':60*f,'power_torque_balance_relative':abs(P-2*math.pi*f*abs(T))/max(P,1e-20)}
        check('Момент Лоренца и WST f=%d'%f,relative(T,Tw),.025)
        check('Вращающееся поле: P=omega_sync T f=%d'%f,row['metrics']['power_torque_balance_relative'],.025)
        finish(row)
    base=[r for r in ROWS if r['id']=='sc-rotating-f400'][0]
    for suffix,limit in [('-fine',.015),('-source96',.015),('-domain',.04)]:
        refined=[r for r in ROWS if r['id']=='sc-rotating-f400'+suffix][0]
        check('Вращающееся поле: независимое уточнение '+suffix,relative(base['metrics']['Lorentz_torque_Nm_per_m'],refined['metrics']['Lorentz_torque_Nm_per_m']),limit)

def permeance():
    lab='slots-permeance-matrix'
    # T-область с идеальными железными поверхностями. Токовый слой внизу
    # конечной толщины t; энергия равномерной части соответствует h-2t/3.
    for depth,gap in [(d,g) for g in [5,20] for d in [10,20,30,50]]:
        for h in ([.6,.3,.15] if depth==20 and gap==5 else [.3]):
            start();slot_domain(width=20,depth=depth,gap=gap,pitch=200,h=h,B=0)
            mi_addmaterial('Обмотка токового слоя',1,1,0,0,0);mi_addcircprop('Ток паза',1000,1)
            rectangle(-10,-depth,10,-depth+.5,h)
            label(0,-depth+.25,'Обмотка токового слоя',min(h,.15),10,'Ток паза',1)
            label(0,2,'Воздух',h,1)
            row=solve('sc-toothtip-g%02d-d%02d-h%03d'%(gap,depth,int(h*100)),lab,{'slot_depth_mm':depth,'slot_width_mm':20,'gap_mm':gap,'span_mm':200,'layer_thickness_mm':.5,'current_A':1000,'mesh_mm':h,'reference_depth_mm':depth-1})
            W=real(total(2));lam=2*W/(MU0*1e6);L=2*W/1e6
            href=depth-1
            aleak=real(mo_getpointvalues(0,-href)[0])-real(mo_getpointvalues(0,gap)[0])
            lambda_leak=aleak/(MU0*1000)
            rho=20/(2*gap)
            lv=(math.log((gap/20)**2+.25)+4*gap/20*math.atan(rho))/(2*math.pi)
            li=1/(2*rho)-(math.log(2/math.sqrt(1+1/rho**2))+math.atan(1/rho)/rho)/math.pi
            row['metrics']={'energy_J':W,'total_permeance_H':L,'lambda_total':lam,'lambda_leakage':lambda_leak,'lambda_neck':href/20,'lambda_tip_signed':lambda_leak-href/20,'formula_Voldek':lv,'formula_Ivanov_Smolensky_beta1':li,'leakage_flux_Wb':aleak,'circuit_linkage_Wb':real(mo_getcircuitproperties('Ток паза')[2])}
            check('Энергия и потокосцепление '+row['id'],relative(L,real(mo_getcircuitproperties('Ток паза')[2])/1000),.003)
            finish(row)
    a=[r for r in ROWS if r['id']=='sc-toothtip-g05-d20-h030'][0];b=[r for r in ROWS if r['id']=='sc-toothtip-g05-d20-h015'][0]
    check('Проводимость коронки: сгущение',abs(a['metrics']['lambda_tip_signed']-b['metrics']['lambda_tip_signed']),.008)
    for gap in [1,2,4]:
        for excitation in ['A','B']:
            start();mi_addmaterial('Сталь линейная мю 10000',10000,10000,0,0,0);mi_addmaterial('Обмотка без вихревых токов',1,1,0,0,0)
            mi_addcircprop('Катушка A',1 if excitation=='A' else 0,1);mi_addcircprop('Катушка B',1 if excitation=='B' else 0,1)
            bound('Удалённая граница A=0');rectangle(-65,-55,65,25,3,'Удалённая граница A=0')
            # Зубцы отдельны от ярма лишь геометрическими линиями; один материал.
            rectangle(-45,-40,45,-30,2);label(0,-35,'Сталь линейная мю 10000',2,2)
            for cx in [-15,15]:rectangle(cx-5,-30,cx+5,-gap,.5);label(cx,-15,'Сталь линейная мю 10000',.5,2)
            rectangle(-45,0,45,8,1);label(0,4,'Сталь линейная мю 10000',1,2)
            for index,cx in enumerate([-15,15]):
                for sign in [-1,1]:
                    xc=cx+sign*8;rectangle(xc-1,-27,xc+1,-24,.5)
                    label(xc,-25.5,'Обмотка без вихревых токов',.4,10+index*2+(sign+1)//2,'Катушка '+('A' if index==0 else 'B'),sign)
            label(55,10,'Воздух',2,1)
            row=solve('sc-matrix-g%02d-%s'%(gap,excitation),lab,{'gap_mm':gap,'excitation':excitation,'current_A':1,'turns':1,'mu_r':10000})
            psi=[real(mo_getcircuitproperties('Катушка '+p)[2]) for p in ['A','B']]
            weighted=[]
            for idx in [0,1]:
                means=[real(integral(10+idx*2+j,1))/real(integral(10+idx*2+j,10)) for j in [0,1]]
                weighted.append(means[1]-means[0])
            row['metrics']={'psi_A_Wb':psi[0],'psi_B_Wb':psi[1],'energy_J':real(total(2)),'area_weighted_A_linkage_A_Wb':weighted[0],'area_weighted_A_linkage_B_Wb':weighted[1]}
            check('Согласование среднего A и потокосцепления '+row['id'],max(abs(weighted[j]-psi[j]) for j in [0,1])/max(abs(z) for z in psi),.002)
            finish(row)
    for g in [1,2,4]:
        a=[r for r in ROWS if r['id']=='sc-matrix-g%02d-A'%g][0];b=[r for r in ROWS if r['id']=='sc-matrix-g%02d-B'%g][0]
        check('Взаимность L12=L21, g=%d'%g,relative(a['metrics']['psi_B_Wb'],b['metrics']['psi_A_Wb']),.002)
    mutual_position()

def wrap(x):return (x+15)%30-15

def mutual_position():
    # Один период двух зубчатых сердечников; равные встречные токи в каждом
    # зубцовом контуре. Проводимость определяется двумя независимыми МДС.
    for shift in range(0,61,3):
        for excitation in ['A','B']:
            start();mi_addmaterial('Сталь линейная мю 10000',10000,10000,0,0,0)
            mats=[]
            for g in range(4):
                amp=.25*(1 if g%2 else -1) if (g<2 and excitation=='A') or (g>=2 and excitation=='B') else 0
                name='Сторона зубцового контура '+str(g+1);mi_addmaterial(name,1,1,0,amp,0);mats.append(name)
            xset=[-15,15,-9,-7,-5,5,7,9]
            for x in [-9,-7,-5,5,7,9]:xset.append(wrap(x+shift))
            xs=sorted(set(round(x,8) for x in xset));ys=[-40,-30,-26,-24,-2,0,4,6,10,20]
            bound('Внешняя поверхность ярма A=0')
            for iy in range(len(ys)-1):
                y=(ys[iy]+ys[iy+1])/2
                pname='Антипериодические стороны '+str(iy+1)
                mi_addboundprop(pname,0,0,0,0,0,0,0,0,5)
                for ix in range(len(xs)-1):
                    x=(xs[ix]+xs[ix+1])/2;u=wrap(x-shift)
                    rectangle(xs[ix],ys[iy],xs[ix+1],ys[iy+1],.6)
                    material='Воздух';group=1
                    if y<-30 or (-30<y<-2 and abs(x)<5) or y>10 or (0<y<10 and abs(u)<5):material='Сталь линейная мю 10000';group=2
                    if -26<y<-24 and 7<abs(x)<9:group=10+(1 if x>0 else 0);material=mats[group-10]
                    if 4<y<6 and 7<abs(u)<9:
                        image=int(round((x-shift-u)/30))
                        physical_sign=(1 if u>0 else -1)*(-1 if image%2 else 1)
                        group=12+(1 if physical_sign>0 else 0);material=mats[group-10]
                    label(x,y,material,.5 if group!=1 else .65,group)
                segment((-15,ys[iy]),(-15,ys[iy+1]),.5,pname)
                segment((15,ys[iy]),(15,ys[iy+1]),.5,pname)
            for ix in range(len(xs)-1):
                segment((xs[ix],-40),(xs[ix+1],-40),.6,'Внешняя поверхность ярма A=0')
                segment((xs[ix],20),(xs[ix+1],20),.6,'Внешняя поверхность ярма A=0')
            row=solve('sc-mutual-x%02d-%s'%(shift,excitation),'slots-permeance-matrix',{'displacement_mm':shift,'tooth_pitch_mm':30,'magnetic_period_mm':60,'gap_mm':2,'tooth_width_mm':10,'excitation':excitation,'current_A':1,'mu_r':10000,'boundary':'anti-periodic: neighbouring tooth contour has opposite MMF'})
            psiA=.25e6*real(optional_integral(11,1)-optional_integral(10,1));psiB=.25e6*real(optional_integral(13,1)-optional_integral(12,1))
            row['metrics']={'psi_A_Wb':psiA,'psi_B_Wb':psiB,'energy_J':real(total(2))}
            finish(row)
    for x in range(0,61,3):
        a=[r for r in ROWS if r['id']=='sc-mutual-x%02d-A'%x][0];b=[r for r in ROWS if r['id']=='sc-mutual-x%02d-B'%x][0]
        check('Взаимность противоположных контуров x=%d'%x,abs(a['metrics']['psi_B_Wb']-b['metrics']['psi_A_Wb'])/max(abs(a['metrics']['psi_A_Wb']),abs(b['metrics']['psi_B_Wb'])),.002)
    a=[r for r in ROWS if r['id']=='sc-mutual-x00-A'][0];b=[r for r in ROWS if r['id']=='sc-mutual-x60-A'][0]
    check('Период взаимной проводимости 60 мм',relative(a['metrics']['psi_B_Wb'],b['metrics']['psi_B_Wb']),.002)

def stress_x(a,b,n,offset,count=180):
    dx=b[0]-a[0];dy=b[1]-a[1];length=math.sqrt(dx*dx+dy*dy)*.001;answer=0
    for k in range(count):
        t=(k+.5)/count;x=a[0]+dx*t+n[0]*offset;y=a[1]+dy*t+n[1]*offset
        p=mo_getpointvalues(x,y);bx=real(p[1]);by=real(p[2])
        answer+=((bx*bx-by*by)*n[0]/2+bx*by*n[1])/MU0
    return answer*length/count

def forces():
    lab='slots-force-partition'
    for d,h in [(d,.12) for d in [.7,2,3.5,3.6,4,5,6,7,12,17]]+[(7,.24),(7,.06)]:
        for current in [-1000,1000]:
            start();slot_domain(h=h)
            mi_addmaterial('Медный проводник без вытеснения тока',1,1,0,0,0);mi_addcircprop('Ток проводника',current,1)
            rectangle(-3,-d-.3,3,-d+.3,.12);label(0,-d,'Медный проводник без вытеснения тока',.12,10,'Ток проводника',1);label(0,1,'Воздух',h,1)
            row=solve('sc-force-d%03d-h%03d-%s'%(int(d*10),int(h*100),'plus' if current>0 else 'minus'),lab,{'center_depth_mm':d,'current_A':current,'Bmean_T':.8,'mesh_mm':h,'copper_size_mm':[6,.6]})
            mo_smooth('on')
            cu=real(integral(10,11))
            # Боковые сечения входят в полный контур. На них задан A, а не нулевая тяга.
            whole=stress_x((-20,1),(20,1),(0,1),0,700)+stress_x((-20,0),(-20,1),(-1,0),-.005,250)+stress_x((20,0),(20,1),(1,0),-.005,250)
            # Непрерывная ступенчатая линия лежит в воздухе и обходит медь;
            # стыки около углов включены в контур без разрывов.
            e=.35;x=4-e;y=-20+e
            wall=stress_x((-20,e),(-x,e),(0,1),0,350)+stress_x((-x,e),(-x,y),(1,0),0,500)+stress_x((-x,y),(x,y),(0,1),0,350)+stress_x((x,y),(x,e),(-1,0),0,500)+stress_x((x,e),(20,e),(0,1),0,350)+stress_x((-20,0),(-20,e),(-1,0),-.005,200)+stress_x((20,0),(20,e),(1,0),-.005,200)
            row['metrics']={'copper_Fx_N':cu,'assembly_stress_Fx_N':whole,'wall_stress_Fx_N':wall,'tooth_mid_B_T':real(mo_getpointvalues(14,.15)[2])}
            row['field_samples']=[{'point_mm':[x,y],'A_Wb_per_m':real(mo_getpointvalues(x,y)[0]),'Bx_T':real(mo_getpointvalues(x,y)[1]),'By_T':real(mo_getpointvalues(x,y)[2])} for x,y in [(0,-7),(0,-1),(0,1),(14,1)]]
            finish(row)
    for d,h in [(d,.12) for d in [.7,2,3.5,3.6,4,5,6,7,12,17]]+[(7,.24),(7,.06)]:
        prefix='sc-force-d%03d-h%03d-'%(int(d*10),int(h*100))
        p=[r for r in ROWS if r['id']==prefix+'plus'][0];m=[r for r in ROWS if r['id']==prefix+'minus'][0]
        cross={key:(p['metrics'][key]-m['metrics'][key])/2 for key in ['copper_Fx_N','assembly_stress_Fx_N','wall_stress_Fx_N']}
        cross['copper_fraction']=abs(cross['copper_Fx_N']/cross['assembly_stress_Fx_N'])
        cross['analytic_assembly_Fx_N']=-1000*(p['metrics']['tooth_mid_B_T']+m['metrics']['tooth_mid_B_T'])/2
        p['interaction']=cross
        cross['force_balance_relative']=abs(cross['copper_Fx_N']+cross['wall_stress_Fx_N']-cross['assembly_stress_Fx_N'])/abs(cross['assembly_stress_Fx_N'])
        check('Независимый баланс медь+стенки '+prefix,cross['force_balance_relative'],.025)
        check('Тяга ячейки Bзубца I '+prefix,relative(cross['assembly_stress_Fx_N'],cross['analytic_assembly_Fx_N']),.025)
    a=[r for r in ROWS if r['id']=='sc-force-d070-h012-plus'][0];b=[r for r in ROWS if r['id']=='sc-force-d070-h006-plus'][0]
    check('Сила меди: сеточная сходимость',relative(a['interaction']['copper_Fx_N'],b['interaction']['copper_Fx_N']),.01)
    force_superposition()
    save()

def force_superposition():
    for kind,B,I in [('even',.8,0),('odd',0,1000)]:
        start();slot_domain(h=.12,B=B)
        mi_addmaterial('Медный проводник без вытеснения тока',1,1,0,0,0);mi_addcircprop('Ток проводника',I,1)
        rectangle(-3,-7.3,3,-6.7,.12);label(0,-7,'Медный проводник без вытеснения тока',.12,10,'Ток проводника',1);label(0,1,'Воздух',.12,1)
        row=solve('sc-force-'+kind,'slots-force-partition',{'center_depth_mm':7,'current_A':I,'Bmean_T':B,'mesh_mm':.12})
        mo_smooth('on');row['metrics']={'energy_J':real(total(2))}
        row['field_samples']=[{'point_mm':[x,y],'A_Wb_per_m':real(mo_getpointvalues(x,y)[0]),'Bx_T':real(mo_getpointvalues(x,y)[1]),'By_T':real(mo_getpointvalues(x,y)[2])} for x,y in [(0,-7),(0,-1),(0,1),(14,1)]]
        finish(row)
    even=[r for r in ROWS if r['id']=='sc-force-even'][0];odd=[r for r in ROWS if r['id']=='sc-force-odd'][0]
    for kind,sign in [('plus',1),('minus',-1)]:
        row=[r for r in ROWS if r['id']=='sc-force-d070-h012-'+kind][0]
        err=max(abs(row['field_samples'][k][key]-even['field_samples'][k][key]-sign*odd['field_samples'][k][key]) for k in range(4) for key in ['Bx_T','By_T'])
        check('Независимая суперпозиция чётного и нечётного полей '+kind,err,.00001)

def conductor_domain(h,slotwidth=8,depth=40,gap=5):
    slot_domain(width=slotwidth,depth=depth,gap=gap,pitch=60,h=h,B=0)

def conductors():
    lab='conductors-ac-transposition'
    for f in [0,50,400]:
        for h in ([.7,.35,.175] if f==50 else [.35]):
            start(f);conductor_domain(h)
            mi_addmaterial('Проводник 20.5 МСм на метр',1,1,0,0,20.5);mi_addcircprop('Стержень',1280,1)
            segment((-4,0),(4,0),h);label(0,-20,'Проводник 20.5 МСм на метр',h,10,'Стержень',1);label(0,2,'Воздух',h,1)
            row=solve('sc-emde-f%03d-h%03d'%(f,int(h*1000)),lab,{'frequency_Hz':f,'conductivity_MS_per_m':20.5,'width_mm':8,'height_mm':40,'current_peak_A':1280,'mesh_mm':h})
            props=mo_getcircuitproperties('Стержень');Z=complex(props[1])/complex(props[0]);P=real(integral(10,4));Rdc=1/(20.5e6*.008*.04)
            xi=.04*math.sqrt(math.pi*MU0*20.5e6*f) if f else 0
            kr=xi*(sinh(2*xi)+math.sin(2*xi))/(cosh(2*xi)-math.cos(2*xi)) if f else 1
            kx=3/(2*xi)*(sinh(2*xi)-math.sin(2*xi))/(cosh(2*xi)-math.cos(2*xi)) if f else 1
            row['metrics']={'loss_W_per_m':P,'Rdc_ohm_per_m':Rdc,'Rac_ohm_per_m':2*P/1280**2 if f else P/1280**2,'Rac_Rdc':(2 if f else 1)*P/(1280**2*Rdc),'Z_ohm_per_m':pair(Z),'Emde_kr':kr,'Emde_kx':kx,'skin_depth_mm':40/xi if f else 0,'internal_magnetic_energy_J':real(integral(10,2))}
            row['profile']=[{'depth_below_top_mm':d,'J_A_per_mm2':pair(mo_getpointvalues(0,-d)[7])} for d in [.1,.5,1,2,4,6,10,15,20,25,30,35,39,39.9]]
            finish(row)
    fine=[r for r in ROWS if r['id']=='sc-emde-f050-h175'][0];medium=[r for r in ROWS if r['id']=='sc-emde-f050-h350'][0]
    check('Стержень: сгущение сетки',relative(medium['metrics']['Rac_Rdc'],fine['metrics']['Rac_Rdc']),.01)
    check('Стержень: формула Эмде для активного сопротивления',relative(fine['metrics']['Rac_Rdc'],fine['metrics']['Emde_kr']),.03)
    dc=[r for r in ROWS if r['id']=='sc-emde-f000-h350'][0]
    for r in ROWS:
        if r['id'].startswith('sc-emde') and r['parameters']['frequency_Hz']:
            r['metrics']['internal_reactance_factor']=2*r['metrics']['internal_magnetic_energy_J']/dc['metrics']['internal_magnetic_energy_J']
    check('Стержень: формула Эмде для внутренней реактивности',relative(fine['metrics']['internal_reactance_factor'],fine['metrics']['Emde_kx']),.02)
    # Постоянны площадь меди 304 мм², пакет 8×42 мм, Ipeak=6000 А, T=120 °C.
    # σ(T)=58/(1+0.00393*(120-20)) МСм/м; изоляционные зазоры не проводят.
    sigma=58/(1+.00393*100)
    for n in [1,2,4,8,16,32,64]:
        for connection in ['parallel','transposed']:
            start(400);conductor_domain(.4,9,46,5)
            mi_addmaterial('Медь 120 градусов',1,1,0,0,sigma)
            ie=6000/(2*n) if connection=='transposed' else 6000
            mi_addcircprop('Параллельные ветви' if connection=='parallel' else 'Равные токи транспонированных жил',ie,0 if connection=='parallel' else 1)
            cname='Параллельные ветви' if connection=='parallel' else 'Равные токи транспонированных жил'
            for j in range(n):
                yc=-44+(j+.5)*42/n;hh=20/n
                for col in [-1,1]:
                    xc=col*2.1;rectangle(xc-1.9,yc-hh,xc+1.9,yc+hh,min(.35,hh/3))
                    label(xc,yc,'Медь 120 градусов',min(.35,hh/3),10,cname,1)
            label(0,2,'Воздух',.4,1)
            row=solve('sc-strands-n%02d-%s'%(n,connection),lab,{'frequency_Hz':400,'temperature_C':120,'conductivity_MS_per_m':sigma,'rows':n,'columns':2,'strand_count':2*n,'connection':connection,'circuit_current_peak_A':ie,'total_current_peak_A':6000,'copper_area_mm2':304,'pack_width_mm':8,'pack_height_mm':42,'strand_width_mm':3.8,'strand_height_mm':40/n})
            P=real(integral(10,4));Rdc=1/(sigma*1e6*304e-6);Rac=2*P/6000**2
            props=mo_getcircuitproperties(cname)
            row['metrics']={'loss_W_per_m':P,'Rac_Rdc':Rac/Rdc,'Rac_ohm_per_m':Rac,'Rdc_ohm_per_m':Rdc,'circuit_voltage_V':pair(props[1]),'circuit_current_A':pair(props[0]),'copper_total_current_A':pair(integral(10,7))}
            Zeq=complex(props[1])/(6000*(2*n if connection=='transposed' else 1))
            row['metrics']['equivalent_Z_ohm_per_m']=pair(Zeq)
            check('Мощность и действительная часть Z '+row['id'],relative(Zeq.real,Rac),.002)
            row['profile']=[{'x_mm':2.1,'y_mm':-44+(j+.5)*42/n,'J_A_per_mm2':pair(mo_getpointvalues(2.1,-44+(j+.5)*42/n)[7])} for j in range(n)]
            strand_currents=[]
            for j in range(n):
                yc=-44+(j+.5)*42/n
                for xc in [-2.1,2.1]:
                    mo_clearblock();mo_selectblock(xc,yc);strand_currents.append(pair(mo_blockintegral(7)));mo_clearblock()
            row['strand_currents_A']=strand_currents
            check('Суммарный ток меди '+row['id'],relative(complex(integral(10,7)),6000),.002)
            finish(row)
    stirrer()

def ring():
    lab='ring-ac-shield'
    for f in [10,50,200]:
        for h in ([1,.5,.25] if f==50 else ([.5,.25] if f==200 else [.5])):
            for sigma in ([0,3] if h==.5 else [3]):
                start(f,True);mi_addmaterial('Стальной замкнутый виток',500,500,0,0,sigma);mi_addmaterial('Обмотка возбуждения',1,1,0,0,0)
                mi_addcircprop('Катушка возбуждения',100,1)
                bound('Внешняя граница A=0');rectangle(0,-400,700,400,24,'Внешняя граница A=0')
                rectangle(100,-30,110,30,2);label(105,0,'Обмотка возбуждения',2,20,'Катушка возбуждения',200)
                rectangle(140,-10,160,10,h);label(150,0,'Стальной замкнутый виток',h,10)
                label(230,80,'Воздух',8,1)
                row=solve('sc-ring-f%03d-s%d-h%03d'%(f,sigma,int(h*100)),lab,{'frequency_Hz':f,'ring_mu_r':500,'ring_conductivity_MS_per_m':sigma,'ring_mean_radius_mm':150,'ring_cross_section_mm':[20,20],'excitation_current_peak_A':100,'turns':200,'mesh_ring_mm':h,'domain_mm':[0,-400,700,400]})
                # Поток через окружность r=130 мм: Phi=2*pi*r*Aphi.
                phi=2*math.pi*.13*complex(mo_getpointvalues(130,0)[0]);scale=.0015/abs(phi)
                it=complex(integral(10,7));P=real(integral(10,4));Rdc=2*math.pi*.15/(3e6*.02*.02)
                row['metrics']={'flux_130mm_Wb':pair(phi),'current_ring_peak_A':pair(it),'loss_W':P,'Rac_Rdc':2*P/(abs(it)**2*Rdc) if abs(it)>1e-12 else 0,'Rdc_ohm':Rdc,'coil_current_for_1p5mWb_A':100*scale,'loss_at_1p5mWb_W':P*scale**2,'ring_current_at_1p5mWb_A':abs(it)*scale,'skin_depth_mm':1000/math.sqrt(math.pi*MU0*500*3e6*f),'Bcenter_T':pair(mo_getpointvalues(.01,0)[2])}
                row['radial_profile']=[{'r_mm':r,'J_A_per_mm2':pair(mo_getpointvalues(r,0)[7]),'Bz_T':pair(mo_getpointvalues(r,0)[2])} for r in [140.05,140.2,140.5,141,142,144,146,148,150,152,154,156,158,159,159.5,159.8,159.95]]
                finish(row)
    a=[r for r in ROWS if r['id']=='sc-ring-f050-s3-h100'][0];b=[r for r in ROWS if r['id']=='sc-ring-f050-s3-h050'][0]
    check('Замкнутый виток: сгущение',relative(a['metrics']['loss_at_1p5mWb_W'],b['metrics']['loss_at_1p5mWb_W']),.04)
    ring_ports()

def ring_ports():
    base=[r for r in ROWS if r['id']=='sc-ring-f050-s3-h050'][0]
    configurations=[('sc-ring-port-f%03d-%s'%(f,e),f,1 if e=='A' else 0,0 if e=='A' else 1,2,False,False) for f in [10,50,200] for e in ['A','B']]
    configurations += [('sc-ring-flux-selected',50,base['metrics']['coil_current_for_1p5mWb_A'],None,2,False,False),('sc-ring-domain-large',50,100,None,4,False,False)]
    configurations += [('sc-ring-nonlinear-I%04d'%i,50,i,None,2,True,False) for i in [100,500]]
    configurations += [('sc-ring-linear-I0500',50,500,None,2,False,False)]
    for ident,f,coilcurrent,ringcurrent,domain,nonlinear,unused in configurations:
        start(f,True);mi_addmaterial('Стальной замкнутый виток',500,500,0,0,3);mi_addmaterial('Обмотка возбуждения',1,1,0,0,0)
        if nonlinear:
            for b,h in [(0,0),(.4,.4/(MU0*500)),(.8,.8/(MU0*500)),(1.1,2400),(1.3,6000),(1.5,15000),(1.7,50000)]:mi_addbhpoint('Стальной замкнутый виток',b,h)
        mi_addcircprop('Катушка возбуждения',coilcurrent,1)
        if ringcurrent is not None:mi_addcircprop('Порт замкнутого витка',ringcurrent,1)
        bound('Внешняя граница A=0');rectangle(0,-200*domain,350*domain,200*domain,12*domain,'Внешняя граница A=0')
        rectangle(100,-30,110,30,2);label(105,0,'Обмотка возбуждения',2,20,'Катушка возбуждения',200)
        rectangle(140,-10,160,10,.5);label(150,0,'Стальной замкнутый виток',.5,10,'Порт замкнутого витка' if ringcurrent is not None else '<None>',1 if ringcurrent is not None else 0)
        label(230,80,'Воздух',8,1)
        row=solve(ident,'ring-ac-shield',{'frequency_Hz':f,'ring_mu_r':500,'ring_conductivity_MS_per_m':3,'excitation_current_peak_A':coilcurrent,'ring_port_current_peak_A':ringcurrent,'nonlinear_illustrative_curve':nonlinear,'domain_factor':domain,'mesh_ring_mm':.5,'ring_mean_radius_mm':150,'ring_cross_section_mm':[20,20],'turns':200})
        phi=2*math.pi*.13*complex(mo_getpointvalues(130,0)[0]);it=complex(integral(10,7));P=real(integral(10,4))
        row['metrics']={'flux_130mm_Wb':pair(phi),'current_ring_peak_A':pair(it),'loss_W':P,'coil_voltage_V':pair(mo_getcircuitproperties('Катушка возбуждения')[1]),'ring_voltage_V':pair(mo_getcircuitproperties('Порт замкнутого витка')[1]) if ringcurrent is not None else [0,0]}
        finish(row)
    for f in [10,50,200]:
        a=[r for r in ROWS if r['id']=='sc-ring-port-f%03d-A'%f][0];b=[r for r in ROWS if r['id']=='sc-ring-port-f%03d-B'%f][0]
        get=lambda row,key:complex(row['metrics'][key][0],row['metrics'][key][1])
        zba=get(a,'ring_voltage_V');zab=get(b,'coil_voltage_V');zbb=get(b,'ring_voltage_V')
        reference=[r for r in ROWS if r['id']=='sc-ring-f%03d-s3-h050'%f][0]
        measured=get(reference,'current_ring_peak_A');predicted=-100*zba/zbb
        b['impedance']={'Zaa_ohm':pair(get(a,'coil_voltage_V')),'Zab_ohm':pair(zab),'Zba_ohm':pair(zba),'Zbb_ohm':pair(zbb),'predicted_shorted_ring_current_at_100A':pair(predicted),'measured_shorted_ring_current_at_100A':pair(measured)}
        check('Взаимность портовых сопротивлений f=%d'%f,relative(zab,zba),.003)
        check('Короткозамкнутый виток: I2=-Z21 I1/Z22 f=%d'%f,relative(predicted,measured),.003)
    selected=[r for r in ROWS if r['id']=='sc-ring-flux-selected'][0]
    large=[r for r in ROWS if r['id']=='sc-ring-domain-large'][0]
    check('Заданный поток 1.5 мВб: прямой пересчёт',relative(abs(complex(*selected['metrics']['flux_130mm_Wb'])),.0015),.002)
    check('Внешняя область витка: двойные размеры',relative(large['metrics']['loss_W'],base['metrics']['loss_W']),.08)
    save()

def conformal():
    lab='slots-force-partition';r0=5;rr=20;A=.001;scale=10
    for kind in ['physical','mapped']:
        start();bound('Внутренняя поверхность A задано',A);bound('Внешняя поверхность A=0',0)
        if kind=='physical':
            # Дуги задаются точно, без многоугольной аппроксимации.
            for r in [r0,rr]:node(r,0);node(0,r)
            mi_addarc(rr,0,0,rr,90,2);mi_addarc(r0,0,0,r0,90,2)
            mi_selectarcsegment(rr/math.sqrt(2),rr/math.sqrt(2));mi_setarcsegmentprop(2,'Внешняя поверхность A=0',0,1);mi_clearselected()
            mi_selectarcsegment(r0/math.sqrt(2),r0/math.sqrt(2));mi_setarcsegmentprop(2,'Внутренняя поверхность A задано',0,1);mi_clearselected()
            segment((r0,0),(rr,0),.35);segment((0,r0),(0,rr),.35);label(8,8,'Воздух',.35,1)
        else:
            rectangle(0,0,scale*math.log(rr/r0),scale*math.pi/2,.35)
            segment((0,0),(0,scale*math.pi/2),.35,'Внутренняя поверхность A задано')
            segment((scale*math.log(rr/r0),0),(scale*math.log(rr/r0),scale*math.pi/2),.35,'Внешняя поверхность A=0')
            label(5,5,'Воздух',.35,1)
        row=solve('sc-conformal-'+kind,lab,{'mapping':'z=5 mm * exp(w/10 mm)','inner_radius_mm':5,'outer_radius_mm':20,'potential_inner_Wb_per_m':A,'potential_outer_Wb_per_m':0,'coordinate_plane':kind})
        values=[]
        for r,theta in [(6,.2),(7,.4),(8,.6),(10,.8),(12,1),(15,1.2),(18,1.4)]:
            x,y=(r*math.cos(theta),r*math.sin(theta)) if kind=='physical' else (scale*math.log(r/r0),scale*theta)
            p=mo_getpointvalues(x,y);bv=math.sqrt(abs(complex(p[1]))**2+abs(complex(p[2]))**2)
            exactA=A*math.log(rr/r)/math.log(rr/r0);exactB=A/((r if kind=='physical' else scale)*.001*math.log(rr/r0))
            values.append({'physical_r_mm':r,'theta_rad':theta,'A_Wb_per_m':real(p[0]),'analytic_A_Wb_per_m':exactA,'B_T':bv,'analytic_B_T':exactB})
        row['samples']=values;row['metrics']={'max_A_absolute_error_Wb_per_m':max(abs(p['A_Wb_per_m']-p['analytic_A_Wb_per_m']) for p in values),'max_B_relative_error':max(relative(p['B_T'],p['analytic_B_T']) for p in values)}
        check('Конформное отображение '+kind,row['metrics']['max_B_relative_error'],.03)
        finish(row)

if MODE in ['all','gap']:gap_carter()
if MODE in ['all','permeance']:permeance()
if MODE in ['all','force']:forces()
if MODE in ['all','ac']:conductors()
if MODE in ['all','ring']:ring()
if MODE in ['all','conformal']:conformal()
save()
print('SC_COMPLETE',MODE,len(ROWS),all(c['passed'] for c in CHECKS))