"""Зубцовый момент PM-генератора: встроенный Python FEMM 5.

Токи равны нулю. Период 10 механических градусов соответствует 36 пазам
и четырём полюсам. Расчёт выполняется на двух сетках; на мелкой сетке
сравниваются 24 и 48 угловых отсчётов. Момент — интеграл взвешенного
тензора Максвелла по всему ротору, без подшипников и механической нагрузки.
"""
import os,sys,json,math
saved=globals().get('__file__','') or (sys.argv[0] if sys.argv else '')
BASE=globals().get('BASE',os.path.dirname(os.path.abspath(saved)) if saved else os.getcwd())
OUTPUT=globals().get('OUTPUT',BASE)
n=1
while os.path.exists(OUTPUT+'/Расчет-%02d'%n):n+=1
OUT=OUTPUT+'/Расчет-%02d'%n;os.mkdir(OUT)
series=[];checks=[]
def harmonic(rows,order=1):
    n=len(rows)
    c=2*sum(rows[k]['torque_Nm']*math.cos(2*math.pi*order*k/n) for k in range(n))/n
    s=2*sum(rows[k]['torque_Nm']*math.sin(2*math.pi*order*k/n) for k in range(n))/n
    return math.sqrt(c*c+s*s)
for mesh,count in [('Сетка-05',24),('Сетка-025',48)]:
    opendocument(BASE+'/PM-'+mesh+'.fem')
    rows=[];angle=0
    for k in range(count):
        target=10*k/count
        if k:
            mi_selectgroup(20);mi_moverotate(0,0,target-angle);mi_clearselected()
        angle=target
        path=OUT+'/PM-'+mesh+'-%02d'%k+'.fem'
        mi_saveas(path);mi_close();opendocument(path);mi_analyze(1);mi_loadsolution()
        mo_clearblock();mo_groupselectblock(20)
        torque=float(complex(mo_blockintegral(22)).real);mo_clearblock()
        flux=[float(complex(mo_getcircuitproperties('Фаза '+p)[2]).real) for p in 'ABC']
        row={'angle_deg':target,'torque_Nm':torque,'phase_flux_Wb':flux,'elements':mo_numelements(),'model':'PM-'+mesh+'-%02d'%k+'.fem','solution':'PM-'+mesh+'-%02d'%k+'.ans'}
        mo_close();mi_close();opendocument(path);mi_loadsolution()
        mo_clearblock();mo_groupselectblock(20);again=float(complex(mo_blockintegral(22)).real);mo_clearblock()
        row['roundtrip_passed']=abs(torque-again)<1e-10
        checks.append({'id':mesh+'-roundtrip-'+str(k),'passed':row['roundtrip_passed']})
        mo_close()
        rows.append(row)
        if (k+1)%6==0:print('PM_COGGING',mesh,k+1,count)
    mi_close()
    mean=sum(v['torque_Nm'] for v in rows)/count
    item={'mesh':mesh,'samples':count,'rows':rows,'mean_Nm':mean,'peak_to_peak_Nm':max(v['torque_Nm'] for v in rows)-min(v['torque_Nm'] for v in rows),'harmonic_36_Nm':harmonic(rows),'harmonics_36_multiple_Nm':[harmonic(rows,h) for h in range(1,7)]}
    series.append(item)
    with open(OUT+'/Результаты.json','w') as f:json.dump({'series':series,'checks':checks},f)
fine=series[-1];coarse=series[0]
angular=abs(harmonic(fine['rows'][::2])-fine['harmonic_36_Nm'])/max(fine['harmonic_36_Nm'],1e-20)
mesherr=abs(coarse['harmonic_36_Nm']-fine['harmonic_36_Nm'])/max(fine['harmonic_36_Nm'],1e-20)
meanrel=abs(fine['mean_Nm'])/max(fine['peak_to_peak_Nm']/2,1e-20)
checks.extend([{'id':'angular-harmonic36','passed':angular<.02,'value':angular,'limit':.02},{'id':'mesh-harmonic36','passed':mesherr<.10,'value':mesherr,'limit':.10},{'id':'period-mean','passed':meanrel<.05,'value':meanrel,'limit':.05}])
result={'passed':all(c['passed'] for c in checks),'series':series,'checks':checks,'angular_relative':angular,'mesh_relative':mesherr,'mean_to_peak_relative':meanrel,'scope':'2D cogging, zero winding currents; 10deg mechanical period; 36th mechanical torque harmonic; no loaded torque, skew, end effects or mechanical losses'}
with open(OUT+'/Результаты.json','w') as f:json.dump(result,f)
with open(OUT+'/Зубцовый-момент.csv','w') as f:
    f.write('mesh;samples;angle_deg;torque_Nm\n')
    for item in series:
        for row in item['rows']:f.write(';'.join(str(v) for v in [item['mesh'],item['samples'],row['angle_deg'],row['torque_Nm']])+'\n')
print('PM_COGGING_DONE',result['passed'],angular,mesherr,meanrel)