"""Неявное интегрирование в потокосцеплениях с обращением нелинейной карты."""
import math
from numerics import solve_linear

def rhs_and_jacobian(magnetic_map,flux,guess,rs,rf,omega,load_resistance,field_voltage):
    result=magnetic_map.inverse(flux,guess,[1.,1.,1.],1e-10);current=result['x']
    value,jac=magnetic_map.evaluate(current)
    columns=[solve_linear(jac,[1. if j==k else 0. for j in range(3)]) for k in range(3)]
    inverse=[[columns[j][i] for j in range(3)] for i in range(3)]
    resistance=[rs+load_resistance,rs+load_resistance,rf]
    rhs=[-resistance[0]*current[0]+omega*flux[1],-resistance[1]*current[1]-omega*flux[0],field_voltage-rf*current[2]]
    derivative=[[-resistance[i]*inverse[i][j] for j in range(3)] for i in range(3)]
    derivative[0][1]+=omega;derivative[1][0]-=omega
    return rhs,derivative,current

def implicit_flux_step(magnetic_map,old_flux,old_current,step,rs,rf,omega,load_resistance,field_voltage):
    if step<=0:raise ValueError('Неположительный шаг')
    flux=old_flux[:];guess=old_current[:]
    for iteration in range(15):
        rhs,jac,current=rhs_and_jacobian(magnetic_map,flux,guess,rs,rf,omega,load_resistance,field_voltage)
        residual=[flux[j]-old_flux[j]-step*rhs[j] for j in range(3)]
        norm=max(abs(residual[j])/max(abs(flux[j]),1.) for j in range(3))
        if norm<1e-10:return flux,current,iteration+1
        matrix=[[(1. if i==j else 0.)-step*jac[i][j] for j in range(3)] for i in range(3)]
        delta=solve_linear(matrix,[-v for v in residual]);factor=1.;accepted=False
        for attempt in range(20):
            trial=[flux[j]+factor*delta[j] for j in range(3)]
            try:
                rr,_,ii=rhs_and_jacobian(magnetic_map,trial,current,rs,rf,omega,load_resistance,field_voltage)
                error=max(abs(trial[j]-old_flux[j]-step*rr[j])/max(abs(trial[j]),1.) for j in range(3))
                if error<norm:accepted=True;break
            except ValueError:pass
            factor*=.5
        if not accepted:raise ValueError('Неявный шаг вышел из карты или не уменьшает невязку')
        flux=trial;guess=ii
    raise ValueError('Неявный шаг не сошёлся')

def open_field_step(magnetic_map,old_flux_f,old_current_f,step,rf,voltage):
    current=old_current_f
    for iteration in range(20):
        flux,jac=magnetic_map.evaluate([0.,0.,current]);residual=flux[2]-old_flux_f-step*(voltage-rf*current)
        if abs(residual)<1e-9:return flux,[0.,0.,current]
        current-=residual/(jac[2][2]+step*rf)
        if not magnetic_map.axes[2][0]<=current<=magnetic_map.axes[2][-1]:raise ValueError('Возбуждение вышло за карту')
    raise ValueError('Не сошёлся шаг возбуждения')