"""Потери массивного ротора: дискретизация возбуждения и гармоники."""
exec(open('common.py').read())

def source_amp(rows,order):
    a,b=spatial(rows,order);return max(abs((a+1j*b)/2),abs((a-1j*b)/2))

def case(name,frequency=50,harmonics=None,count=96,mesh=.3,mu=500,nonlinear=False,standing=False):
    if harmonics is None:harmonics=[(2,1)]
    setup(frequency);mi_smartmesh(0)
    if nonlinear:nonlinear_steel();body='Сталь M-43';mi_modifymaterial(body,5,2)
    else:body='Массивная сталь';material(body,mu,2)
    circle(60,1,.3 if mesh<.2 else .6);label(0,0,body,mesh,1)
    label(62,1,'Воздух',.5);label(67,1,'Воздух',1);outer(80)
    # Одинаковая первая пространственная гармоника при изменении числа участков.
    for k in range(count):
        a=360*k/count;b=360*(k+1)/count;c=math.radians((a+b)/2);j=0j
        for order,amp in harmonics:
            sinc=math.sin(order*math.pi/count)/(order*math.pi/count)
            j+=amp*(math.cos(order*c) if standing else complex(math.cos(order*c),math.sin(order*c)))/sinc
        prop='Возбуждающий слой '+str(k+1);material(prop,j=j)
        sector(64,65,a,b,10,.8);label(64.5,(a+b)/2,prop,.35,10)
    r=solve(name,{'family':'massive-rotor','frequency_Hz':frequency,'harmonics':[[h,pair(j)] for h,j in harmonics],'segments':count,'mesh_body_mm':mesh,'mu_r':mu,'nonlinear_background':nonlinear,'conditions':'Массивный цилиндр R60мм, длина100мм, σ2МСм/м; токовый слой64..65мм; внешняя границаR80мм A=0. '+('Статика M-43 для оценки рабочего смещения; проводимость не создает DC потерь.' if frequency==0 else 'Гармоническое поле в системе ротора; амплитудные фазоры; потери только вихретоковые.')},62,body_groups=[1])
    r['Pbody_W']=r['bodies']['1']['P_W'];r['Bsurface_T']=source_amp(r['ring'],harmonics[0][0])
    opendocument(name+'.fem');mi_loadsolution();x,y=polar(40,45);v=mo_getpointvalues(x,y)
    B=math.sqrt(abs(v[1])**2+abs(v[2])**2);H=math.sqrt(abs(v[5])**2+abs(v[6])**2)
    r['probe']={'r_mm':40,'theta_deg':45,'B_T':B,'H_Apm':H,'mu_sec_r':B/(MU0*H) if H else 0}
    close();save_result();return r

cal=case('Ротор-Калибровка');amp=.5/cal['Bsurface_T']
for count in [24,48,96]:case('Ротор-Возбуждение-'+str(count),harmonics=[(2,amp)],count=count)
case('Ротор-Сетка',harmonics=[(2,amp)],mesh=.15)
for f in [10,25,100,300]:case('Ротор-Частота-'+str(f),frequency=f,harmonics=[(2,amp)])
for order in [10,14]:case('Ротор-Гармоника-'+str(order),frequency=300,harmonics=[(order,amp*.2)])
case('Ротор-Гармоники-вместе',frequency=300,harmonics=[(10,amp*.2),(14,amp*.2)])
case('Ротор-Гармоники-сетка',frequency=300,harmonics=[(10,amp*.2),(14,amp*.2)],mesh=.125)
# Локальная секущая проницаемость служит явно ограниченной моделью смещения.
for target in [.3,1.5]:
    current=amp
    for it in range(5):
        bg=case('Ротор-Смещение-'+str(target)+'-'+str(it),frequency=0,harmonics=[(2,current)],nonlinear=True,standing=True)
        err=abs(bg['probe']['B_T']/target-1)
        if err<.01:break
        current*=target/bg['probe']['B_T']
    sec=bg['probe']['mu_sec_r']
    rr=case('Ротор-Смещение-потери-'+str(target),frequency=300,harmonics=[(10,amp*.2),(14,amp*.2)],mu=sec)
    rr['frozen_background_probe_T']=bg['probe']['B_T'];rr['background_model']=bg['model'];rr['frozen_permeability_definition']='локальная секущая B/(μ0H),r40мм,θ45°';save_result()
print('MASSIVE_ROTOR_COMPLETE',len(records))