"""Нагрузка PM-генератора: встроенный Python FEMM 5.

Проверяется основная гармоника при синусоидальном токе, постоянной
скорости1500об/мин и заданной активной нагрузке на фазу. Высшие гармоники
напряжения не замыкаются на нагрузку. Работа магнитов оценивается по
обратимой линии μrec=1,05,Br=μ0μrecHc; необратимое размагничивание требует
отдельной температурной характеристики изготовителя.
"""
import os,sys,math,json
saved=globals().get('__file__','') or (sys.argv[0] if sys.argv else '')
BASE=globals().get('BASE',os.path.dirname(os.path.abspath(saved)) if saved else os.getcwd())
OUTPUT=globals().get('OUTPUT',BASE)
n=1
while os.path.exists(OUTPUT+'/Расчет-%02d'%n):n+=1
OUT=globals().get('REUSE_OUTPUT','') or OUTPUT+'/Расчет-%02d'%n
if not os.path.exists(OUT):os.mkdir(OUT)
reference=json.load(open(BASE+'/Опорный-ХХ.json'))
PHI=reference['flux_phase_rad'];PSIPM=reference['psi_peak_Wb']
MU0=4*math.pi*1e-7;MUREC=1.05;BR=MU0*MUREC*900000
OMEGA=2*math.pi*50;OMEGAM=2*math.pi*1500/60;RPHASE=1.0;COUNT=24
NAMES=['Фаза A','Фаза B','Фаза C'];RUNS=[];CHECKS=[]
def real(x):return float(complex(x).real)
def basis(theta):
    return [math.cos(theta+PHI+2*math.pi*k/3) for k in range(3)],[-math.sin(theta+PHI+2*math.pi*k/3) for k in range(3)]
def field(tag,angle,id_peak,iq_peak):
    path=OUT+'/'+tag+'.fem'
    cached=open(path,'rb').read() if os.path.exists(path) and os.path.exists(OUT+'/'+tag+'.ans') else None
    opendocument(BASE+'/PM-Нагрузка.fem')
    if angle:mi_selectgroup(20);mi_moverotate(0,0,angle);mi_clearselected()
    d,q=basis(2*math.radians(angle));currents=[id_peak*d[k]+iq_peak*q[k] for k in range(3)]
    for name,value in zip(NAMES,currents):mi_modifycircprop(name,1,value)
    mi_saveas(path);mi_close();opendocument(path)
    if cached is None or cached!=open(path,'rb').read():mi_analyze(1)
    mi_loadsolution()
    flux=[real(mo_getcircuitproperties(name)[2]) for name in NAMES]
    pd=2*sum(flux[k]*d[k] for k in range(3))/3;pq=2*sum(flux[k]*q[k] for k in range(3))/3
    mo_clearblock();mo_groupselectblock(20);torque=real(mo_blockintegral(22));mo_clearblock()
    bpar=[];hpar=[];samples=[];minimum=None;mu_deviation=0
    for pole in range(4):
        sign=1 if pole%2==0 else -1
        for r in [45,45.5,48.525,51.5,52.05]:
            for offset in [-33,-30,-25,-15,0,15,25,30,33]:
                theta=math.radians(pole*90+angle+offset)
                values=mo_getpointvalues(r*math.cos(theta),r*math.sin(theta))
                bp=sign*(real(values[1])*math.cos(theta)+real(values[2])*math.sin(theta));hr=(bp-BR)/(MU0*MUREC)
                bpar.append(bp);hpar.append(hr);samples.append([pole,r,offset,bp,hr])
                mu_deviation=max(mu_deviation,abs(real(values[9])-MUREC),abs(real(values[10])-MUREC))
                if minimum is None or hr<minimum['Hrecoil_Apm']:minimum={'pole':pole,'radius_mm':r,'offset_deg':offset,'Bparallel_T':bp,'Hrecoil_Apm':hr}
    row={'tag':tag,'angle_deg':angle,'id_peak_A':id_peak,'iq_peak_A':iq_peak,'phase_currents_peak_A':currents,'flux_Wb':flux,'psi_d_Wb':pd,'psi_q_Wb':pq,'torque_Nm':torque,'Bparallel_min_T':min(bpar),'Bparallel_max_T':max(bpar),'Hparallel_min_Apm':min(hpar),'Hparallel_max_Apm':max(hpar),'magnet_samples':samples,'sample_columns':['pole','radius_mm','offset_deg','Bparallel_T','Hrecoil_Apm'],'minimum_recoil_position':minimum,'sample_count':180,'elements':mo_numelements(),'model':tag+'.fem','solution':tag+'.ans'}
    mo_close();mi_close();opendocument(path);mi_loadsolution()
    again=[real(mo_getcircuitproperties(name)[2]) for name in NAMES]
    row['roundtrip_passed']=max(abs(a-b) for a,b in zip(again,flux))<1e-10
    mo_close();mi_close();RUNS.append(row)
    CHECKS.append({'id':tag+'-roundtrip','passed':row['roundtrip_passed']})
    CHECKS.append({'id':tag+'-magnet-region','passed':mu_deviation<1e-9,'value':mu_deviation,'limit':1e-9})
    return row

