Дифференциальная пара: матрица ёмкостей и квазистатическое сопротивление

"""Дифференциальная пара: матрица ёмкостей и квазистатическое сопротивление.
Среда: встроенный MicroPython FEMM 5; команды FEMM доступны глобально.
Геометрические размеры: миллиметры; поля и интегралы возвращаются в СИ.
Все материалы и проводники имеют русские имена. Каждая модель сохраняется,
закрывается и повторно открывается до нового расчёта, включая кириллицу.
Запуск создаёт новую папку calculated-NN рядом с программой; опубликованные
модели и предыдущие расчёты не перезаписываются. Все необходимые
исходные данные включены. Условия задачи и ограничения описаны в README.md.
"""
import os
import sys
import math
import json

saved = globals().get('__file__', '') or (sys.argv[0] if sys.argv else '')
base = os.path.dirname(os.path.abspath(saved)) if saved else os.getcwd()
index = 1
while os.path.exists(base + '/calculated-%02d' % index):
    index += 1
directory = base + '/calculated-%02d' % index
os.mkdir(directory)
records = []
checks = []
EPS0 = 8.85418781762e-12


def real(value):
    if isinstance(value, (list, tuple)):
        return real(value[0])
    return complex(value).real


def check(name, value, limit):
    checks.append({'name': name, 'value': value, 'limit': limit,
                   'passed': math.isfinite(value) and value <= limit})


def label(x, y, material, mesh, group=1):
    ei_addblocklabel(x, y)
    ei_selectlabel(x, y)
    ei_setblockprop(material, 0, mesh, group)
    ei_clearselected()


def circle(x, radius, conductor, angle):
    ei_addnode(x - radius, 0)
    ei_addnode(x + radius, 0)
    ei_addarc(x - radius, 0, x + radius, 0, 180, angle)
    ei_addarc(x + radius, 0, x - radius, 0, 180, angle)
    ei_selectarcsegment(x, radius)
    ei_selectarcsegment(x, -radius)
    ei_setarcsegmentprop(angle, '<None>', 0, 0, conductor)
    ei_clearselected()


def solve(model):
    ei_saveas(model)
    ei_close()
    opendocument(model)
    ei_analyze(1)
    ei_loadsolution()
    eo_smooth('off')


def electric_result(mesh, model, input_name, output_name, exact_C):
    voltage_in, charge_in = eo_getconductorproperties(input_name)
    voltage_out, charge_out = eo_getconductorproperties(output_name)
    voltage = real(voltage_in) - real(voltage_out)
    charge = real(charge_in)
    eo_groupselectblock(1)
    energy = real(eo_blockintegral(0))
    eo_clearblock()
    capacitance = charge / voltage
    energy_C = 2 * energy / voltage ** 2
    row = {'mesh_mm': mesh, 'model': os.path.basename(model),
           'nodes': eo_numnodes(), 'elements': eo_numelements(),
           'voltage_V': voltage, 'charge_C': charge, 'energy_J': energy,
           'capacitance_F': capacitance, 'energy_capacitance_F': energy_C,
           'analytic_capacitance_F': exact_C,
           'analytic_C_relative_error': abs(capacitance / exact_C - 1),
           'energy_charge_relative_error': abs(energy_C / capacitance - 1)}
    records.append(row)
    check('Энергия и заряд: ' + os.path.basename(model),
          row['energy_charge_relative_error'], 1e-5)
    return row


def circle_xy(cx, cy, radius, conductor, angle):
    ei_addnode(cx - radius, cy)
    ei_addnode(cx + radius, cy)
    ei_addarc(cx - radius, cy, cx + radius, cy, 180, angle)
    ei_addarc(cx + radius, cy, cx - radius, cy, 180, angle)
    ei_selectarcsegment(cx, cy + radius)
    ei_selectarcsegment(cx, cy - radius)
    ei_setarcsegmentprop(angle, '<None>', 0, 0, conductor)
    ei_clearselected()


