"""Встроенный Python FEMM 5. Нелинейная карта SRM и производная коэнергии.
Рядом должна находиться SRM-основа.fem с русскими материалами и фазами.
Угол задаётся механическими градусами; производная энергии — по радианам.
Построены статические характеристики; питание преобразователем не моделируется.
"""
import os,sys,math,json
script=globals().get('__file__','') or (sys.argv[0] if sys.argv else '')
base=os.path.dirname(os.path.abspath(script)) if script else os.getcwd()
index=1
while 'calculated-%02d'%index in os.listdir(base):index+=1
directory=base+'/calculated-%02d'%index;os.mkdir(directory)
seed=base+'/SRM-основа.fem'
records=[];checks=[]
def scalar(v):return complex(v).real
def check(name,value,limit):checks.append({'name':name,'value':value,'limit':limit,'passed':math.isfinite(value) and value<=limit})
def solve(angle,current):
    opendocument(seed)
    mi_modifycircprop('Фаза А',1,current);mi_modifycircprop('Фаза В',1,0);mi_modifycircprop('Фаза С',1,0)
    mi_selectgroup(2);mi_moverotate(0,0,angle-15);mi_clearselected()
    name='SRM-угол%s-ток%s'%(str(angle).replace('.','p'),str(current).replace('.','p'))
    path=directory+'/'+name+'.fem';mi_saveas(path);mi_close();opendocument(path);mi_analyze(1);mi_loadsolution()
    flux=scalar(mo_getcircuitproperties('Фаза А')[2])
    mo_groupselectblock(2);torque=scalar(mo_blockintegral(22));mo_clearblock();mo_groupselectblock()
    energy=scalar(mo_blockintegral(2));coenergy=scalar(mo_blockintegral(17));mo_clearblock()
    row={'angle_deg':angle,'current_A':current,'flux_Wb':flux,'Lsec_H':flux/current if current else None,'energy_J':energy,'coenergy_J':coenergy,'torque_Nm':torque,'nodes':mo_numnodes(),'elements':mo_numelements(),'model':name+'.fem','solution':name+'.ans'}
    mo_close();mi_close();records.append(row)
    with open(directory+'/results.json','w') as stream:json.dump({'records':records,'checks':checks},stream)
    return row
lookup={}
for angle in [0,15,30,45]:
    for current in [1.5,3,6]:lookup[(angle,current)]=solve(angle,current)
for current in [0,0.5,1,2,2.5,3.5,4,4.5,5,5.5]:lookup[(15,current)]=solve(15,current)
derivatives=[]
for step in [1,0.5]:
    lower=solve(15-step,3);upper=solve(15+step,3)
    derivative=(upper['coenergy_J']-lower['coenergy_J'])/(2*step*math.pi/180)
    derivatives.append({'step_deg':step,'torque_from_coenergy_Nm':derivative,'relative_to_stress':abs(derivative-lookup[(15,3)]['torque_Nm'])/max(abs(derivative),1e-12)})
check('SRM: момент по коэнергии и напряжениям Максвелла',derivatives[-1]['relative_to_stress'],0.06)
check('SRM: уточнение угловой производной',abs(derivatives[-1]['torque_from_coenergy_Nm']-derivatives[0]['torque_from_coenergy_Nm'])/abs(derivatives[-1]['torque_from_coenergy_Nm']),0.03)
current_series=sorted([r for r in records if r['angle_deg']==15],key=lambda r:r['current_A'])
integral=0
for lower,upper in zip(current_series,current_series[1:]):integral+=(upper['current_A']-lower['current_A'])*(upper['flux_Wb']+lower['flux_Wb'])/2
exact=lookup[(15,6)]['coenergy_J']
check('SRM: интеграл Ψ di и коэнергия поля',abs(integral-exact)/abs(exact),0.03)
identities=[abs(r['energy_J']+r['coenergy_J']-r['current_A']*r['flux_Wb'])/max(abs(r['current_A']*r['flux_Wb']),1e-10) for r in records if r['current_A']]
check('SRM: W+Wprime=IΨ',max(identities),0.002)
result={'records':records,'checks':checks,'derivatives':derivatives,'integrated_coenergy_J':integral,'exact_coenergy_J':exact,'current_series':current_series,'passed':all(c['passed'] for c in checks),'scope':'Нелинейная B(H) и одна фаза, статические точки; токи 0..6 А, углы0..45 градусов. Коммутация, механическая динамика, взаимная динамика фаз и потери не рассчитаны.'}
with open(directory+'/results.json','w') as stream:json.dump(result,stream)
print('RELUCTANCE_COMPLETE',result['passed'],len(records),len(checks),directory)