"""FEMM 5 MicroPython. Собственный учебный магнитостатический RapidIM-опыт.
Токи стержней заданы по первой пространственной гармонике. Новые FE-решения.
"""
import math,os,json,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)))

try:os.mkdir('rapidim-output')
except OSError:pass
def xy(r,a):return r*math.cos(a),r*math.sin(a)
def circle(r):
 for a in [0,math.pi]:
  x,y=xy(r,a);u,v=xy(r,a+math.pi);mi_addnode(x,y);mi_addnode(u,v);mi_addarc(x,y,u,v,180,3)
def sector(r1,r2,a,w):
 a-=w/2;b=a+w
 for r in [r1,r2]:
  x,y=xy(r,a);u,v=xy(r,b);mi_addnode(x,y);mi_addnode(u,v);mi_addarc(x,y,u,v,w*180/math.pi,2)
 for t in [a,b]:
  x,y=xy(r1,t);u,v=xy(r2,t);mi_addsegment(x,y,u,v)
def label(r,a,material,circuit='',group=0,h=2):
 x,y=xy(r,a);mi_addblocklabel(x,y);mi_selectlabel(x,y);mi_setblockprop(material,0,h,circuit,0,group,1 if circuit else 0);mi_clearselected()
def close():
 for n in ['mo_close','mi_close']:
  try:globals()[n]()
  except Exception:pass
def geometry(nonlinear):
 close();newdocument(0);mi_probdef(0,'millimeters','planar',1e-9,80,25)
 for mat,mu in [('Air',1),('Winding',1),('Steel',1000 if not nonlinear else 1)]:mi_addmaterial(mat,mu,mu,0,0,0)
 if nonlinear:
  for b,h in [(0,0),(.25,100),(.5,200),(.9,400),(1.2,800),(1.45,2000),(1.6,5000),(1.7,10000)]:mi_addbhpoint('Steel',b,h)
 for r in [18,58,60,100,150]:circle(r)
 mi_addboundprop('A0',0,0,0,0,0,0,0,0,0)
 for a in [math.pi/2,3*math.pi/2]:
  x,y=xy(150,a);mi_selectarcsegment(x,y);mi_setarcsegmentprop(3,'A0',0,0);mi_clearselected()
 label(0,0,'Air',group=2);label(30,0,'Steel',group=2);label(59,.1,'Air',h=.45);label(95,0,'Steel');label(125,0,'Air',h=8)
 sources=[]
 for family,count,r1,r2,half,coef,group in [('S',36,60.5,82,1.8,40.,0),('R',24,48,57.5,2.5,60.,2)]:
  for k in range(count):
   angle=2*math.pi*(k+.5)/count;name=family+str(k);mi_addcircprop(name,0,1);sector(r1,r2,angle,2*half*math.pi/180);label((r1+r2)/2,angle,'Winding',name,group,1.5)
   sources.append({'name':name,'family':family,'d':coef*math.cos(2*angle),'q':coef*math.sin(2*angle)})
 return sources
def solve(sources,sd,sq,rq,name):
 for c in sources:mi_modifycircprop(c['name'],1,c['d']*sd+c['q']*sq if c['family']=='S' else c['q']*rq)
 path='rapidim-output/'+name+'.fem';mi_saveas(path);mi_analyze(1);mi_loadsolution()
 flux={'sd':0.,'sq':0.,'rd':0.,'rq':0.}
 for c in sources:
  psi=float(mo_getcircuitproperties(c['name'])[2]);tag='s' if c['family']=='S' else 'r';flux[tag+'d']+=2/3*c['d']*psi;flux[tag+'q']+=2/3*c['q']*psi
 mo_clearblock();mo_groupselectblock(2);T=float(mo_blockintegral(22));mo_clearblock()
 flux.update({'irq_A':rq,'isd_A':sd,'isq_A':sq,'T_WST_Nm':T,'T_dq_Nm':3*(flux['sd']*sq-flux['sq']*sd),'elements':mo_numelements(),'model':path})
 mo_close();return flux
results=[]
for nonlinear,sd,sq in [(False,2.,4.),(True,2.,4.),(True,6.,8.)]:
 sources=geometry(nonlinear);tag='nonlinear-' if nonlinear else 'linear-';tag+=str(sd)
 a=solve(sources,sd,sq,-sq,tag+'-0');r1=a['irq_A'];f1=a['rq']
 # Общая секущая: безопасная проба не требует неоднозначных знаков формул оригинала.
 r2=0.;b=solve(sources,sd,sq,r2,tag+'-1');history=[a,b]
 assert f1*b['rq']<0,'Нужен интервал смены знака λrq'
 for it in range(2,25):
  r=r2-b['rq']*(r2-r1)/(b['rq']-f1)
  lo=min(r1,r2);hi=max(r1,r2)
  if min(r-lo,hi-r)<.05*(hi-lo):r=(r1+r2)/2
  c=solve(sources,sd,sq,r,tag+'-'+str(it));history.append(c)
  if abs(c['rq'])<=1e-5*max(abs(c['rd']),1e-8):break
  if c['rq']*f1>0:r1,f1=r,c['rq']
  else:r2,b=r,c
 last=history[-1];last['omega_sl_for_Rr_0_2_rad_s']=-.2*last['irq_A']/last['rd']
 results.append({'nonlinear':nonlinear,'history':history,'passed':abs(last['rq'])<1e-5*abs(last['rd'])})
 if not nonlinear:
  # M и Lr определяем отдельными линейными расчётами; не подгоняем к итерации.
  M=solve(sources,0,1,0,'linear-M')['rq'];Lr=solve(sources,0,0,1,'linear-Lr')['rq'];results[-1]['linear_independent']={'M_H':M,'Lr_H':Lr,'expected_irq_A':-M/Lr*sq,'actual_irq_A':last['irq_A'],'relative_error':abs(last['irq_A']+M/Lr*sq)/abs(last['irq_A'])}
close();out={'scope':'Учебная собственная 36-пазовая статорная и 24-стержневая роторная модель, p=2, глубина80мм, воздушный зазор2мм. Эквивалентные токовые координаты нормированы по виртуальной мощности 3/2; геометрия и B–H не являются паспортом промышленного АД. Равенство эффективных ампервитков обеспечено отношением коэффициентов60/40=36/24. Кольца/скин-эффект не решаются; Rr=0.2Ом — отдельно заданный учебный приведённый параметр.','cases':results,'passed':all(c['passed'] for c in results)}
with open('rapidim-result.json','w') as f:json.dump(out,f)
print('RapidIM',out['passed'],[(c['nonlinear'],len(c['history']),c['history'][-1]['irq_A'],c['history'][-1]['rq']) for c in results])
assert out['passed']