def polygon(points, conductor='<None>', mesh=0):
    for x, y in points:
        ei_addnode(x, y)
    for k, (x, y) in enumerate(points):
        x2, y2 = points[(k + 1) % len(points)]
        ei_addsegment(x, y, x2, y2)
        if conductor != '<None>':
            ei_selectsegment((x + x2) / 2, (y + y2) / 2)
            ei_setsegmentprop('<None>', mesh, 0, 0, 0, conductor)
            ei_clearselected()


def terminal(x, y, name, mesh):
    ei_selectsegment(x, y)
    ei_setsegmentprop('<None>', mesh, 0, 0, 0, name)
    ei_clearselected()


def summary(mesh, model, source, ground, exact_C=None):
    v1, q1 = eo_getconductorproperties(source)
    v0, q0 = eo_getconductorproperties(ground)
    v1, v0, q1, q0 = real(v1), real(v0), real(q1), real(q0)
    voltage = v1 - v0
    eo_groupselectblock(1)
    energy = real(eo_blockintegral(0))
    eo_clearblock()
    capacitance = q1 / voltage
    row = {'mesh_mm': mesh, 'model': os.path.basename(model),
           'nodes': eo_numnodes(), 'elements': eo_numelements(),
           'voltage_V': voltage, 'charge_C': q1, 'energy_J': energy,
           'capacitance_F': capacitance,
           'energy_capacitance_F': 2 * energy / voltage ** 2}
    if exact_C is not None:
        row['analytic_capacitance_F'] = exact_C
        row['analytic_C_relative_error'] = abs(capacitance / exact_C - 1)
    records.append(row)
    check('Энергия и терминальный заряд: ' + os.path.basename(model),
          abs(2 * energy / voltage ** 2 / capacitance - 1), 1e-5)
    return row


def field(x, y):
    values = eo_getpointvalues(x, y)
    return math.sqrt(real(values[3]) ** 2 + real(values[4]) ** 2)


def rectangle_field(mesh, filename, epsx=2, epsy=2, voltage=100, width=10, height=10, depth=100):
    newdocument(1)
    ei_probdef('millimeters', 'planar', 1e-10, depth, 30)
    ei_addmaterial('Диэлектрик учебного образца', epsx, epsy, 0)
    ei_addconductorprop('Питающий электрод', voltage, 0, 1)
    ei_addconductorprop('Заземлённый электрод', 0, 0, 1)
    polygon(((0, 0), (width, 0), (width, height), (0, height)))
    terminal(0, height / 2, 'Питающий электрод', mesh)
    terminal(width, height / 2, 'Заземлённый электрод', mesh)
    label(width / 2, height / 2, 'Диэлектрик учебного образца', mesh)
    model = directory + '/' + filename + '.fee'
    solve(model)
    return model


MODE='pair'

if MODE=='rogowski':
    voltage=200;gap=10;umin=-5;umax=1.5
    def mapped(u,v):
        scale=gap/math.pi
        return (scale*(u+math.exp(u)*math.cos(v)),scale*(v+math.exp(u)*math.sin(v)))
    for level,mesh in enumerate((.8,.4,.2)):
        newdocument(1);ei_probdef('millimeters','planar',1e-11,1000,30)
        ei_addmaterial('Воздух между профилями',1,1,0)
        ei_addconductorprop('Верхний профиль Роговского',voltage/2,0,1)
        ei_addconductorprop('Нижний профиль Роговского',-voltage/2,0,1)
        count=40*(2**level)
        sides=[([mapped(umin+(umax-umin)*j/count,math.pi/2) for j in range(count+1)],'Верхний профиль Роговского'),([mapped(umax,math.pi/2-math.pi*j/count) for j in range(count+1)],None),([mapped(umax-(umax-umin)*j/count,-math.pi/2) for j in range(count+1)],'Нижний профиль Роговского'),([mapped(umin,-math.pi/2+math.pi*j/count) for j in range(count+1)],None)]
        for points,name in sides:
            for x,y in points:ei_addnode(x,y)
            for j,(x,y) in enumerate(points[:-1]):
                x2,y2=points[j+1];ei_addsegment(x,y,x2,y2)
                if name:terminal((x+x2)/2,(y+y2)/2,name,mesh)
        label(0,0,'Воздух между профилями',mesh)
        model=directory+'/rogowski-%02d.fee'%(level+1);solve(model)
        exact=EPS0*(umax-umin)/math.pi
        row=summary(mesh,model,'Верхний профиль Роговского','Нижний профиль Роговского',exact)
        row['level']=level+1
        check('Конформная аналитическая ёмкость '+os.path.basename(model),row['analytic_C_relative_error'],.004)
        eo_smooth('on')
        samples=[]
        for u in (-3,0,1):
            x,y=mapped(u,0);value=field(x,y);theory=voltage/(gap*.001)/(1+math.exp(u));samples.append([u,x,value,theory])
            check('Поле конформного профиля u='+str(u)+' '+os.path.basename(model),abs(value/theory-1),.025/(2**level))
        row['midplane_samples_u_mm_V_m']=samples
        eo_close();ei_close()
    check('Сходимость ёмкости Роговского',abs(records[-1]['capacitance_F']/records[-2]['capacitance_F']-1),.002)

