# Run inside native FEMM5. Every path is relative to this supplied folder.
import os
def _path_exists(path):
    try:
        os.stat(path)
        return True
    except OSError:
        return False
if _path_exists('files/reproduce.py'):
    os.chdir('files')
elif '__file__' in globals() and '/' in __file__.replace('\\','/'):
    os.chdir(__file__.replace('\\','/').rsplit('/',1)[0])
import json,os,sys,time,math,cmath
import numpy as np
CONFIG={'output': 'replay-native-result.json', 'directories': {'practice-ii08-1': '.'}, 'models': {'./compact-source-model.fem': {'problem': 'planar', 'depth': 80, 'angle': 25, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}}}
DATA={'runtime':sys.version,'implementation':str(sys.implementation),'experiments':{},'field_solutions':0}
assert sys.implementation.name=='micropython'
def save():
    with open(CONFIG['output'],'w') as f:json.dump(DATA,f)
def serial(v):
    if isinstance(v,complex):return [v.real,v.imag]
    if isinstance(v,(tuple,list)):return [serial(x) for x in v]
    if isinstance(v,dict):return {k:serial(x) for k,x in v.items()}
    return v
def close():
    for n in ('mo_close','mi_close'):
        try:globals()[n]()
        except Exception:pass
def start(name):
    close();d=CONFIG['directories'][name];os.chdir(d);return d
def model(name):return CONFIG['directories'][name]+'/source-model.fem'
def solve(path,tag,currents=None,frequency=None,rotate=0,group=2,mesh=None,material=None):
    close();opendocument(path);mi_saveas(os.getcwd()+'/'+tag+'.fem')
    if frequency is not None:
        info=CONFIG['models'][path]
        mi_probdef(frequency,'millimeters',info['problem'],1e-8,info['depth'],info['angle'])
    if rotate:
        mi_selectgroup(group);mi_moverotate(0,0,rotate);mi_clearselected()
    if material:
        x,y,material_name=material;mi_selectlabel(x,y);mi_setblockprop(material_name,0,1,'',0,0,0);mi_clearselected()
    if mesh:
        for x,y,mat,sz,circ,mag,g,turns in mesh:
            mi_selectlabel(x,y);mi_setblockprop(mat,0,sz,circ,mag,g,turns);mi_clearselected()
    if currents:
        for name,value in currents:mi_modifycircprop(name,1,value)
    mi_analyze(1);mi_loadsolution();DATA['field_solutions']+=1
    info={'file':os.getcwd()+'/'+tag+'.fem','nodes':mo_numnodes(),'elements':mo_numelements(),'circuits':{}}
    for name in CONFIG['models'][path]['circuits']:
        info['circuits'][name]=list(mo_getcircuitproperties(name))
    mo_clearblock();mo_groupselectblock();info['energy_J']=mo_blockintegral(2);info['coenergy_J']=mo_blockintegral(17)
    info['loss_magnetic_W']=mo_blockintegral(3);info['loss_resistive_W']=mo_blockintegral(4);info['loss_total_W']=mo_blockintegral(6)
    mo_clearblock()
    try:
        mo_groupselectblock(group);info['torque_Nm']=mo_blockintegral(22);info['force_x_N']=mo_blockintegral(18);info['force_y_N']=mo_blockintegral(19);mo_clearblock()
    except Exception as error:info['force_error']=str(error)
    return info
def real(x):return complex(x).real
def add(name,value):DATA['experiments'][name]=serial(value);save()
def z_matrix(path,prefix,names):
    columns=[];evidence=[]
    for j in range(len(names)):
        row=solve(path,prefix+'-basis-'+str(j),[(n,1 if j==k else 0) for k,n in enumerate(names)])
        columns.append([complex(row['circuits'][n][1]) for n in names]);evidence.append(row)
    return [[columns[j][i] for j in range(len(names))] for i in range(len(names))],evidence
