"""Встроенный Python FEMM 5: подключение нагрузки к полевой d/q-модели.
Параметры Ld,Lq,Ψf читаются из field-parameters.json рядом с программой.
Это редуцированная электрическая динамика при фиксированной скорости;
пространственное поле в каждом временном шаге повторно не решается.
"""
import os,sys,math,json
script=globals().get('__file__','') or (sys.argv[0] if sys.argv else '')
base=os.path.dirname(os.path.abspath(script)) if script else os.getcwd()
with open(base+'/field-parameters.json') as stream:parameters=json.load(stream)
Ld=parameters['Ld_H'];Lq=parameters['Lq_H'];psif=parameters['psi_f_Wbturn'];Rs=parameters['Rphase_Ohm'];p=parameters['pole_pairs'];we=2*math.pi*parameters['frequency_Hz'];omega=we/p
index=1
while 'calculated-%02d'%index in os.listdir(base):index+=1
directory=base+'/calculated-%02d'%index;os.mkdir(directory)
records=[];checks=[]
def check(name,value,limit):checks.append({'name':name,'value':value,'limit':limit,'passed':math.isfinite(value) and value<=limit})
def derivative(state,load):
    id,iq=state[0],state[1]
    did=(-(Rs+load)*id+we*Lq*iq)/Ld
    diq=(-(Rs+load)*iq-we*(Ld*id+psif))/Lq
    torque=1.5*p*(psif*iq+(Ld-Lq)*id*iq)
    drive=-torque*omega;pl=1.5*load*(id*id+iq*iq);pcu=1.5*Rs*(id*id+iq*iq)
    return [did,diq,drive,pl,pcu]
def transient(load,dt):
    state=[0.0]*5;samples=[];maximum_balance=0;power_identity=0
    for step in range(int(round(0.01/dt))):
        a=derivative(state,load);b=derivative([v+dt*dv/2 for v,dv in zip(state,a)],load)
        c=derivative([v+dt*dv/2 for v,dv in zip(state,b)],load);d=derivative([v+dt*dv for v,dv in zip(state,c)],load)
        state=[v+dt*(v1+2*v2+2*v3+v4)/6 for v,v1,v2,v3,v4 in zip(state,a,b,c,d)]
        id,iq=state[0],state[1];time=(step+1)*dt
        energy=.75*(Ld*id*id+Lq*iq*iq)
        balance=state[2]-state[3]-state[4]-energy;maximum_balance=max(maximum_balance,abs(balance))
        abc=[id*math.cos(we*time-2*math.pi*k/3)-iq*math.sin(we*time-2*math.pi*k/3) for k in range(3)]
        pabc=load*sum(v*v for v in abc);pdq=1.5*load*(id*id+iq*iq)
        power_identity=max(power_identity,abs(pabc-pdq))
        if (step+1)%max(1,int(round(.0001/dt)))==0:
            samples.append({'time_s':time,'id_peak_A':id,'iq_peak_A':iq,'storage_J':energy,'drive_energy_J':state[2],'load_energy_J':state[3],'copper_energy_J':state[4],'balance_J':balance})
    iq_exact=-we*psif/((Rs+load)+we*we*Ld*Lq/(Rs+load));id_exact=we*Lq*iq_exact/(Rs+load)
    end=derivative(state,load)
    return {'Rload_Ohm':load,'dt_s':dt,'samples':samples,'end_id_A':state[0],'end_iq_A':state[1],'steady_id_A':id_exact,'steady_iq_A':iq_exact,'steady_relative_error':math.sqrt((state[0]-id_exact)**2+(state[1]-iq_exact)**2)/math.sqrt(id_exact**2+iq_exact**2),'max_energy_balance_J':maximum_balance,'relative_energy_balance':maximum_balance/max(abs(state[2]),1e-12),'abc_dq_power_error_W':power_identity,'end_Pdrive_W':end[2],'end_Pload_W':end[3],'end_Pcu_W':end[4]}
for load in [50,25]:
    series=[transient(load,dt) for dt in [20e-6,10e-6,5e-6]]
    records.extend(series)
    for row in series:
        check('Установившийся ток R=%g dt=%g'%(load,row['dt_s']),row['steady_relative_error'],1e-7)
        check('Баланс энергии R=%g dt=%g'%(load,row['dt_s']),row['relative_energy_balance'],1e-5)
        check('Мощность abc/dq R=%g'%load,row['abc_dq_power_error_W'],1e-8)
    fine=series[-1]['samples'];coarse=series[-2]['samples']
    difference=max(math.sqrt((a['id_peak_A']-b['id_peak_A'])**2+(a['iq_peak_A']-b['iq_peak_A'])**2) for a,b in zip(fine,coarse))
    check('Уточнение шага тока R=%g'%load,difference,1e-5)
result={'parameters':parameters,'records':records,'checks':checks,'passed':all(r['passed'] for r in checks),'scope':'Фиксированная скорость; постоянное возбуждение; линейные Ld/Lq из трёх независимых полевых проб. Rфазы задано2Ом, активная нагрузка25/50Ом. Не transientFE и не динамика привода/насыщения.'}
with open(directory+'/results.json','w') as stream:json.dump(result,stream)
print('DQ_TRANSIENT_COMPLETE',result['passed'],len(records),len(checks),directory)