elif MODE=='pair':
    vacuum=[]
    for level,mesh in enumerate((.4,.2,.1)):
        for eps in (1,4.2):
            matrix=[[0,0],[0,0]];rows=[]
            for excitation in (0,1):
                newdocument(1);ei_probdef('millimeters','planar',1e-11,1000,30)
                ei_addmaterial('Воздух над парой',1,1,0);ei_addmaterial('Подложка пары',eps,eps,0)
                names=('Первая сигнальная полоса','Вторая сигнальная полоса')
                for i,name in enumerate(names):ei_addconductorprop(name,1 if i==excitation else 0,0,1)
                ei_addconductorprop('Опорный проводник и экран',0,0,1)
                polygon(((-10,-1.6),(10,-1.6),(10,0),(10,10),(-10,10),(-10,0)),'Опорный проводник и экран',mesh)
                xs=(-10,-1.25,-.25,.25,1.25,10)
                for x in xs:ei_addnode(x,0)
                for j in range(len(xs)-1):ei_addsegment(xs[j],0,xs[j+1],0)
                terminal(-.75,0,names[0],mesh/4);terminal(.75,0,names[1],mesh/4)
                # Сгущение обоих диэлектриков около краёв идеальных полос.
                for x in (-3,3):
                    for y in (-1.6,0,2):ei_addnode(x,y)
                    ei_addsegment(x,-1.6,x,0);ei_addsegment(x,0,x,2)
                ei_addsegment(-3,2,3,2)
                label(0,-.8,'Подложка пары',mesh/10);label(0,1,'Воздух над парой',mesh/10)
                label(-6,-.8,'Подложка пары',mesh/2);label(6,-.8,'Подложка пары',mesh/2);label(0,5,'Воздух над парой',mesh)
                model=directory+'/pair-eps%s-e%d-%02d.fee'%(str(eps).replace('.','p'),excitation+1,level+1);solve(model)
                charges=[real(eo_getconductorproperties(n)[1]) for n in names]
                for k in range(2):matrix[k][excitation]=charges[k]
                qground=real(eo_getconductorproperties('Опорный проводник и экран')[1])
                eo_groupselectblock(1);energy=real(eo_blockintegral(0));eo_clearblock()
                row={'model':os.path.basename(model),'level':level+1,'mesh_mm':mesh,'eps_r':eps,'excitation':excitation+1,'nodes':eo_numnodes(),'elements':eo_numelements(),'energy_J':energy,'first_charge_C':charges[0],'second_charge_C':charges[1]}
                check('Баланс трёх зарядов '+os.path.basename(model),abs(sum(charges)+qground)/max(abs(q) for q in charges),1e-6)
                check('Энергия и заряд матричного опыта '+os.path.basename(model),abs(2*energy/charges[excitation]-1),1e-5)
                records.append(row);rows.append(row);eo_close();ei_close()
            cd=(matrix[0][0]+matrix[1][1]-matrix[0][1]-matrix[1][0])/4
            cm=sum(sum(r) for r in matrix)
            check('Взаимность ёмкостей ε='+str(eps)+' сетка'+str(level),abs(matrix[0][1]/matrix[1][0]-1),1e-5)
            check('Пассивная ёмкостная матрица ε='+str(eps)+' сетка'+str(level),max(0,matrix[0][1]**2/(matrix[0][0]*matrix[1][1])-1),0)
            for row in rows:row['differential_C_F_per_m']=cd;row['common_C_F_per_m']=cm
            if eps==1:vacuum.append(cd)
            else:
                effective=cd/vacuum[-1];zd=1/(299792458*math.sqrt(cd*vacuum[-1]))
                for row in rows:row['effective_eps_differential']=effective;row['differential_Z_ohm']=zd
                check('Эффективная проницаемость между1иεr',max(0,1-effective,effective-eps),0)
    for eps in (1,4.2):
        selected=[r for r in records if r['eps_r']==eps and r['excitation']==1]
        check('Сходимость дифференциальной ёмкости ε='+str(eps),abs(selected[-1]['differential_C_F_per_m']/selected[-2]['differential_C_F_per_m']-1),.006)