zero=field('Опорное-поле',0,0,0)
dp=field('Приращение-d-плюс',0,1,0);dm=field('Приращение-d-минус',0,-1,0)
qp=field('Приращение-q-плюс',0,0,1);qm=field('Приращение-q-минус',0,0,-1)
LD=(dp['psi_d_Wb']-dm['psi_d_Wb'])/2;LQ=(qp['psi_q_Wb']-qm['psi_q_Wb'])/2
CHECKS.append({'id':'positive-dq-L','passed':LD>0 and LQ>0,'value':[LD,LQ]})
loads=[]
for resistance in [400,200,100]:
    a=resistance+RPHASE
    iq=-OMEGA*PSIPM/(a+OMEGA*OMEGA*LD*LQ/a);id=OMEGA*LQ*iq/a
    history=[]
    for iteration in range(3):
        rows=[]
        for k in range(COUNT):
            row=field('R%d-итерация%d-%02d'%(resistance,iteration,k),180*k/COUNT,id,iq);rows.append(row)
        pd=sum(r['psi_d_Wb'] for r in rows)/COUNT;pq=sum(r['psi_q_Wb'] for r in rows)/COUNT
        vd=RPHASE*id-OMEGA*pq;vq=RPHASE*iq+OMEGA*pd
        fd=vd+resistance*id;fq=vq+resistance*iq
        scale=OMEGA*PSIPM;residual=math.sqrt(fd*fd+fq*fq)/scale
        history.append({'iteration':iteration,'id_peak_A':id,'iq_peak_A':iq,'psi_d_Wb':pd,'psi_q_Wb':pq,'residual_relative':residual})
        print('PM_LOAD',resistance,iteration,residual)
        if residual<.005:break
        if iteration<2:
            det=a*a+OMEGA*OMEGA*LD*LQ
            id+=(-a*fd-OMEGA*LQ*fq)/det
            iq+=(OMEGA*LD*fd-a*fq)/det
    irms=math.sqrt(id*id+iq*iq)/math.sqrt(2);urms=math.sqrt(vd*vd+vq*vq)/math.sqrt(2)
    pout=1.5*resistance*(id*id+iq*iq);pcu=1.5*RPHASE*(id*id+iq*iq)
    torque_raw=sum(r['torque_Nm'] for r in rows)/COUNT
    baseline=sum(reference['zero_current_torque_Nm'])/COUNT
    # Separate the current-produced torque from the paired zero-current field.
    torque=torque_raw-baseline
    pin=-torque*OMEGAM;balance=(pin-pout-pcu)/max(pout+pcu,1)
    # The same balanced current series gives a three-phase fundamental power.
    pem=1.5*OMEGA*(pd*iq-pq*id)
    field_balance=(-pem-pout-pcu)/max(pout+pcu,1)
    row={'Rload_ohm':resistance,'id_peak_A':id,'iq_peak_A':iq,'Iphase_rms_A':irms,'Uphase_rms_V':urms,'Pload_W':pout,'Pcu_W':pcu,'Tmean_Nm':torque,'Tmean_raw_Nm':torque_raw,'zero_current_torque_sample_mean_Nm':baseline,'Pmechanical_input_W':pin,'mechanical_power_residual_relative':balance,'field_terminal_power_residual_relative':field_balance,'voltage_residual_relative':residual,'Bparallel_min_T':min(r['Bparallel_min_T'] for r in rows),'Hparallel_min_Apm':min(r['Hparallel_min_Apm'] for r in rows),'history':history,'rows':rows}
    loads.append(row)
    CHECKS.extend([{'id':'load-%d-terminal'%resistance,'passed':residual<.005,'value':residual,'limit':.005},{'id':'load-%d-power'%resistance,'passed':abs(balance)<.03,'value':balance,'limit':.03},{'id':'load-%d-field-power'%resistance,'passed':abs(field_balance)<.01,'value':field_balance,'limit':.01}])
    with open(OUT+'/Результаты.json','w') as f:json.dump({'loads':loads,'runs':RUNS,'checks':CHECKS},f)

