# 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-power-controls-native-result.json', 'directories': {'practice-ii03-3': '.', 'practice-ii10-1': '.'}, 'models': {'./practice-ii03-3-source-model.fem': {'problem': 'planar', 'depth': 40.0, 'angle': 25.0, 'circuits': ['Первичная обмотка', 'Вторичная обмотка']}, './strict-nonlinear-source-model.fem': {'problem': 'planar', 'depth': 40.0, 'angle': 25.0, 'circuits': ['Первичная обмотка', 'Вторичная обмотка']}, './practice-ii10-1-source-model.fem': {'problem': 'planar', 'depth': 40.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 c(v):return complex(v[0],v[1]) if isinstance(v,list) else complex(v)
try:
    with open('native-experiments-25.json') as f:newton=json.load(f)['experiments']['nonlinear_newton_validation']
    name='practice-ii03-3';start(name);path=model(name);names=CONFIG['models'][path]['circuits'];close();opendocument(path);mi_probdef(50,'millimeters','planar',1e-10,40,25);strict=os.getcwd()+'/strict-nonlinear-source-model.fem';mi_saveas(strict);close();CONFIG['models'][strict]=dict(CONFIG['models'][path]);x=newton['final_components_A'];r=solve(strict,'strict-nonlinear-final',[(names[0],complex(x[0],x[1])),(names[1],complex(x[2],x[3]))],material=(0,40,'Нелинейная сталь'));I=[complex(x[0],x[1]),complex(x[2],x[3])];r['residual_V']=[complex(r['circuits'][names[0]][1])+.8*I[0]-100,complex(r['circuits'][names[1]][1])+10*I[1]];r['precision']=1e-10;add('nonlinear_stricter_field_tolerance',r)
    name='practice-ii10-1';start(name);path=model(name);names=CONFIG['models'][path]['circuits']
    with open('native-experiments-02.json') as f:Z=[[c(z) for z in row] for row in json.load(f)['experiments']['transformer_linear']['Z_ohm']]
    rows=[]
    for T in (20.,100.):
        scale=1+.00393*(T-20);rp=(.8+.2)*scale;rs=(.1+.05)*scale
        for load in (5.,10.):
            matrix=[[Z[0][0]+rp,Z[0][1]],[Z[1][0],Z[1][1]+rs+load]];I=[complex(v) for v in np.linalg.solve(np.array(matrix),np.array([100+0j,0j]))];r=solve(path,'power-T'+str(T)+'-R'+str(load),list(zip(names,I)));v=[complex(r['circuits'][n][1]) for n in names];Pin=real(100*I[0].conjugate())/2;Pload=load*abs(I[1])**2/2;Pcu=(rp*abs(I[0])**2+rs*abs(I[1])**2)/2;r.update({'temperature_C':T,'load_ohm':load,'Rprimary_active_ohm':.8*scale,'Rprimary_end_ohm':.2*scale,'Rsecondary_active_ohm':.1*scale,'Rsecondary_end_ohm':.05*scale,'input_power_amplitude_W':Pin,'input_power_rms_W':real((100/math.sqrt(2))*(I[0]/math.sqrt(2)).conjugate()),'load_power_W':Pload,'external_copper_W':Pcu,'external_end_copper_W':(.2*scale*abs(I[0])**2+.05*scale*abs(I[1])**2)/2,'power_balance_W':Pin-Pload-Pcu-real(r['loss_total_W']),'native_component_sum_error_W':real(r['loss_total_W']-r['loss_magnetic_W']-r['loss_resistive_W']),'residual_V':[v[0]+rp*I[0]-100,v[1]+(rs+load)*I[1]],'scope':'Source coils have zero conductivity; no windingR is included in native field voltage. Explicit externalR fourpieces are assumptions at20C,not measured. Source transformer has no rotor,endrings,inverter,mechanicaloutput. All externalR added once.'});rows.append(r);add('power_temperature_two_loads',rows)
except Exception as error:
    DATA['fatal_error']=str(error);DATA['fatal_error_type']=type(error).__name__;save()
finally:close()
save()