elif MODE=='ring':
    major=20;outer=1000
    for radius in (.25,.5,1):
      for level,mesh in enumerate((1,.5,.25)):
        newdocument(1);ei_probdef('millimeters','axi',1e-11,1,30)
        ei_addmaterial('Воздух вокруг кольца',1,1,0)
        ei_addconductorprop('Проволочное кольцо',1,0,1);ei_addconductorprop('Удалённая сферическая граница',0,0,1)
        circle_xy(major,0,radius,'Проволочное кольцо',4/(2**level));label(major,0,'<No Mesh>',radius/5,0)
        ei_addnode(0,-outer);ei_addnode(0,outer);ei_addarc(0,-outer,0,outer,180,4/(2**level))
        ei_selectarcsegment(outer,0);ei_setarcsegmentprop(4/(2**level),'<None>',0,0,'Удалённая сферическая граница');ei_clearselected()
        ei_addsegment(0,-outer,0,outer)
        polygon(((0,-25),(50,-25),(50,25),(0,25)))
        label(10,0,'Воздух вокруг кольца',mesh/2);label(200,0,'Воздух вокруг кольца',40/(2**level))
        model=directory+'/ring-a%s-%02d.fee'%(str(radius).replace('.','p'),level+1);solve(model)
        # Ориентир тонкого кольца с ведущей монопольной поправкой конечной сферы.
        isolated=4*math.pi**2*EPS0*(major*.001)/math.log(8*major/radius)
        exact=1/(1/isolated-1/(4*math.pi*EPS0*outer*.001))
        row=summary(mesh,model,'Проволочное кольцо','Удалённая сферическая граница',exact)
        row.update({'level':level+1,'wire_radius_mm':radius,'slender_isolated_C_F':isolated,'outer_radius_mm':outer})
        check('Тонкое кольцо и поправка удалённой сферы '+os.path.basename(model),row['analytic_C_relative_error'],.015)
        eo_close();ei_close()
    for radius in (.25,.5,1):
        rr=[r for r in records if r['wire_radius_mm']==radius]
        check('Сходимость ёмкости кольца a='+str(radius),abs(rr[-1]['capacitance_F']/rr[-2]['capacitance_F']-1),.005)

