"""Встроенный Python FEMM5. Три учебные геометрии на основе Katona, независимый адаптер.
Новые серии: 32 положения ХХ и 32 положения нагрузки на зубцовый период.
Все токи пиковые. Nф=40 последовательных витков, Iq=5А, d=40мм.
"""
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)))

try:os.mkdir('katona-output')
except OSError:pass
def close():
 for n in ['mo_close','mi_close']:
  try:globals()[n]()
  except Exception:pass
def labels(path):
 with open(path,'r') as f:lines=f.read().split('\n')
 names=[];rows=[];inlabels=False;count=0
 for line in lines:
  if '<BlockName>' in line:names.append(line.split('"')[1])
  if '[NumBlockLabels]' in line:inlabels=True;count=int(line.split('=')[1]);continue
  if inlabels and len(rows)<count:rows.append(line.split())
 return names,rows
def load(case):
 path='katona-'+str(case)+'-course.fem';close();opendocument(path);mi_probdef(0,'millimeters','planar',1e-8,40,25)
 materials,rows=labels(path);counts=[0,0,0]
 for row in rows:
  name=materials[int(row[2])-1]
  if name in ['A','B','C','a','b','c']:counts['ABC'.index(name.upper())]+=1
 assert len(set(counts))==1
 # Каждая сторона катушки содержит N, число последовательных витков фазы = сумма |N|/2.
 turns=80/counts[0]
 for name in ['A','B','C']:mi_addcircprop('Фаза '+name,0,1)
 mi_addmaterial('Идеальная обмотка',1,1,0,0,0)
 for row in rows:
  mat=materials[int(row[2])-1]
  if mat in ['A','B','C','a','b','c']:
   x,y=float(row[0]),float(row[1]);mi_selectlabel(x,y);mi_setblockprop('Идеальная обмотка',0,1,'Фаза '+mat.upper(),0,0,turns if mat.isupper() else -turns);mi_clearselected()
 # Группа ротора заранее назначена по радиусам обоих концов каждой линии/дуги.
 # Выбор окружностью опасен для длинной внешней дуги, чья хорда проходит через ротор.
 # Воздушная метка зазора не входит в интегрируемое тело ротора.
 for row in rows:
  mat=materials[int(row[2])-1]
  if mat=='Воздух':
   x,y=float(row[0]),float(row[1]);mi_selectlabel(x,y);mi_setblockprop(mat,1,0,'',0,0,0);mi_clearselected()
 return {'coil_side_count_per_phase':counts,'turns_per_side':turns,'series_turns_per_phase':40}
def sample(case,name,theta,beta,offset):
 load(case);p=2 if case<3 else 4
 mi_selectgroup(2);mi_moverotate(0,0,theta);mi_clearselected()
 angle=math.radians(p*theta)+offset
 currents=[0.,0.,0.] if name=='open' else [-5*math.sin(angle-k*2*math.pi/3) for k in range(3)]
 for k,n in enumerate(['A','B','C']):mi_modifycircprop('Фаза '+n,1,currents[k])
 path='katona-output/case-'+str(case)+'-'+name+'-'+str(round(theta,4))+'.fem';mi_saveas(path);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();row={'mechanical_angle_deg':theta,'currents_A':currents,'flux_Wb':psi,'T_Nm':T,'elements':mo_numelements()};close();return row
N=2 if len(sys.argv)>1 and sys.argv[1]=='quick' else 32
out={'source':'Katona three documented FEM geometries; independent FEMM5 embedded-Python adapter, not the original external py2femm program.','scope':'Сравнение топологий с одинаковыми Nф=40, Iq=5А и глубиной40мм. Пазовые области идеализированы как заданные токовые обмотки с нулевой проводимостью; потери меди не решаются. Геометрии исходника, число полюсов, площадь меди и длина лобовых частей различаются: это не изолированный опыт укорочения шага и не сравнение КПД. Расчётные числа не являются числами публикации Katona.','angles_per_period':N,'cases':[]}
for case in [1,2,3]:
 meta=load(case);close();zero=sample(case,'open',0,0,0)
 alpha=2/3*(zero['flux_Wb'][0]-.5*zero['flux_Wb'][1]-.5*zero['flux_Wb'][2]);beta=2/3*math.sqrt(3)/2*(zero['flux_Wb'][1]-zero['flux_Wb'][2]);offset=math.atan2(beta,alpha)
 period=30 if case<3 else 15;series={'open':[],'load':[]}
 for state in ['open','load']:
  for k in range(N):series[state].append(sample(case,state,k*period/N,0,offset))
 stats={}
 for state,rows in series.items():
  v=[r['T_Nm'] for r in rows];coarse=v[::2];stats[state]={'mean_T_Nm':sum(v)/len(v),'ripple_pp_Nm':max(v)-min(v),'coarse_16_mean_Nm':sum(coarse)/len(coarse),'coarse_16_ripple_pp_Nm':max(coarse)-min(coarse)}
  if state=='open':
   vd=vq=0
   for row in rows:
    angle=math.radians((2 if case<3 else 4)*row['mechanical_angle_deg'])+offset
    vd+=2/3*sum(row['flux_Wb'][k]*math.cos(angle-k*2*math.pi/3) for k in range(3))
    vq-=2/3*sum(row['flux_Wb'][k]*math.sin(angle-k*2*math.pi/3) for k in range(3))
   stats[state]['pm_linkage_fundamental_projection_Wb']=math.sqrt(vd*vd+vq*vq)/len(rows)
 out['cases'].append({'case':case,'p':2 if case<3 else 4,'mechanical_period_deg':period,'electrical_offset_rad':offset,'winding':meta,'series':series,'statistics':stats})
with open('katona-control-result.json' if N<32 else 'katona-result.json','w') as f:json.dump(out,f)
print('Katona comparison',[(c['case'],c['statistics']) for c in out['cases']])