# 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 json,os,sys,time,math,cmath
import numpy as np
CONFIG={'output': 'postprocess-full-native-result.json', 'directories': {'practice-ii01-2': '.'}, 'models': {'./practice-ii01-2-source-model.fem': {'problem': 'planar', 'depth': 80.0, 'angle': 25.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './sector-source-model.fem': {'problem': 'planar', 'depth': 80.0, 'angle': 25.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 read(attempt):
    with open('native-experiments-'+attempt+'.json') as f:return json.load(f)['experiments']
def z(v):return complex(v[0],v[1]) if isinstance(v,list) else complex(v)
def spectrum(values,harmonics):
    N=len(values);F=np.fft.fft(np.array(values));return {'N':N,'mean':real(F[0])/N,'harmonics':[{'order':k,'amplitude':2*abs(complex(F[k]))/N,'phase_rad':cmath.phase(complex(F[k]))} for k in harmonics]}
try:
    angles=read('03');outputs={}
    for assignment in ('practice-ii01-2','practice-ii07-3'):
        original=angles[assignment]['rows'];table=[]
        counts=(48,96,192,128,256,768) if assignment=='practice-ii01-2' else (128,256)
        for N in counts:
            stride=len(original)//N;rows=original[::stride]
            assert len(rows)==N
            names=list(rows[0]['circuits']);phase=[]
            for name in names:
                linkage=spectrum([real(z(r['circuits'][name][2])) for r in rows],[1,3,5,7,11])
                for harmonic in linkage['harmonics']:
                    k=harmonic['order'];harmonic['E_rms_50Hz_V']=2*math.pi*50*k*harmonic['amplitude']/math.sqrt(2);harmonic['E_rms_100Hz_V']=2*math.pi*100*k*harmonic['amplitude']/math.sqrt(2)
                phase.append({'name':name,'linkage_Wb':linkage})
            torque=[real(z(r['torque_Nm'])) for r in rows]
            table.append({'N':N,'stride':stride,'first_angle_deg':rows[0]['mechanical_angle_deg'],'last_angle_deg':rows[-1]['mechanical_angle_deg'],'repeated_endpoint':False,'unique_count':len(set(r['mechanical_angle_deg'] for r in rows)),'phases':phase,'torque_Nm':spectrum(torque,[1,6,12,18]),'torque_peak_to_peak_Nm':max(torque)-min(torque)})
        outputs[assignment]=table
    add('field_spectra',outputs)
    add('periods',{'actual_source_models':{'pole_pairs':2,'slots':12,'flux_linkage_mechanical_period_deg':180.,'cogging_period_deg':30.},'separate_assignment_Q12_2p10':{'pole_pairs':5,'slots':12,'lcm_slots_poles':60,'cogging_period_deg':6.,'N_cogging':64,'mechanical_step_deg':6/64,'phase_emf_period_deg':72.,'phase_samples_for_same_step':768},'separate_intro_example_p4_speed1500':{'pole_pairs':4,'speed_rpm':1500,'electrical_frequency_Hz':100,'flux_linkage_period_deg':90},'actual_spm_at_1500rpm':{'pole_pairs':2,'speed_rpm':1500,'electrical_frequency_Hz':50,'flux_linkage_period_deg':180}})
    # Nonlinear matrix values and independent derivative/interpolation validation.
    matrix=read('02')['nonlinear_matrix'];points=[]
    for p in matrix['points']:
        steps=[]
        for s in p['steps']:
            L=s['L_H'];trace=L[0][0]+L[1][1];disc=math.sqrt((L[0][0]-L[1][1])**2+4*L[0][1]*L[1][0]);steps.append({'h_A':s['h_A'],'L_H':L,'eigenvalues_H':[(trace-disc)/2,(trace+disc)/2],'reciprocity_abs_H':abs(L[0][1]-L[1][0]),'reciprocity_relative':abs(L[0][1]-L[1][0])/max(abs(L[0][1]),abs(L[1][0]),1e-20)})
        a,b=steps[0]['L_H'],steps[1]['L_H'];points.append({'base_A':p['base_A'],'steps':steps,'step_difference_relative':max(abs(a[i][j]-b[i][j])/max(abs(b[i][j]),1e-12) for i in range(2) for j in range(2))})
    def interpolate(i1,i2):
        if not (0<=i1<=4 and 0<=i2<=2):raise ValueError('Ток вне рассчитанной карты [0,4]×[0,2]A')
        j=min(1,int(i1/2));k=min(1,int(i2));t=(i1-2*j)/2;u=i2-k;weights=[((j,k),(1-t)*(1-u)),((j+1,k),t*(1-u)),((j,k+1),(1-t)*u),((j+1,k+1),t*u)];names=list(matrix['map'][0]['circuits']);return [sum(real(z(matrix['map'][a*3+b]['circuits'][name][2]))*w for (a,b),w in weights) for name in names]
    try:interpolate(5,1);outside='FAILED'
    except ValueError as error:outside=str(error)
    add('nonlinear_matrix_checks',{'points':points,'interpolation_at_1A_0p5A_Wb':interpolate(1,.5),'out_of_map_error':outside,'reciprocity_policy':'Original off-diagonal values retained; no symmetrisation hides their discrepancy.'})
    mesh=read('15')['mesh_angular_series'];convergence=[]
    for group in mesh:
        rows=group['rows'];T=[real(z(r['torque_Nm'])) for r in rows];names=list(rows[0]['circuits']);flux={n:spectrum([real(z(r['circuits'][n][2])) for r in rows],[1]) for n in names};convergence.append({'h_mm':group['h_mm'],'angle_count':len(rows),'mean_torque_Nm':sum(T)/len(T),'torque_peak_to_peak_Nm':max(T)-min(T),'phase_linkage_spectra_Wb':flux,'elements_min':min(r['elements'] for r in rows),'elements_max':max(r['elements'] for r in rows)})
    add('spm_mesh_convergence',{'rows':convergence,'predefined_phase_for_flux':'Фаза В: source reference |PsiB|≈0.29Wb, while PsiA≈0.0018Wb is a small cancellation quantity; use absolute1e-4Wb for the latter.','criteria':{'phase_flux_fundamental_relative':.01,'mean_torque_relative':.02},'single_point_external_domain':read('12')['spm_external_domain']})
except Exception as error:
    DATA['fatal_error']=str(error);DATA['fatal_error_type']=type(error).__name__;save()
finally:close()
save()
print('NATIVE_POSTPROCESS_COMPLETE')