elif MODE=='coupling':
    frequency=1000;eps=3.5;sigma=2e-8
    names=('Жила A — источник','Жила B — приёмник','Жила C — защитная')
    def current_circle(x,y,r,name,angle):
        ci_addnode(x-r,y);ci_addnode(x+r,y);ci_addarc(x-r,y,x+r,y,180,angle);ci_addarc(x+r,y,x-r,y,180,angle)
        ci_selectarcsegment(x,y+r);ci_selectarcsegment(x,y-r);ci_setarcsegmentprop(angle,'<None>',0,0,name);ci_clearselected()
    def current_label(x,y,name,mesh,group=1):
        ci_addblocklabel(x,y);ci_selectlabel(x,y);ci_setblockprop(name,0,mesh,group);ci_clearselected()
    def cable(voltages,mesh,level,stem):
        newdocument(3);ci_probdef('millimeters','planar',frequency,1e-11,1000,30)
        ci_addmaterial('Изоляция с диэлектрическими потерями',sigma,sigma,eps,eps,0,0)
        ci_addconductorprop('Заземлённая оболочка',0,0,1)
        current_circle(0,0,10,'Заземлённая оболочка',3/(2**level))
        for k,(x,y) in enumerate(((-4,0),(0,4),(4,0))):
            ci_addconductorprop(names[k],voltages[k],0,1);current_circle(x,y,1.5,names[k],3/(2**level));current_label(x,y,'<No Mesh>',mesh,0)
        current_label(0,-4,'Изоляция с диэлектрическими потерями',mesh)
        model=directory+'/'+stem+'.fec';ci_saveas(model);ci_close();opendocument(model);ci_analyze(1);ci_loadsolution();co_smooth('off')
        currents=[complex(co_getconductorproperties(n)[1]) for n in names]
        currents.append(complex(co_getconductorproperties('Заземлённая оболочка')[1]))
        co_groupselectblock(1);power=real(co_blockintegral(0));co_clearblock()
        inputpower=.5*sum((complex(voltages[k])*currents[k].conjugate()).real for k in range(3))
        row={'model':os.path.basename(model),'level':level+1,'mesh_mm':mesh,'nodes':co_numnodes(),'elements':co_numelements(),'dielectric_loss_W':power,'terminal_power_W':inputpower}
        check('Баланс токов четырёх электродов '+stem,abs(sum(currents))/max(abs(i) for i in currents),1e-6)
        check('Баланс диэлектрической мощности '+stem,abs(power-inputpower)/max(power,1e-20),1e-5)
        co_close();ci_close()
        return row,currents
    for level,mesh in enumerate((.5,.25,.125)):
        matrix=[[0j]*3 for _ in range(3)]
        for exc in range(3):
            vv=[0,0,0];vv[exc]=1
            row,currents=cable(vv,mesh,level,'Y%d-%02d'%(exc+1,level+1));row['excitation']=exc+1
            for k in range(3):matrix[k][exc]=currents[k]
            records.append(row)
        for i in range(3):
          for j in range(3):check('Взаимность комплексной матрицы Y%d%d сетка%d'%(i,j,level+1),abs(matrix[i][j]-matrix[j][i])/max(abs(matrix[i][i]),1e-20),1e-5)
        for load in (1e6,1e7,1e8):
            induced=-matrix[1][0]*100/(matrix[1][1]+1/load)
            row,currents=cable([100,induced,0],mesh,level,'load%s-%02d'%(str(int(load)),level+1))
            row.update({'load_ohm':load,'receiver_voltage_real_V':induced.real,'receiver_voltage_imag_V':induced.imag,'receiver_voltage_amplitude_V':abs(induced),'load_power_W':abs(induced)**2/(2*load),'source_power_W':.5*(100*currents[0].conjugate()).real})
            check('Узел нагруженной жилы '+row['model'],abs(currents[1]+induced/load)/max(abs(induced/load),1e-20),1e-5)
            check('Источник питает потери изоляции и нагрузку '+row['model'],abs(row['source_power_W']-row['dielectric_loss_W']-row['load_power_W'])/row['source_power_W'],1e-5)
            records.append(row)
    for load in (1e6,1e7,1e8):
        rr=[r for r in records if r.get('load_ohm')==load]
        check('Сходимость наведённого напряжения R='+str(load),abs(rr[-1]['receiver_voltage_amplitude_V']/rr[-2]['receiver_voltage_amplitude_V']-1),.005)

result={'runtime':'FEMM 5; embedded MicroPython','records':records,'checks':checks,'passed':all(c['passed'] for c in checks),'scope':'Самостоятельная учебная электротехническая модель; численные проверки относятся к точным сохранённым файлам.'}
with open(directory+'/result.json','w') as f:json.dump(result,f)
print(json.dumps(result));assert result['passed']