"""Полюсные коэффициенты и генератор: встроенный Python FEMM 5."""
exec(open('common.py').read())

def rotor_poles():
    for k in range(4):
        c=90*k
        line(35,c-20,54,c-20,2);arc(54,c-32.5,c-20,2)
        line(54,c-32.5,60,c-32.5,2);arc(60,c-32.5,c+32.5,2,.6)
        line(60,c+32.5,54,c+32.5,2);arc(54,c+20,c+32.5,2)
        line(54,c+20,35,c+20,2);arc(35,c+20,c+70,2)

def stator(slotted=False):
    if not slotted:circle(63,1,.6)
    else:
        for k in range(36):
            c=5+10*k
            arc(63,c-8.8,c-1.2,1,.6)
            line(63,c-1.2,66,c-1.2,1);line(66,c-1.2,69,c-2.5,1)
            line(69,c-2.5,78,c-2.5,1);arc(78,c-2.5,c+2.5,1)
            line(78,c+2.5,69,c+2.5,1);line(69,c+2.5,66,c+1.2,1);line(66,c+1.2,63,c+1.2,1)
    circle(90,1);label(83,0,'Сталь M-43' if slotted else 'Сталь линейная',1,1)

def build(field=0,d=0,q=0,smooth=False,count=96,mesh=.4,mu=10000,slotted=False,N=100):
    setup();mi_smartmesh(0)
    steel(mu=mu)
    if slotted:nonlinear_steel()
    stator(slotted);outer(100);label(61,0,'Воздух',mesh)
    if smooth:circle(60,2,.6)
    else:rotor_poles()
    label(0,0,'Сталь линейная',1.2,2)
    if not smooth:
        for k in range(4):
            c=90*k;sign=1 if k%2==0 else -1
            for side,a,b,sgn in [('Левая',-40,-25,-1),('Правая',25,40,1)]:
                name='Возбуждение '+str(k+1)+' '+side;material(name,j=field*sign*sgn)
                sector(36,52,c+a,c+b,20,.8);label(44,c+(a+b)/2,name,.65,20)
    # Идеальная распределенная трехфазная обмотка: 100 витков/фазу.
    label(62.95,1,'Воздух',.08)
    area=.5*((62.8e-3)**2-(62.2e-3)**2)*2*math.pi/count
    weights=[]
    for k in range(count):
        a=360*k/count;b=360*(k+1)/count;c=(a+b)/2;t=math.radians(2*c)
        n=[N*.5*math.sin(t-ph)*2*math.pi/count for ph in [0,2*math.pi/3,-2*math.pi/3]]
        currents=[d,-.5*d+math.sqrt(3)*.5*q,-.5*d-math.sqrt(3)*.5*q]
        J=sum(n[i]*currents[i] for i in range(3))/area/1e6
        name='Обмотка участок '+str(k+1);material(name,j=J)
        sector(62.2,62.8,a,b,100+k,.8);label(62.5,c,name,.18,100+k);weights.append(n)
    return weights

def point_H(r,a):
    x,y=polar(r,a);v=mo_getpointvalues(x,y);return scalar(v[5]),scalar(v[6])

def radial_H(r1,r2,a,n=80):
    t=math.radians(a);dr=(r2-r1)*.001/n;v=0
    for k in range(n):
        h=point_H(r1+(k+.5)*(r2-r1)/n,a);v+=(h[0]*math.cos(t)+h[1]*math.sin(t))*dr
    return v

def yoke_H():
    v=radial_H(78,84,0)-radial_H(78,84,90)
    for k in range(180):
        a=(k+.5)*.5;t=math.radians(a);h=point_H(84,a)
        v+=(-h[0]*math.sin(t)+h[1]*math.cos(t))*.084*math.pi/360
    return v

def case(name,field=0,d=0,q=0,smooth=False,count=96,mesh=.4,mu=10000,slotted=False,flux=False):
    weights=build(field,d,q,smooth,count,mesh,mu,slotted)
    r=solve(name,{'family':'pole-generator','field_J_MApm2':field,'id_peak_A':d,'iq_peak_A':q,'smooth_rotor':smooth,'segments':count,'mu_reference':mu,'slotted_nonlinear_stator':slotted,'conditions':'Четырехполюсная машина: Rр60мм, Rс63мм, длина100мм; заданные токи, статика; идеальная синусоидальная обмотка100витков/фазу. Линейный ротор, '+('пазовый нелинейный статор M-43.' if slotted else 'гладкий статор с μr='+str(mu)+'.')},61.5,[2,20],list(range(100,100+count)) if flux else None)
    r['B1_T']=real_spatial_amp(r['ring'],2)
    if flux:
        psi=[0,0,0]
        for k in range(count):
            body=r['bodies'][str(100+k)];avg=scalar(complex(*body['A_integral']))/body['area_m2']
            for i in range(3):psi[i]+=avg*weights[k][i]
        r['flux_abc_Wbturn']=psi;r['psi_d_Wbturn']=(2*psi[0]-psi[1]-psi[2])/3;r['psi_q_Wbturn']=(psi[1]-psi[2])/math.sqrt(3)
    if slotted:
        opendocument(name+'.fem');mi_loadsolution();mo_smooth('off')
        fg=radial_H(60.01,62.99,0)-radial_H(60.01,62.99,90)
        ft=radial_H(63.01,77.99,0)-radial_H(63.01,77.99,90);fy=yoke_H()
        r['mmf_At']={'gap':fg,'teeth':ft,'yoke':fy,'kz':(fg+ft+fy)/fg};close()
    save_result();return r


cal=case('Генератор-Калибровка',field=1)
jf=.85/abs(complex(*cal['ring'][0][0]).real)
main=case('Генератор-Возбуждение',field=jf,flux=True)
rd=case('Генератор-Ось-d',d=10,flux=True)
rq=case('Генератор-Ось-q',q=10,flux=True)
pf=main['psi_d_Wbturn'];ld=rd['psi_d_Wbturn']/10;lq=rq['psi_q_Wbturn']/10
omega=2*math.pi*50;rs=2
for rl in [200,100,50,25]:
    rr=rs+rl;iq=-omega*pf/(rr+omega*omega*ld*lq/rr);id=omega*lq*iq/rr
    r=case('Генератор-Нагрузка-'+str(rl),field=jf,d=id,q=iq,flux=True)
    r.update({'Rphase_ohm':rs,'Rload_ohm':rl,'Ld_H':ld,'Lq_H':lq,'psi_f_Wbturn':pf,'frequency_Hz':50,'Pload_W':1.5*rl*(id*id+iq*iq),'Pcu_W':1.5*rs*(id*id+iq*iq),'Uphase_RMS_V':rl*math.sqrt((id*id+iq*iq)/2),'Iphase_RMS_A':math.sqrt((id*id+iq*iq)/2),'T_dq_Nm':3*(r['psi_d_Wbturn']*iq-r['psi_q_Wbturn']*id)})
    save_result()
    if rl==25:
        fine=case('Генератор-Нагрузка-сетка',field=jf,d=id,q=iq,flux=True,count=192,mesh=.2)
        fine.update({'Rload_ohm':rl,'T_dq_Nm':3*(fine['psi_d_Wbturn']*iq-fine['psi_q_Wbturn']*id)});save_result()
print('GENERATOR_COMPLETE',len(records))