"""FEMM5 MicroPython. Независимый контролируемый опыт изменения дуги PM.
Геометрия строится заново; оригинальные файлы курса не изменяются.
"""
import math,json,os,sys
assert sys.implementation.name=='micropython'

# FEMM CLI и команда запуска сохранённого файла передают его путь в sys.argv[0].
# Локальные модели/таблицы ищутся возле программы, а не в случайной рабочей папке.
import os
_saved_file=globals().get('__file__') or (sys.argv[0] if sys.argv else '')
if not _saved_file or not os.path.isfile(_saved_file):
 raise RuntimeError('Запустите сохранённый файл из полного распакованного комплекта командой FEMM Run Script; не вставляйте его частями в консоль.')
os.chdir(os.path.dirname(os.path.abspath(_saved_file)))

MU0=4e-7*math.pi
nodes={};edges=set();arcs=set();global_angles=[];pm_label_specs=[]
try:os.mkdir('magnet-arc-output')
except OSError:pass
def node(x,y):
    key=(round(x,7),round(y,7))
    if key not in nodes:mi_addnode(x,y);nodes[key]=1

def line(x1,y1,x2,y2):
    node(x1,y1);node(x2,y2)
    key=tuple(sorted(((round(x1,7),round(y1,7)),(round(x2,7),round(y2,7)))))
    if key not in edges:mi_addsegment(x1,y1,x2,y2);edges.add(key)

def rect(x1,y1,x2,y2):
    line(x1,y1,x2,y1);line(x2,y1,x2,y2);line(x2,y2,x1,y2);line(x1,y2,x1,y1)

def polar(r,a):
    return r*math.cos(math.radians(a)), r*math.sin(math.radians(a))

def arc(r,a,b):
    a=a%360; span=(b-a)%360
    if span<1e-7:return
    cuts=[a]+sorted([a+(t-a)%360 for t in global_angles if 1e-6<(t-a)%360<span-1e-6])+[a+span]
    for left,right in zip(cuts[:-1],cuts[1:]):
        x1,y1=polar(r,left);x2,y2=polar(r,right);node(x1,y1);node(x2,y2)
        key=(round(r,7),round(left%360,7),round(right%360,7))
        if key not in arcs:mi_addarc(x1,y1,x2,y2,right-left,4);arcs.add(key)

def circle(r):
    angles=sorted(set([0,90,180,270]+[round(a%360,7) for a in global_angles]))
    for a,b in zip(angles,angles[1:]+[angles[0]+360]):arc(r,a,b)

def sector(r1,r2,a,b):
    arc(r1,a,b);arc(r2,a,b)
    x1,y1=polar(r1,a);x2,y2=polar(r2,a);line(x1,y1,x2,y2)
    x1,y1=polar(r1,b);x2,y2=polar(r2,b);line(x1,y1,x2,y2)

def label(x,y,material,circuit='',turns=0,mesh=2,group=0,magdir=0):
    mi_addblocklabel(x,y);mi_selectlabel(x,y)
    mi_setblockprop(material,0,mesh,circuit,magdir,group,turns);mi_clearselected()
    if material=='Magnet':
        pm_label_specs.append({'x':x,'y':y,'direction':magdir,
                               'numeric_angle_deg':0 if isinstance(magdir,str) else magdir})

def boundary_segment(x,y):
    mi_selectsegment(x,y);mi_setsegmentprop('A0',0,1,0,0);mi_clearselected()

def setup(frequency=0,kind='planar',depth=40):
    global nodes,edges,arcs,global_angles,pm_label_specs
    nodes={};edges=set();arcs=set();global_angles=[];pm_label_specs=[]
    newdocument(0);mi_probdef(frequency,'millimeters',kind,1e-8,depth,25)
    mi_addmaterial('Air',1,1,0,0,0)
    mi_addmaterial('Winding',1,1,0,0,0)
    mi_addmaterial('LinearSteel',1000,1000,0,0,0)
    mi_addmaterial('NonlinearSteel',1,1,0,0,0)
    for H,B in [(0,0),(100,.25),(200,.5),(400,.9),(800,1.2),(2000,1.45),(5000,1.6),(10000,1.7)]:mi_addbhpoint('NonlinearSteel',B,H)
    mi_addmaterial('Magnet',1.05,1.05,1.2/(MU0*1.05),0,0)
    mi_addboundprop('A0',0,0,0,0,0,0,0,0,0)

def outer_box(size):
    rect(-size,-size,size,size)
    for x,y in [(0,-size),(size,0),(0,size),(-size,0)]:boundary_segment(x,y)
    label(size*.85,size*.85,'Air',mesh=size/20)

