# Execute with FEMM5 native embeddedMicroPython from this extracted folder.
import os
def _exists(p):
    try:os.stat(p);return True
    except OSError:return False
if _exists('files/reproduce.py'):os.chdir('files')
elif '__file__' in globals() and '/' in __file__.replace('\\','/'):os.chdir(__file__.replace('\\','/').rsplit('/',1)[0])

import os
_inputs=['gap-refined-source-model.fem', 'native-experiments-02.json', 'practice-ii05-1-source-model.fem', 'ring-effective-source-model.fem', 'uniformT10.0-source-model.fem', 'uniformT100.0-source-model.fem', 'uniformT110.0-source-model.fem', 'uniformT20.0-source-model.fem', 'uniformT30.0-source-model.fem', 'uniformT90.0-source-model.fem']
os.mkdir('rebuild-identification-output')
for _name in _inputs:
 with open(_name,'rb') as _f:_bytes=_f.read()
 with open('rebuild-identification-output/'+_name,'wb') as _f:_f.write(_bytes)
os.chdir('rebuild-identification-output')
import json,os,sys,time,math,cmath
import numpy as np
CONFIG={'output': 'rebuild-identification-native-result.json', 'directories': {'practice-ii05-1': '.'}, 'models': {'./practice-ii05-1-source-model.fem': {'problem': 'planar', 'depth': 100.0, 'angle': 20.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './gap-refined-source-model.fem': {'problem': 'planar', 'depth': 100.0, 'angle': 20.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './ring-effective-source-model.fem': {'problem': 'planar', 'depth': 100.0, 'angle': 20.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './uniformT10.0-source-model.fem': {'problem': 'planar', 'depth': 100.0, 'angle': 20.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './uniformT100.0-source-model.fem': {'problem': 'planar', 'depth': 100.0, 'angle': 20.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './uniformT110.0-source-model.fem': {'problem': 'planar', 'depth': 100.0, 'angle': 20.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './uniformT20.0-source-model.fem': {'problem': 'planar', 'depth': 100.0, 'angle': 20.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './uniformT30.0-source-model.fem': {'problem': 'planar', 'depth': 100.0, 'angle': 20.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './uniformT90.0-source-model.fem': {'problem': 'planar', 'depth': 100.0, 'angle': 20.0, '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 './'+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 c(v):return complex(v[0],v[1]) if isinstance(v,list) else complex(v)
try:
    name='practice-ii05-1';start(name);path=model(name);names=CONFIG['models'][path]['circuits']
    with open('native-experiments-02.json') as f:old=json.load(f)['experiments']['induction_frequencies']['rows']
    measured=[]
    for row in old:
        Z=sum(c(v[1])*c(v[0]).conjugate()/3 for v in row['circuits'].values());measured.append((row['frequency_Hz'],Z))
    training=[(f,Z) for f,Z in measured if f in (1.,3.)];heldout=[(f,Z) for f,Z in measured if f in (2.,4.)]
    def response(p,f):
        Ls,Lm,Rr=[math.exp(v) for v in p[:3]];R=Rr*(1+math.exp(p[3])*f*f) if len(p)==4 else Rr;w=2*math.pi*f;return 1j*w*Ls+1j*w*Lm*R/(R+1j*w*Lm)
    def fit(dim):
        p=[math.log(.005),math.log(.08),math.log(.5)]+([math.log(.0001)] if dim==4 else []);log=[]
        def residual(p):
            out=[]
            for f,Z in training:
                e=(response(p,f)-Z)/abs(Z);out.extend((e.real,e.imag))
            return out
        for iteration in range(30):
            r=residual(p);norm=sum(v*v for v in r);log.append({'iteration':iteration,'parameters':[math.exp(v) for v in p],'normalized_SSE':norm})
            if norm<1e-20:break
            h=1e-4;columns=[]
            for j in range(dim):
                pp=list(p);pm=list(p);pp[j]+=h;pm[j]-=h;rp=residual(pp);rm=residual(pm);columns.append([(rp[k]-rm[k])/(2*h) for k in range(4)])
            J=[[columns[j][i] for j in range(dim)] for i in range(4)];normal=[[sum(J[k][i]*J[k][j] for k in range(4))+(1e-14 if i==j else 0) for j in range(dim)] for i in range(dim)];rhs=[-sum(J[k][i]*r[k] for k in range(4)) for i in range(dim)];step=np.linalg.solve(np.array(normal),np.array(rhs));alpha=1.;improved=False
            while alpha>=1/64:
                candidate=[p[j]+alpha*float(step[j]) for j in range(dim)]
                if all(-20<v<5 for v in candidate) and sum(v*v for v in residual(candidate))<norm:p=candidate;improved=True;break
                alpha/=2
            if not improved or max(abs(float(v)*alpha) for v in step)<1e-8:break
        return {'parameters':{'L_leak_sector_H':math.exp(p[0]),'L_magnetizing_sector_H':math.exp(p[1]),'R_rotor_referred_sector_ohm':math.exp(p[2]),'R_frequency_coefficient_per_Hz2':math.exp(p[3]) if dim==4 else 0},'log':log,'training_Hz':[1,3],'heldout':[{'frequency_Hz':f,'field_Z_ohm':Z,'predicted_Z_ohm':response(p,f),'relative_prediction_error':abs(response(p,f)-Z)/abs(Z)} for f,Z in heldout],'physical_scope':'Low-frequency reduced rotor branch; positive Rr(f)=Rr0(1+k f²) is an explicit approximate low-frequency skin/distribution correction, validated only on held-out2/4Hz, not physical validation of the motor.'}
    constant=fit(3);corrected=fit(4)
    add('induction_identification',{'measured_positive_sequence_Z':measured,'constant_rotor_branch':constant,'frequency_corrected_rotor_branch':corrected,'leave_one_frequency_out_note':'A single complex frequency gives2real equations for3branch parameters, so the complete field branch is unidentifiable without another constraint. Do not disguise regularisation as unique identification. The separate analytic two-parameter RL example remains identifiable from one frequency.','analytic_RL':{'R_ohm':.5,'L_H':.02,'training_Hz':[1,3],'artificial_R_frequency_law_ohm':'0.5*(1+0.02*f_Hz**2)','constant_R_fit_ohm':.55,'residual_at_1Hz_ohm':-.04,'residual_at_3Hz_ohm':.04}})
    # Refine conductor mesh at held-out2Hz without recalibrating either model.
    points=[(.9,33.5),(7.3,32),(15,29.6),(20.2,25.8),(25.8,20.6),(29,14),(31.2,6.7),(31.6,.3)];mesh=[(x,y,'Алюминий',.5,'',0,4,1) for x,y in points];balanced=[cmath.exp(2j*math.pi*k/3) for k in range(3)];refined=solve(path,'heldout-2Hz-conductor-refined',list(zip(names,balanced)),2.,mesh=mesh);refined['positive_sequence_Z_ohm']=sum(complex(v[1])*complex(v[0]).conjugate()/3 for v in refined['circuits'].values());add('induction_heldout_mesh',refined)
    # Inspect the true annular air region before using Maxwell torque.
    solve(path,'fresh-ring-inspection',list(zip(names,balanced)),1.);mo_smooth('off');probes=[]
    for radius in (40.,40.125,40.25,40.375,40.5,40.625,40.75,40.875):
        for N in (256,1024):
            total=0.;nonair=0;mu_min=1e99;mu_max=0.
            for k in range(N):
                a=(k+.5)*(math.pi/2)/N;co=math.cos(a);si=math.sin(a);v=mo_getpointvalues(radius*co,radius*si);bx=complex(v[1]);by=complex(v[2]);mu=real(v[9]);mu_min=min(mu_min,mu);mu_max=max(mu_max,mu)
                if abs(mu-1)>1e-6:nonair+=1
                total+=real((bx*co+by*si)*(-bx*si+by*co).conjugate())*.5
            torque=4*.1*(radius*.001)**2/(4*math.pi*1e-7)*total*(math.pi/2)/N;probes.append({'radius_mm':radius,'N_quarter':N,'non_air_samples':nonair,'mu_relative_min':mu_min,'mu_relative_max':mu_max,'full_torque_Nm':torque,'smoothing':'off'})
    mo_smooth('on');add('induction_airgap_inspection',probes)
except Exception as error:
    DATA['fatal_error']=str(error);DATA['fatal_error_type']=type(error).__name__;save()
finally:close()
save()
print('NATIVE_INDUCTION_VALIDATION_COMPLETE',DATA['field_solutions'])

assert 'fatal_error' not in DATA,DATA.get('fatal_error')