"""Согласованные единицы и нормировка; встроенный Python FEMM 5.

Токи положительны внутрь машины. q-строка Парка = -sin. Напряжение и
потокосцепление имеют ту же полярность, что и ток. Положительная отдаваемая
мощность равна минус электрической мощности на зажимах.
"""
import math


def abc_to_dq0(abc, theta):
    angles = [theta, theta-2*math.pi/3, theta+2*math.pi/3]
    return [2/3*sum(abc[k]*math.cos(angles[k]) for k in range(3)),
            -2/3*sum(abc[k]*math.sin(angles[k]) for k in range(3)), sum(abc)/3]


def dq0_to_abc(dq0, theta):
    d,q,zero=dq0
    return [d*math.cos(a)-q*math.sin(a)+zero
            for a in (theta,theta-2*math.pi/3,theta+2*math.pi/3)]


def power_dq0(voltage, current):
    return 1.5*(voltage[0]*current[0]+voltage[1]*current[1])+3*voltage[2]*current[2]


def torque_dq(flux, current, pole_pairs=1):
    """Н·м, полный активный участок. Только усреднённая dq-модель."""
    return 1.5*pole_pairs*(flux[0]*current[1]-flux[1]*current[0])


def flux_derivative(flux, currents, voltages, rs, rf, omega_e):
    """dλ/dt, В; порядок d,q,f. Скорость задана, вихревые контуры исключены."""
    return [voltages[0]-rs*currents[0]+omega_e*flux[1],
            voltages[1]-rs*currents[1]-omega_e*flux[0],
            voltages[2]-rf*currents[2]]


def harmonics(samples, maximum):
    """Равномерный полный период без повторной конечной точки; пиковые a,b."""
    n=len(samples)
    if n<3 or maximum<0 or maximum>=n/2:
        raise ValueError('Недостаточно отсчётов: порядок должен быть меньше N/2')
    rows=[{'order':0,'cos':sum(samples)/n,'sin':0.0,'peak':abs(sum(samples)/n)}]
    for h in range(1,maximum+1):
        a=2/n*sum(samples[k]*math.cos(2*math.pi*h*k/n) for k in range(n))
        b=2/n*sum(samples[k]*math.sin(2*math.pi*h*k/n) for k in range(n))
        rows.append({'order':h,'cos':a,'sin':b,'peak':math.sqrt(a*a+b*b)})
    return rows


def induced_emf(harmonic, omega_e):
    """e=-dλ/dt, при λ=a cos(hθ)+b sin(hθ); RMS отдельно от пика."""
    h=harmonic['order']; a=-h*omega_e*harmonic['sin']; b=h*omega_e*harmonic['cos']
    return {'order':h,'cos':a,'sin':b,'rms':math.sqrt(a*a+b*b)/math.sqrt(2)}


def relative_change(fine, coarse, absolute_scale):
    if absolute_scale<=0:raise ValueError('Нужна положительная абсолютная шкала')
    return abs(fine-coarse)/max(abs(fine),absolute_scale)


def maxwell_circle(samples, radius_m, length_m, axis_m=(0.0,0.0)):
    """Контур в воздухе вокруг всей подвижной части; Bx,By в серединах интервалов.

    samples = [(angle_rad,Bx_T,By_T), ...], равномерный полный круг.
    Возвращает Fx,Fy и момент относительно заданной оси, все в СИ.
    """
    n=len(samples)
    if n<8:raise ValueError('Нужен полный контур с достаточным числом отсчётов')
    for k in range(1,n):
        if abs((samples[k][0]-samples[k-1][0])-2*math.pi/n)>1e-8:
            raise ValueError('Контур должен быть равномерным полным кругом')
    mu0=4e-7*math.pi; fx=fy=torque=0.0; ds=radius_m*length_m*2*math.pi/n
    for angle,bx,by in samples:
        nx,ny=math.cos(angle),math.sin(angle); bn=bx*nx+by*ny; b2=bx*bx+by*by
        dx=(bn*bx-b2*nx/2)/mu0*ds; dy=(bn*by-b2*ny/2)/mu0*ds
        fx+=dx; fy+=dy
        torque+=(radius_m*nx-axis_m[0])*dy-(radius_m*ny-axis_m[1])*dx
    return [fx,fy,torque]


def winding_factor(winding, slot_count, span, harmonic=1):
    """Комплексная сумма всех сторон одной фазы, нормированная их числом.

    Геометрический угол — электрический для p=1. Здесь проверяется таблица,
    а не предполагается синусоидальная форма распределения.
    """
    result={}
    for phase in 'ABC':
        re=im=0.0; count=0
        for row in winding:
            angle=2*math.pi*(row['slot']-1)/slot_count
            for layer in ('lower','upper'):
                if row[layer+'_phase']==phase:
                    sign=row[layer+'_sign']; re+=sign*math.cos(harmonic*angle)
                    im+=sign*math.sin(harmonic*angle); count+=1
        result[phase]={'factor':math.sqrt(re*re+im*im)/count,'sides':count,'cos':re/count,'sin':im/count}
    return result