def stator(p=2,current=2,angle=0):
    q=12//(6*p);mapping=[('A',1),('C',-1),('B',1),('A',-1),('C',1),('B',-1)]
    for c,a in [('A',0),('B',-120),('C',120)]:mi_addcircprop(c,current*math.cos(math.radians(angle+a)),1)
    circle(60);circle(100);label(0,94,'LinearSteel',mesh=3,group=1)
    for k in range(12):
        center=k*30;sector(60,85,center-7,center+7);x,y=polar(73,center);c,sign=mapping[(k//q)%6];label(x,y,'Winding',c,sign*40,mesh=2,group=1)
def close():
 for n in ['mo_close','mi_close']:
  try:globals()[n]()
  except Exception:pass
def one(half,theta,current,gamma,name):
 close();setup(depth=80)
 global global_angles
 global_angles=[0,90,180,270]
 for k in range(12):global_angles.extend([k*30-7,k*30+7])
 for k in range(4):global_angles.extend([theta+k*90-half,theta+k*90+half])
 outer_box(170);stator(2,current,gamma);circle(18);circle(48);circle(57)
 label(*polar(32,theta+45),'LinearSteel',mesh=2,group=2)
 for k in range(4):
  center=theta+k*90;sector(48,57,center-half,center+half)
  label(*polar(52,center),'Magnet',mesh=1.5,group=2,magdir='theta' if k%2==0 else 'theta+180')
  label(*polar(52,center+45),'Air',mesh=1)
 label(0,0,'Air',mesh=2);label(*polar(59,17),'Air',mesh=.35)
 mi_selectcircle(0,0,58.01,4);mi_setgroup(2);mi_clearselected()
 mi_saveas('magnet-arc-output/'+name+'.fem');mi_analyze(1);mi_loadsolution()
 psi=[float(mo_getcircuitproperties(n)[2]) for n in ['A','B','C']]
 assert all(math.isfinite(v) for v in psi)
 mo_clearblock();mo_groupselectblock(2);T=float(mo_blockintegral(22));mo_clearblock()
 return {'theta_mechanical_deg':theta,'flux_Wb':psi,'T_Nm':T,'elements':mo_numelements()}
N=4 if len(sys.argv)>1 and sys.argv[1]=='quick' else 32
results=[]
for half in [40.,36.]:
 first=one(half,0,0,0,'cal-'+str(half));psi=first['flux_Wb'];alpha=2/3*(psi[0]-.5*psi[1]-.5*psi[2]);beta=math.sqrt(3)/3*(psi[1]-psi[2]);offset=math.atan2(beta,alpha)
 series={};stats={}
 for state in ['open','load']:
  rows=[]
  for k in range(N):
   theta=30*k/N;gamma=2*theta+math.degrees(offset)+90
   rows.append(one(half,theta,0 if state=='open' else 5,gamma,('control-' if N<32 else '')+str(half)+'-'+state+'-'+str(k)))
  series[state]=rows;v=[r['T_Nm'] for r in rows];stats[state]={'mean_T_Nm':sum(v)/len(v),'ripple_pp_Nm':max(v)-min(v),'coarse16_mean_Nm':sum(v[::2])/(N//2),'coarse16_ripple_pp_Nm':max(v[::2])-min(v[::2])}
  if state=='open':
   pd=pq=0.
   for row in rows:
    angle=math.radians(2*row['theta_mechanical_deg'])+offset
    pd+=2/3*sum(row['flux_Wb'][j]*math.cos(angle-j*2*math.pi/3) for j in range(3))
    pq-=2/3*sum(row['flux_Wb'][j]*math.sin(angle-j*2*math.pi/3) for j in range(3))
   amp=math.sqrt(pd*pd+pq*pq)/N;stats[state]['pm_flux_fundamental_Wb']=amp;stats[state]['phase_emf_rms_V_at_1000rpm']=2*2*math.pi*1000/60*amp/math.sqrt(2)
 results.append({'full_magnet_arc_deg':2*half,'magnet_volume_relative':half/40,'series':series,'statistics':stats})
close()
out={'scope':'Собственная SPM Q12/2p4, p2, глубина80мм, PM48..57мм, статор60..100мм, одинаковые B–H/Br/токIамп5А/сеточные настройки/угловая сетка32положения. Изменена только дуга80→72градуса; объём PMуменьшаетсяна10%. Фактические сетки имеют разное число элементов вследствие изменения геометрии, локальные размеры одинаковы. Сталь линейнаяμ1000; потери и нагрев не рассчитываются. ЭДСосновнойгармоникиполучена проекцией потокосцепленийХХ на вращающийсябазис с циклическойсимметрией.','cases':results}
with open('magnet-arc-control-result.json' if N<32 else 'magnet-arc-result.json','w') as f:json.dump(out,f)
print('Magnet arc',[(r['full_magnet_arc_deg'],r['statistics']) for r in results])