def linear_load(path,prefix,names,Z,rloads,rprimary,source,neutral=False):
    n=len(names);m=n+1 if neutral else n
    matrix=[[0j for j in range(m)] for i in range(m)]
    rhs=[0j]*m
    for i in range(n):
        for j in range(n):matrix[i][j]=Z[i][j]
        if i<n//2:matrix[i][i]+=rprimary;rhs[i]=source[i]
        else:
            matrix[i][i]+=rloads[i-n//2]
            if neutral:matrix[i][n]=1
    if neutral:
        for j in range(n//2,n):matrix[n][j]=1
    answer=np.linalg.solve(np.array(matrix),np.array(rhs));currents=[complex(answer[j]) for j in range(n)]
    row=solve(path,prefix,[(n,currents[j]) for j,n in enumerate(names)])
    voltages=[complex(row['circuits'][n][1]) for n in names]
    residual=[]
    for i in range(n):
        r=voltages[i]+(rprimary if i<n//2 else rloads[i-n//2])*currents[i]-(source[i] if i<n//2 else 0)
        if neutral and i>=n//2:r+=complex(answer[n])
        residual.append(r)
    row['current_A']=currents;row['voltage_V']=voltages;row['residual_V']=residual
    row['source_power_W']=sum(real(source[i]*currents[i].conjugate())/2 for i in range(n//2))
    row['external_primary_copper_W']=sum(rprimary*abs(currents[i])**2/2 for i in range(n//2))
    row['load_power_W']=sum(rloads[i-n//2]*abs(currents[i])**2/2 for i in range(n//2,n))
    row['neutral_current_A']=-sum(currents[n//2:]);row['neutral_potential_V']=complex(answer[n]) if neutral else 0
    row['power_balance_W']=row['source_power_W']-row['external_primary_copper_W']-row['load_power_W']-real(row['loss_total_W'])
    return row

def z(v):return complex(v[0],v[1]).real if isinstance(v,list) else real(v)
try:
    name='practice-ii08-1';start(name)
    with open('source-map.json') as f:source=json.load(f)['experiments']
    table=source['srm_map'];assert table['status']=='complete'
    lookup={(r['mechanical_angle_deg'],r['requested_current_A']):r for r in table['rows']}
    all_currents=table['currents_A'];fine_angles=table['angles_deg']
    def linkage(angle,current):return z(lookup[(angle,current)]['circuits']['Фаза А'][2])
    derived=[]
    for angle in fine_angles:
        q={'angle_deg':angle,'currents':[]}
        for current in (2.,8.):
            psi=linkage(angle,current);steps=[]
            for h in (.1,.05):steps.append({'h_A':h,'L_dynamic_H':(linkage(angle,current+h)-linkage(angle,current-h))/(2*h)})
            q['currents'].append({'I_A':current,'psi_Wb':psi,'L_static_H':psi/current,'steps':steps})
        derived.append(q)
    add('srm_inductance_curves',{'rows':derived,'source_period_deg':90,'source_rotor_shift_deg':-15,'definition':'L_static=Psi/I (no PM); L_dynamic=central dPsi/dI, fixed theta; all original off-grid perturbations independently solved.'})
    # Interpolation and energy are derived from our OWN complete FEMM Psi map.
    class Map:
        def __init__(self,coarse=False):
            self.angles=[5.*k for k in range(19)] if coarse else fine_angles
            self.currents=[0.,2.,4.,6.,8.,10.] if coarse else all_currents
        def interval(self,current,theta):
            if current<-1e-10 or current>10.+1e-10:raise ValueError('Current outside measured [0,10] A; extrapolation forbidden')
            current=max(0.,min(10.,current));theta=theta%90
            a=0
            while a<len(self.angles)-2 and theta>=self.angles[a+1]:a+=1
            j=0
            while j<len(self.currents)-2 and current>=self.currents[j+1]:j+=1
            lo,hi=self.angles[a],self.angles[a+1];t=(theta-lo)/(hi-lo)
            return current,lo,hi,t,j
        def at_angle(self,current,angle,j):
            values=[linkage(angle,i) for i in self.currents];W=0.
            for k in range(j):W+=(values[k]+values[k+1])*.5*(self.currents[k+1]-self.currents[k])
            a,b=self.currents[j],self.currents[j+1];psi=values[j]+(values[j+1]-values[j])*(current-a)/(b-a);W+=(values[j]+psi)*.5*(current-a)
            return psi,W
        def state(self,current,theta):
            current,a,b,t,j=self.interval(current,theta);pa,wa=self.at_angle(current,a,j);pb,wb=self.at_angle(current,b,j);psi=pa*(1-t)+pb*t;coenergy=wa*(1-t)+wb*t
            return {'psi_Wb':psi,'coenergy_J':coenergy,'stored_energy_J':current*psi-coenergy,'torque_Nm':(wb-wa)/math.radians(b-a)}
        def inverse(self,psi,theta):
            if psi<-1e-10:raise ValueError('Negative flux for one-directional phase')
            if psi>self.state(10.,theta)['psi_Wb']+1e-9:raise ValueError('Flux above measured map at current angle; no extrapolation')
            lo,hi=0.,10.
            for _ in range(32):
                mid=(lo+hi)/2
                if self.state(mid,theta)['psi_Wb']<psi:lo=mid
                else:hi=mid
            return (lo+hi)/2
    fine=Map(False);coarse=Map(True);map_tests=[]
    for current in (2.,8.):
        a=22.5;h=2.5;minus=z(lookup[(a-h,current)]['coenergy_J']);plus=z(lookup[(a+h,current)]['coenergy_J']);torque=(plus-minus)/(2*math.radians(h));field=z(lookup[(a,current)]['torque_Nm']);map_tests.append({'I_A':current,'angle_deg':a,'coenergy_derivative_torque_Nm':torque,'FEMM_WST_torque_Nm':field,'relative_difference':abs(torque-field)/max(abs(field),.01)})
    coupling=source['srm_two_phase'];single=lookup[(22.5,2.)];add('srm_map_checks',{'coenergy_torque_comparison':map_tests,'two_phase_control':coupling,'single_phase_linkage_Wb':linkage(22.5,2.),'two_phase_linkage_Wb':z(coupling['circuits']['Фаза А'][2]),'phase_model_scope':'One active phase only; no unsupported independent-phase superposition during overlap.','coarse_grid':{'angles_deg':coarse.angles,'currents_A':coarse.currents},'fine_grid':{'angles_deg':fine.angles,'currents_A':fine.currents}})
    # An analytically verified RL integrator, independent of the FEMM map.
    rl=[]
    for dt in (.0002,.0001,.00005):
        i=0.;steps=int(.01/dt)
        for _ in range(steps):k1=(24-2*i)/.01;predict=i+dt*k1;k2=(24-2*predict)/.01;i+=dt*(k1+k2)/2
        exact=12*(1-math.exp(-2*.01/.01));rl.append({'dt_s':dt,'current_A':i,'exact_A':exact,'absolute_error_A':abs(i-exact)})
    add('RL_integrator',rl)
    def simulate(map,dt,off_angle,load=None,duration=.032):
        R=2.;U=24.;J=.001;t=0.;theta=-20.;omega=100.;psi=0.;I=0.;source_energy=0.;copper=0.;mechanical=0.;returned=0.;rows=[];state='ON';initial_kinetic=.5*J*omega*omega
        while t<duration-1e-12:
            position=(theta+20)%90;off=off_angle+20
            if position<off-1e-9:
                if I>=3.:state='FREEWHEEL'
                elif I<=2.9:state='ON'
                if state not in ('ON','FREEWHEEL'):state='ON'
                u=U if state=='ON' else 0.
            else:state='RETURN' if psi>1e-10 else 'OPEN';u=-U if state=='RETURN' else 0.
            h=min(dt,duration-t)
            # Land exactly on angular on/off events for prescribed-speed pulses;
            # for the coupled model, the correction is accurate to the Heun order.
            event=off-position if position<off-1e-9 else 90-position
            if omega>0 and event>1e-10:h=min(h,math.radians(event)/omega)
            old=map.state(I,theta);rate=u-R*I
            if rate<0 and psi+h*rate<0:h=min(h,max(1e-10,psi/(-rate)))
            accel=(old['torque_Nm']-load)/J if load is not None else 0.
            omega_predict=omega+h*accel;theta_predict=theta+math.degrees(h*(omega+omega_predict)/2)
            psi_predict=max(0.,psi+h*rate);ip=map.inverse(psi_predict,theta_predict) if psi_predict>1e-12 else 0.
            psi_next=max(0.,psi+h*((u-R*I)+(u-R*ip))/2)
            i_next=map.inverse(psi_next,theta_predict) if psi_next>1e-12 else 0.;new=map.state(i_next,theta_predict)
            avgI=(I+i_next)/2;energy=u*avgI*h;source_energy+=energy;copper+=R*(I*I+i_next*i_next)*h/2;dtheta=math.radians(theta_predict-theta);mechanical+=(old['torque_Nm']+new['torque_Nm'])*.5*dtheta
            if u<0:returned-=energy
            if load is not None:omega+=h*((old['torque_Nm']-load)+(new['torque_Nm']-load))/(2*J)
            t+=h;theta=theta_predict;psi=psi_next;I=i_next
            magnetic=new['stored_energy_J'];balance=source_energy-copper-magnetic-mechanical
            rows.append({'time_s':t,'mechanical_angle_deg':theta,'current_A':I,'linkage_Wb':psi,'voltage_V':u,'converter_state':state,'speed_rad_s':omega,'torque_Nm':new['torque_Nm'],'source_energy_J':source_energy,'copper_energy_J':copper,'magnetic_energy_J':magnetic,'mechanical_work_J':mechanical,'energy_residual_J':balance})
        final=rows[-1];kinetic=.5*J*omega*omega-initial_kinetic;load_work=load*math.radians(theta+20) if load is not None else None
        return {'dt_s':dt,'off_angle_deg':off_angle,'load_torque_Nm':load,'duration_s':duration,'rows':rows,'peak_current_A':max(r['current_A'] for r in rows),'returned_energy_J':returned,'final_speed_rad_s':omega,'energy_residual_J':final['energy_residual_J'],'relative_energy_residual':abs(final['energy_residual_J'])/max(abs(source_energy),abs(mechanical),.001),'mechanical_energy_residual_J':mechanical-kinetic-load_work if load is not None else None,'limits':'Ideal one-phase asymmetric bridge; specified24V,2ohm; hysteresis2.9/3A; return voltage−24V. No switching losses or transient eddy-current field memory.'}
    pulses=[]
    for off in (5.,10.):
        for dt in (.0001,.00005):pulses.append({'map':'fine','simulation':simulate(fine,dt,off)})
        pulses.append({'map':'coarse','simulation':simulate(coarse,.00005,off)})
    add('srm_commutation',pulses)
    mechanics=[]
    for load in (.01,.05):
        for dt in (.0001,.00005):mechanics.append({'map':'fine','simulation':simulate(fine,dt,10.,load,.05)})
        mechanics.append({'map':'coarse','simulation':simulate(coarse,.00005,10.,load,.05)})
    add('srm_mechanics',mechanics)
    # Independent fresh fields at three states taken from this coupled simulation.
    chosen=mechanics[1]['simulation']['rows'];active=[r for r in chosen if r['current_A']>.5];replays=[];path=CONFIG['directories'][name]+'/compact-source-model.fem'
    for k,index in enumerate((len(active)//5,len(active)//2,4*len(active)//5)):
        state=active[index];r=solve(path,'dynamic-state-'+str(k),[('Фаза А',state['current_A']),('Фаза В',0),('Фаза С',0)],rotate=state['mechanical_angle_deg']%90);actual=z(r['circuits']['Фаза А'][2]);replays.append({'state':state,'field':r,'field_linkage_Wb':actual,'interpolated_linkage_Wb':state['linkage_Wb'],'relative_interpolation_error':abs(actual-state['linkage_Wb'])/max(abs(actual),1e-9)})
    add('srm_dynamic_field_replay',replays)
except Exception as error:
    DATA['fatal_error']=str(error);DATA['fatal_error_type']=type(error).__name__;save()
finally:close()
save()
print('NATIVE_SRM_DYNAMICS_COMPLETE',DATA['field_solutions'])

# Strict portable completion receipt; process exit alone never proves completion.
_checks=[]
def _ck(label,actual,expected,tolerance):
    ok=abs(actual-expected)<=tolerance
    _checks.append({'label':label,'actual':actual,'expected':expected,'tolerance':tolerance,'passed':ok})
_ck('native controller has no fatal_error',int('fatal_error' in DATA),0,0)
_ck('three independent fresh field states',DATA.get('field_solutions',0),3,0)
_e=DATA.get('experiments',{})
_field=_e.get('srm_dynamic_field_replay',[])
_ck('three field controls recorded',len(_field),3,0)
for _j,_r in enumerate(_field):_ck('field interpolation '+str(_j),_r['relative_interpolation_error'],0,.01)
_rl=_e.get('RL_integrator',[])
_ck('three exactRL step controls',len(_rl),3,0)
_ck('RL errors decrease under step halving',int(len(_rl)==3 and _rl[2]['absolute_error_A']<_rl[1]['absolute_error_A']<_rl[0]['absolute_error_A']),1,0)
for _key in ('srm_commutation','srm_mechanics'):
    _series=_e.get(_key,[])
    _fine=[_r['simulation'] for _r in _series if _r['map']=='fine']
    _ck(_key+' two cases at two time steps',len(_fine),4,0)
    for _j,_s in enumerate(_fine):
        _ck(_key+' energy residual '+str(_j),_s['relative_energy_residual'],0,.02 if _s['dt_s']>.00005 else .01)
    _steps=sorted(set(_s['dt_s'] for _s in _fine))
    _ck(_key+' two distinct time steps',len(_steps),2,0)
for _r in _e.get('srm_map_checks',{}).get('coenergy_torque_comparison',[]):_ck('independent coenergy/WST torque '+str(_r['I_A']),_r['relative_difference'],0,.03)
_passed=all(_q['passed'] for _q in _checks)
with open('replay-result.json','w') as _f:json.dump({'passed':_passed,'checks':_checks,'nativeFreshFieldCount':DATA.get('field_solutions',0)},_f)
assert _passed,'Portable SRM replay failed a numerical or completeness check'