last=loads[-1];fine_rows=[]
for k in range(48):fine_rows.append(field('R100-угол-48-%02d'%k,180*k/48,last['id_peak_A'],last['iq_peak_A']))
pd=sum(r['psi_d_Wb'] for r in fine_rows)/48;pq=sum(r['psi_q_Wb'] for r in fine_rows)/48
fine_torque=sum(r['torque_Nm'] for r in fine_rows)/48-sum(reference['zero_current_torque_48_Nm'])/48
torque_error=abs(fine_torque-last['Tmean_Nm'])/max(abs(fine_torque),1e-20)
fine_voltage=math.sqrt(((RPHASE+100)*last['id_peak_A']-OMEGA*pq)**2+((RPHASE+100)*last['iq_peak_A']+OMEGA*pd)**2)/(OMEGA*PSIPM)
refinement={'samples_coarse':24,'samples_fine':48,'Rload_ohm':100,'Tload_24_Nm':last['Tmean_Nm'],'Tload_48_Nm':fine_torque,'Tload_relative_change':torque_error,'terminal_residual_48_relative':fine_voltage,'rows':fine_rows}
CHECKS.extend([{'id':'loaded-torque-angle-refinement','passed':torque_error<.03,'value':torque_error,'limit':.03},{'id':'loaded-terminal-angle-refinement','passed':fine_voltage<.005,'value':fine_voltage,'limit':.005}])

result={'passed':all(c['passed'] for c in CHECKS),'settings':{'speed_rpm':1500,'frequency_Hz':50,'Rphase_ohm':RPHASE,'samples':COUNT,'model':'nonlinear M43, static synchronous snapshots; sinusoidal currents; fundamental resistive terminal balance','phase_basis_rad':PHI,'psi_PM_fundamental_peak_Wb':PSIPM,'Ld_local_H':LD,'Lq_local_H':LQ,'Br_T':BR,'mu_rec':MUREC},'loads':loads,'refinement':refinement,'runs':RUNS,'checks':CHECKS}
with open(OUT+'/Результаты.json','w') as f:json.dump(result,f)
with open(OUT+'/Нагрузочная-характеристика.csv','w') as f:
    keys=['Rload_ohm','Iphase_rms_A','Uphase_rms_V','Pload_W','Pcu_W','Tmean_Nm','Pmechanical_input_W','voltage_residual_relative','mechanical_power_residual_relative','Bparallel_min_T','Hparallel_min_Apm']
    f.write(';'.join(keys)+'\n')
    for row in loads:f.write(';'.join(str(row[k]) for k in keys)+'\n')
print('PM_LOAD_DONE',result['passed'],len(RUNS),[(r['Rload_ohm'],r['Uphase_rms_V'],r['mechanical_power_residual_relative']) for r in loads])