# 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': 'native-steel-loss-map-native-result.json', 'directories': {'practice-ii07-4': '.'}, 'models': {'./practice-ii07-4-source-model.fem': {'problem': 'planar', 'depth': 80.0, 'angle': 25.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './N35SH-150C-magmesh-0.4-source-model.fem': {'problem': 'planar', 'depth': 80.0, 'angle': 25.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './N35SH-150C-magmesh-0.8-source-model.fem': {'problem': 'planar', 'depth': 80.0, 'angle': 25.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './N35SH-20C-magmesh-0.4-source-model.fem': {'problem': 'planar', 'depth': 80.0, 'angle': 25.0, 'circuits': ['Фаза А', 'Фаза В', 'Фаза С']}, './N35SH-20C-magmesh-0.8-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

try:
    name='practice-ii07-4';start(name)
    with open(CONFIG['directories'][name]+'/materialized-models.json') as f:materials=json.load(f)
    spec=next(s for s in materials['models'] if s['temperature_C']==150 and s['magnet_mesh_mm']==.4);path=spec['path'].replace('\\','/');names=CONFIG['models'][path]['circuits'];N=48;alpha=math.log(2.50/1.05)/math.log(1.5);rows=[]
    for regime,id_,iq in [('nominal',0.,5.),('demagnetizing',-30.,0.)]:
        cells=[];acc={n:[] for n in (12,24,48)};fields=[];probes=[]
        for k in range(N):
            theta=180*k/N;phi=math.radians(90+2*theta);currents=[(n,id_*math.cos(phi-2*math.pi*j/3)-iq*math.sin(phi-2*math.pi*j/3)) for j,n in enumerate(names)];r=solve(path,'loss-angle-transient',currents,rotate=theta);r.pop('file');r['angle_deg']=theta;fields.append(r);mo_smooth('off')
            if not cells:
                for index in range(1,int(mo_numelements())+1):
                    el=mo_getelement(index);x,y,area=float(el[3]),float(el[4]),float(el[5]);v=mo_getpointvalues(x,y)
                    if min(real(v[9]),real(v[10]))>2:
                        rotor=int(el[6])==2;cells.append([x,y,area,rotor]);
                        for n in (12,24,48):acc[n].append([0.,0.,0.,0.,0.,0.])
                DATA['experiments']['steel_loss_progress']={'regime':regime,'cells':len(cells),'angle_index':k};save()
            angle=math.radians(theta);co,si=math.cos(angle),math.sin(angle);phase=2*angle;cc,ss=math.cos(phase),math.sin(phase);probe=[]
            for j,(x,y,area,rotor) in enumerate(cells):
                px,py=(x*co-y*si,x*si+y*co) if rotor else (x,y);v=mo_getpointvalues(px,py);bx,by=real(v[1]),real(v[2])
                if rotor:bx,by=bx*co+by*si,-bx*si+by*co
                for n in (12,24,48):
                    if k%(48//n)==0:
                        a=acc[n][j];a[0]+=bx*cc;a[1]+=bx*ss;a[2]+=by*cc;a[3]+=by*ss;a[4]+=bx;a[5]+=by
                if j%max(1,len(cells)//12)==0:probe.append({'cell':j,'x_mm':px,'y_mm':py,'Bx_T':bx,'By_T':by,'rotor_frame':rotor})
            probes.append({'angle_deg':theta,'time_s':k/(N*50),'probes':probe});mo_smooth('on')
            DATA['experiments']['steel_loss_progress']={'regime':regime,'cells':len(cells),'angle_index':k};save()
        estimates=[]
        for n in (12,24,48):
            output=[];mass_total=0.;mass_outside=0.;power_total=0.;power_outside=0.;rotor_power=0.
            for j,(x,y,area,rotor) in enumerate(cells):
                a=acc[n][j];cx,sx,cy,sy=[2*v/n for v in a[:4]];u=cx*cx+cy*cy;v=sx*sx+sy*sy;dot=cx*sx+cy*sy;peak=math.sqrt(max(0,(u+v+math.sqrt((u-v)**2+4*dot*dot))/2));mass=area*1e-6*.08*7650;loss=1.05*peak**alpha*mass;inside=1<=peak<=1.5
                mass_total+=mass;power_total+=loss
                if not inside:mass_outside+=mass;power_outside+=loss
                if rotor:rotor_power+=loss
                output.append({'x_reference_mm':x,'y_reference_mm':y,'area_mm2':area,'mass_kg':mass,'rotor':rotor,'B50_major_peak_T':peak,'DC_Bx_T':a[4]/n,'DC_By_T':a[5]/n,'P50_estimate_W':loss,'inside_calibration_Brange':inside})
            estimates.append({'N_angles':n,'steel_mass_kg':mass_total,'P50_estimate_W':power_total,'P50_rotor_estimate_W':rotor_power,'mass_outside_calibration_kg':mass_outside,'mass_outside_fraction':mass_outside/mass_total,'P_outside_calibration_W':power_outside,'P_outside_fraction':power_outside/power_total,'cells':output if n==48 else []})
        rows.append({'regime':regime,'id_A':id_,'iq_A':iq,'temperature_C':150,'mechanical_speed_rpm':1500,'electrical_fundamental_Hz':50,'fields':fields,'probes':probes,'estimates':estimates});add('steel_loss_estimate',{'rows':rows,'status':'complete' if len(rows)==2 else 'running','law':{'grade':'voestalpine isovac270-35A','source':'official January2026 product range PDF page2, longitudinal+transverse typical50Hz DINEN60404-2','points_B_T_p_W_kg':[[1,1.05],[1.5,2.50]],'density_kg_m3':7650,'alpha_derived':alpha,'formula':'p50=1.05*B50_major_peak**alpha W/kg','calibration_B_T':[1,1.5]},'limits':'Conditional50Hz fundamental alternating-loss estimate on course8pointBH field, not an isovacBH identity or validated whole-machine loss. Outside1..1.5T power-law extrapolation reported separately. Rotating/elliptical magnetization approximated by major-axis alternating amplitude; DC, higher harmonics, excess losses not separately calibrated; magnet eddy currents, mechanical losses and temperature-dependent steel loss absent. Co-rotating rotor coordinates prevent interpreting DC rotor flux as50Hz flux. NativeDC loss commands separately return0, not this postprocessed estimate.'})
except Exception as error:
    DATA['fatal_error']=str(error);DATA['fatal_error_type']=type(error).__name__;save()
finally:close()
save()