"""Встроенный Python FEMM 5: поле провода, магнитная стрелка и аналогии полей.

Запустите в отдельной рабочей папке. Каждый повтор создаёт calculated-NN.
Геометрия в миллиметрах; интегральные результаты — СИ. Возбуждение — ток,
постоянная намагниченность либо напряжение на электродах; все условия заданы
здесь. Механический и RC-расчёты выполняются явно после статических полей.
"""
import os, sys, math, json

script=globals().get('__file__','') or (sys.argv[0] if sys.argv else '')
base=os.path.dirname(os.path.abspath(script)) if script else os.getcwd()
index=1
while 'calculated-%02d'%index in os.listdir(base):index+=1
directory=base+'/calculated-%02d'%index
os.mkdir(directory)
MU0=4e-7*math.pi
rows=[]
checks=[]
def real(value):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 save():
    with open(directory+'/results.json','w') as stream:json.dump({'records':rows,'checks':checks},stream)
def magnetic_material(name,mu=1,hc=0,j=0):
    mi_addmaterial(name,mu,mu,hc,j,0,0,0,1,0,0,0,0)
def circle(cx,cy,radius,angle=5):
    mi_addnode(cx-radius,cy);mi_addnode(cx+radius,cy)
    mi_addarc(cx-radius,cy,cx+radius,cy,180,angle)
    mi_addarc(cx+radius,cy,cx-radius,cy,180,angle)
def rect(x1,y1,x2,y2):mi_drawrectangle(x1,y1,x2,y2)
def label(x,y,material,mesh,group=0,circuit='<None>',turns=0,magangle=0):
    mi_addblocklabel(x,y);mi_selectlabel(x,y)
    mi_setblockprop(material,0,mesh,circuit,magangle,group,turns);mi_clearselected()
def magnetic_boundary(cx,radius):
    circle(cx,0,radius,5)
    mi_addboundprop('Внешняя граница A=0',0,0,0,0,0,0,0,0,0)
    for y in [-radius,radius]:
        mi_selectarcsegment(cx,y);mi_setarcsegmentprop(5,'Внешняя граница A=0',0,0);mi_clearselected()
def solve_magnetic(name):
    path=directory+'/'+name+'.fem'
    mi_saveas(path);mi_close();opendocument(path);mi_analyze(1);mi_loadsolution()
    return {'model':name+'.fem','solution':name+'.ans','nodes':mo_numnodes(),'elements':mo_numelements()}
def rotated(x,y,angle):
    a=angle*math.pi/180
    return x*math.cos(a)-y*math.sin(a),x*math.sin(a)+y*math.cos(a)
def rectangles_rotated(pieces,angle):
    points={};edges={}
    for piece in pieces:
        x1,y1,x2,y2=piece
        poly=[rotated(x1,y1,angle),rotated(x2,y1,angle),rotated(x2,y2,angle),rotated(x1,y2,angle)]
        for p in poly:
            key=(round(p[0],10),round(p[1],10))
            if key not in points:points[key]=p;mi_addnode(p[0],p[1])
        for i,p in enumerate(poly):
            q=poly[(i+1)%4]
            kp=(round(p[0],10),round(p[1],10));kq=(round(q[0],10),round(q[1],10))
            key=tuple(sorted((kp,kq)))
            if key not in edges:edges[key]=True;mi_addsegment(p[0],p[1],q[0],q[1])

# Эрстед: нечётная по внешнему току часть отделяет взаимодействие от
# собственного поля стрелки и от взаимодействия её поля с внешней границей.
def needle(kind,mu,thickness,angle,current,mesh,radius=200):
    newdocument(0);mi_probdef(0,'millimeters','planar',1e-9,2,30)
    magnetic_material('Воздух')
    magnetic_material('Сердцевина стрелки',mu)
    magnetic_material('Постоянная намагниченность',mu,100000)
    magnetic_material('Левый токовый слой',1,0,100000/(thickness*1e-3)/1e6)
    magnetic_material('Правый токовый слой',1,0,-100000/(thickness*1e-3)/1e6)
    magnetic_material('Внешний провод',1,0,current/(math.pi*0.5**2))
    magnetic_boundary(25,radius)
    circle(25,0,0.5,8);label(25,0,'Внешний провод',0.15,40)
    circle(0,0,15,2);label(0,13,'Воздух',mesh*8)
    if kind=='pm':
        rectangles_rotated([(-1,-10,1,10)],angle)
        label(0,0,'Постоянная намагниченность',mesh,20,magangle=90+angle)
    else:
        t=thickness
        rectangles_rotated([(-1-t/2,-10,-1+t/2,10),(-1+t/2,-10,1-t/2,10),(1-t/2,-10,1+t/2,10)],angle)
        for x,material in [(-1,'Левый токовый слой'),(0,'Сердцевина стрелки'),(1,'Правый токовый слой')]:
            px,py=rotated(x,0,angle);label(px,py,material,min(mesh,t/3),20)
    label(50,30,'Воздух',8)
    name='Стрелка-%s-mu%s-t%s-a%s-I%s-h%s'%(kind,str(mu).replace('.','p'),str(thickness).replace('.','p'),angle,current,str(mesh).replace('.','p'))
    row=solve_magnetic(name)
    mo_groupselectblock(20);row['torque_Nm']=real(mo_blockintegral(22));mo_clearblock()
    mo_groupselectblock(40);row['wire_force_y_N']=real(mo_blockintegral(12));mo_clearblock()
    row['reaction_torque_Nm']=-0.025*row['wire_force_y_N']
    mo_close();mi_close()
    row.update({'experiment':'needle','kind':kind,'mu_rec':mu,'layer_mm':thickness,'angle_deg':angle,'wire_current_A':current,'mesh_mm':mesh})
    rows.append(row);save();return row

def analytical_needle(angle,thickness=0,wire_current=100,segments=160):
    # Тонкие амперовские слои с K=Hc=100000 A/m, глубиной 2 мм.
    # Интегрирование по длине сохраняет полный ток Hc*L=2000 A каждого слоя.
    a=angle*math.pi/180;torque=0;cross_energy=0
    nwidth=1 if thickness==0 else 4
    for side,sign in [(-1,1),(1,-1)]:
        for ix in range(nwidth):
            x0=side*1e-3+(0 if thickness==0 else ((ix+0.5)/nwidth-0.5)*thickness*1e-3)
            for k in range(segments):
                y0=(-10+(k+0.5)*20/segments)*1e-3
                x=x0*math.cos(a)-y0*math.sin(a);y=x0*math.sin(a)+y0*math.cos(a)
                dx=x-0.025;r2=dx*dx+y*y
                bx=-MU0*wire_current*y/(2*math.pi*r2);by=MU0*wire_current*dx/(2*math.pi*r2)
                di=sign*100000*0.02/segments/nwidth
                fx=-di*0.002*by;fy=di*0.002*bx
                torque+=x*fy-y*fx
                potential=-MU0*wire_current/(4*math.pi)*math.log(r2)
                cross_energy+=di*0.002*potential
    return torque,cross_energy

needle_pairs=[]
for mu in [1,6.678]:
    for kind,t in [('pm',0.125),('coil',0.5),('coil',0.25),('coil',0.125)]:
        pair=[needle(kind,mu,t,30,sign*100,0.06) for sign in [1,-1]]
        cross=(pair[0]['torque_Nm']-pair[1]['torque_Nm'])/2
        needle_pairs.append({'kind':kind,'mu_rec':mu,'layer_mm':t,'cross_torque_Nm':cross,'models':[r['model'] for r in pair]})
for angle in [0,60]:
    pair=[needle('pm',1,0.125,angle,sign*100,0.06) for sign in [1,-1]]
    cross=(pair[0]['torque_Nm']-pair[1]['torque_Nm'])/2
    analytical=analytical_needle(angle)[0]
    needle_pairs.append({'kind':'pm','mu_rec':1,'angle_deg':angle,'cross_torque_Nm':cross,'analytical_torque_Nm':analytical})
    check('Момент Эрстеда при угле %d'%angle,abs(cross-analytical)/max(abs(analytical),abs(analytical_needle(30)[0])),0.03)
torque_mesh=[]
for h in [0.035,0.02,0.01]:
    fine=[needle('pm',1,0.125,30,sign*100,h) for sign in [1,-1]]
    torque_mesh.append({'mesh_mm':h,'cross_torque_Nm':(fine[0]['torque_Nm']-fine[1]['torque_Nm'])/2})
fine_torque=torque_mesh[-1]['cross_torque_Nm']
pm30=needle_pairs[0]['cross_torque_Nm'];exact30=analytical_needle(30)[0]
check('Момент Эрстеда, тонкие слои/аналитика',abs(fine_torque-exact30)/abs(exact30),0.03)
check('Уточнение сетки момента Эрстеда',abs(fine_torque-torque_mesh[-2]['cross_torque_Nm'])/abs(fine_torque),0.02)

# Фарадей: пассивная N=100 рамка возле длинного круглого провода.
# Две стороны рамки, длина 100 мм; внешнее сопротивление и ответный ток отсутствуют.
faraday=[]
for mesh in [0.18,0.09,0.045]:
    newdocument(0);mi_probdef(0,'millimeters','planar',1e-10,100,30)
    magnetic_material('Воздух');magnetic_material('Источник тока',1,0,10/(math.pi*0.5**2))
    mi_addcircprop('Измерительная рамка',0,1)
    magnetic_boundary(0,100);circle(0,0,0.5,6);label(0,0,'Источник тока',0.1)
    for x,turns in [(10,100),(20,-100)]:
        rect(x-0.1,-0.1,x+0.1,0.1);label(x,0,'Воздух',mesh,30,'Измерительная рамка',turns)
    label(40,20,'Воздух',mesh*8)
    name='Рамка-Фарадея-h'+str(mesh).replace('.','p')
    row=solve_magnetic(name)
    props=mo_getcircuitproperties('Измерительная рамка')
    row.update({'experiment':'faraday','mesh_mm':mesh,'flux_linkage_Wb':real(props[2]),'mutual_H':real(props[2])/10})
    field_values=mo_getpointvalues(15,0)
    row['B_at_15mm_T']=math.sqrt(real(field_values[1])**2+real(field_values[2])**2)
    mo_close();mi_close();rows.append(row);faraday.append(row);save()
mutual_exact=MU0*100*0.1/(2*math.pi)*math.log(20/10)
check('Рамка: M по полю и аналитике',abs(faraday[-1]['mutual_H']-mutual_exact)/mutual_exact,0.005)
check('Рамка: сеточное уточнение M',abs(faraday[-1]['mutual_H']-faraday[0]['mutual_H'])/mutual_exact,0.005)
emf_peak=2*math.pi*50*10*faraday[-1]['mutual_H']
derivative=[]
for count in [48,96,192]:
    dt=1/(50*count)
    psi=[10*faraday[-1]['mutual_H']*math.sin(2*math.pi*k/count) for k in range(count)]
    emf=[-(psi[(k+1)%count]-psi[(k-1)%count])/(2*dt) for k in range(count)]
    rms=math.sqrt(sum(v*v for v in emf)/count)
    derivative.append({'samples':count,'emf_rms_V':rms,'relative_error':abs(rms-emf_peak/math.sqrt(2))/(emf_peak/math.sqrt(2))})
check('Рамка: производная потокосцепления',derivative[-1]['relative_error'],0.001)

# Прямая и обратная аналогии: одна L-образная область и сопряжённые границы.
POLYGON=[(0,0),(40,0),(40,10),(20,10),(20,30),(0,30)]
EDGES=[((20,0),'Низ'),((40,5),'Право'),((30,10),'Верх'),((20,20),'Верх'),((10,30),'Верх'),((0,15),'Лево')]
def current_case(dual,mesh):
    newdocument(3);ci_probdef('millimeters','planar',0,1e-10,100,30)
    ci_addmaterial('Аналоговая проводящая среда',2,2,1,1,0,0)
    ci_addconductorprop('Входной электрод',1,0,1);ci_addconductorprop('Выходной электрод',0,0,1)
    for x,y in POLYGON:ci_addnode(x,y)
    for i,p in enumerate(POLYGON):
        q=POLYGON[(i+1)%len(POLYGON)];ci_addsegment(p[0],p[1],q[0],q[1])
    for (x,y),edge in EDGES:
        name=('Входной электрод' if edge==('Низ' if dual else 'Лево') else 'Выходной электрод' if edge==('Верх' if dual else 'Право') else None)
        if name:
            ci_selectsegment(x,y);ci_setsegmentprop('<None>',mesh,0,0,0,name);ci_clearselected()
    ci_addblocklabel(10,5);ci_selectlabel(10,5);ci_setblockprop('Аналоговая проводящая среда',0,mesh,1);ci_clearselected()
    name='Аналог-'+('обратный' if dual else 'прямой')+'-h'+str(mesh).replace('.','p')
    ci_saveas(directory+'/'+name+'.fec');ci_close();opendocument(directory+'/'+name+'.fec');ci_analyze(1);ci_loadsolution()
    vin,iin=co_getconductorproperties('Входной электрод');vout,iout=co_getconductorproperties('Выходной электрод')
    conductance=abs(real(iin)/(real(vin)-real(vout)))
    samples=[{'x_mm':x,'y_mm':y,'values':[real(v) for v in co_getpointvalues(x,y)]} for x,y in [(5,5),(15,5),(25,5),(10,15),(10,25)]]
    corner=[]
    if not dual:
        for distance in [0.5,1,2,4]:
            v=co_getpointvalues(20-distance/math.sqrt(2),10-distance/math.sqrt(2))
            corner.append({'distance_mm':distance,'current_density_A_m2':math.sqrt(real(v[1])**2+real(v[2])**2)})
    row={'experiment':'analogy-current','dual':dual,'mesh_mm':mesh,'conductance_S':conductance,'lambda':conductance/(2*0.1),'current_balance_A':abs(real(iin)+real(iout)),'samples':samples,'corner_samples':corner,'model':name+'.fec','solution':name+'.anc','nodes':co_numnodes(),'elements':co_numelements()}
    co_close();ci_close();rows.append(row);save();return row
def magnetic_analogy(mesh):
    newdocument(0);mi_probdef(0,'millimeters','planar',1e-10,100,30);magnetic_material('Немагнитный зазор')
    for x,y in POLYGON:mi_addnode(x,y)
    for i,p in enumerate(POLYGON):
        q=POLYGON[(i+1)%len(POLYGON)];mi_addsegment(p[0],p[1],q[0],q[1])
    for name,value in [('Нижняя линия потока',0),('Верхняя линия потока',1e-4)]:mi_addboundprop(name,value,0,0,0,0,0,0,0,0)
    for (x,y),edge in EDGES:
        if edge in ['Низ','Верх']:
            mi_selectsegment(x,y);mi_setsegmentprop('Нижняя линия потока' if edge=='Низ' else 'Верхняя линия потока',mesh,0,0,0);mi_clearselected()
    label(10,5,'Немагнитный зазор',mesh)
    name='Магнитный-зазор-аналог-h'+str(mesh).replace('.','p');row=solve_magnetic(name)
    mo_addcontour(0,5);mo_addcontour(40,5);integral=mo_lineintegral(1);mo_clearcontour()
    mmf=abs(real(integral[0] if isinstance(integral,(list,tuple)) else integral))
    flux=0.1*1e-4;permeance=flux/mmf
    row.update({'experiment':'analogy-magnetic','mesh_mm':mesh,'mmf_A':mmf,'flux_Wb':flux,'permeance_H':permeance,'lambda':permeance/(MU0*0.1)})
    mo_close();mi_close();rows.append(row);save();return row
analog=[]
for mesh in [1,0.5,0.25]:
    direct=current_case(False,mesh);dual=current_case(True,mesh);mag=magnetic_analogy(mesh)
    analog.append({'mesh_mm':mesh,'direct_lambda':direct['lambda'],'dual_lambda':dual['lambda'],'duality_product':direct['lambda']*dual['lambda'],'magnetic_lambda':mag['lambda'],'magnetic_permeance_H':mag['permeance_H']})
check('Лаплас: произведение сопряжённых проводимостей',abs(analog[-1]['duality_product']-1),0.015)
check('Прямая электрическая и магнитная аналогии',abs(analog[-1]['magnetic_lambda']-analog[-1]['direct_lambda'])/analog[-1]['direct_lambda'],0.02)
check('Аналогия: сеточная устойчивость проводимости',abs(analog[-1]['direct_lambda']-analog[-2]['direct_lambda'])/analog[-1]['direct_lambda'],0.015)
corner=direct['corner_samples']
corner_exponent=math.log(corner[-1]['current_density_A_m2']/corner[1]['current_density_A_m2'])/math.log(corner[-1]['distance_mm']/corner[1]['distance_mm'])
check('Угол 270 градусов: показатель Шварца–Кристоффеля',abs(corner_exponent+1/3),0.08)

def resistor_network(h):
    nx=int(round(40/h));ny=int(round(30/h));cells=[];lookup={}
    for j in range(ny):
        for i in range(nx):
            if (i+0.5)*h<20 or (j+0.5)*h<10:
                lookup[(i,j)]=len(cells);cells.append((i,j))
    neighbors=[];fixed=[];diagonal=[];voltage=[]
    for i,j in cells:
        ns=[];rhs=0;diag=0
        for di,dj in [(1,0),(-1,0),(0,1),(0,-1)]:
            p=(i+di,j+dj)
            if p in lookup:ns.append(lookup[p]);diag+=1
            elif di==-1 and i==0:diag+=2;rhs+=2
            elif di==1 and i==nx-1:diag+=2
        neighbors.append(ns);fixed.append(rhs);diagonal.append(diag);voltage.append(1-(i+0.5)/nx)
    residual=1
    for sweep in range(6000):
        residual=0
        for k in range(len(cells)):
            candidate=(fixed[k]+sum(voltage[n] for n in neighbors[k]))/diagonal[k]
            delta=1.7*(candidate-voltage[k]);voltage[k]+=delta;residual=max(residual,abs(delta))
        if residual<1e-11:break
    factor=sum(2*(1-voltage[k]) for k,(i,j) in enumerate(cells) if i==0)
    return {'cell_mm':h,'nodes':len(cells),'iterations':sweep+1,'residual_V':residual,'lambda':factor,'relative_to_FE':abs(factor-analog[-1]['direct_lambda'])/analog[-1]['direct_lambda']}
networks=[resistor_network(h) for h in [5,2.5,1.25]]
check('Резистивная сеть: итерационная невязка',networks[-1]['residual_V'],1e-9)
check('Резистивная сеть относительно непрерывного FE',networks[-1]['relative_to_FE'],0.035)

# RC-аналог магнитной диффузии: одинаковое уравнение с коэффициентом μσ.
# Физическая схема не изготовляется; решена её комплексная узловая система.
def rc_diffusion(count,frequency):
    conductivity=58e6;delta=math.sqrt(2/(2*math.pi*frequency*MU0*conductivity));length=5*delta;dx=length/count
    r_sheet=1.0;c_area=MU0*conductivity/r_sheet
    s=1j*2*math.pi*frequency*r_sheet*c_area*dx*dx
    a=[-1.0]*(count-1);b=[2+s]*(count-1);c=[-1.0]*(count-1);d=[0j]*(count-1);d[0]=1
    for k in range(1,count-1):
        factor=a[k]/b[k-1];b[k]-=factor*c[k-1];d[k]-=factor*d[k-1]
    v=[0j]*(count+1);v[0]=1;v[-2]=d[-1]/b[-1]
    for k in range(count-3,-1,-1):v[k+1]=(d[k]-c[k]*v[k+2])/b[k]
    q=complex(1,1)/delta;x=length*0.4
    def sinh(z):return (complex(math.e)**z-complex(math.e)**(-z))/2
    exact=sinh(q*(length-x))/sinh(q*length)
    sample=v[int(round(0.4*count))]
    return {'frequency_Hz':frequency,'cells':count,'skin_depth_m':delta,'R_sheet_ohm':r_sheet,'C_per_area_F_m2':c_area,'x_over_delta':2,'transfer_re':sample.real,'transfer_im':sample.imag,'exact_re':exact.real,'exact_im':exact.imag,'relative_error':abs(sample-exact)/abs(exact)}
rc=[rc_diffusion(n,f) for f in [50,500] for n in [20,40,80]]
for r in rc:
    if r['cells']==80:check('RC-диффузия %d Гц'%r['frequency_Hz'],r['relative_error'],0.003)

# Нелинейный обратный аналог: исходные последовательные падения МДС
# превращаются в параллельные токи при общем аналоге потока U.
nonlinear=[]
area=0.02*0.1;steel_length=0.2;gap=0.001;kphi=1000;kmmf=0.01
def h_of_b(b):return 120*b+900*b**3
for mmf in [100,300,600]:
    lo=0;hi=3
    for iteration in range(80):
        mid=(lo+hi)/2
        if steel_length*h_of_b(mid)+gap*mid/MU0>mmf:hi=mid
        else:lo=mid
    direct=(lo+hi)/2;u_lo=0;u_hi=3*area*kphi
    for iteration in range(80):
        u=(u_lo+u_hi)/2;b=u/(kphi*area)
        current=kmmf*(steel_length*h_of_b(b)+gap*b/MU0)
        if current>kmmf*mmf:u_hi=u
        else:u_lo=u
    voltage=(u_lo+u_hi)/2;recovered=voltage/(kphi*area)
    nonlinear.append({'mmf_A':mmf,'B_T':direct,'flux_Wb':direct*area,'analog_voltage_V':voltage,'analog_source_A':kmmf*mmf,'recovered_B_T':recovered,'residual_A':abs(steel_length*h_of_b(direct)+gap*direct/MU0-mmf),'scale_voltage_per_Wb':kphi,'scale_current_per_Aturn':kmmf})
    check('Нелинейная дуальная схема %d А'%mmf,abs(recovered-direct),1e-10)

# Механика стрелки: RK4 по аналитическому моменту, проверенному точками FE.
mass=7500*0.002*0.02*0.002;inertia=mass*(0.02**2+0.002**2)/12
table=[analytical_needle(k*2,segments=80) for k in range(180)]
def interpolate(angle,component):
    x=(angle*180/math.pi)%360/2;k=int(x);u=x-k
    return table[k][component]*(1-u)+table[(k+1)%180][component]*u
def mechanics(dt,damping):
    theta=math.pi/4;omega=0;loss=0;initial=-interpolate(theta,1);max_balance=0
    def acceleration(a,w):return (interpolate(a,0)-damping*w)/inertia
    snapshots=[]
    for k in range(int(round(1/dt))):
        a1=omega;b1=acceleration(theta,omega)
        a2=omega+dt*b1/2;b2=acceleration(theta+dt*a1/2,omega+dt*b1/2)
        a3=omega+dt*b2/2;b3=acceleration(theta+dt*a2/2,omega+dt*b2/2)
        a4=omega+dt*b3;b4=acceleration(theta+dt*a3,omega+dt*b3)
        old=omega;theta+=dt*(a1+2*a2+2*a3+a4)/6;omega+=dt*(b1+2*b2+2*b3+b4)/6
        loss+=damping*(old*old+omega*omega)*dt/2
        balance=0.5*inertia*omega*omega-interpolate(theta,1)+loss-initial
        max_balance=max(max_balance,abs(balance))
        if (k+1)%int(round(0.05/dt))==0:snapshots.append({'time_s':(k+1)*dt,'angle_deg':theta*180/math.pi,'omega_rad_s':omega})
    return {'dt_s':dt,'damping_Nm_s':damping,'mass_kg':mass,'inertia_kg_m2':inertia,'end_angle_deg':theta*180/math.pi,'end_speed_rad_s':omega,'maximum_energy_balance_error_J':max_balance,'relative_energy_error':max_balance/max(abs(initial),1e-12),'samples':snapshots}
dynamic=[mechanics(dt,damping) for damping in [0,1e-7] for dt in [0.001,0.0005]]
for r in dynamic:check('Механика: баланс энергии dt=%g c=%g'%(r['dt_s'],r['damping_Nm_s']),r['relative_energy_error'],0.004)

result={'records':rows,'checks':checks,'needle_pairs':needle_pairs,'torque_mesh':torque_mesh,'analytical_torque_30_Nm':exact30,'fine_torque_30_Nm':fine_torque,'faraday_mutual_exact_H':mutual_exact,'faraday_emf_peak_V':emf_peak,'faraday_derivative':derivative,'analogies':analog,'corner_current_exponent':corner_exponent,'resistor_networks':networks,'rc_diffusion':rc,'nonlinear_analog':nonlinear,'mechanics':dynamic,'passed':all(c['passed'] for c in checks)}
with open(directory+'/results.json','w') as stream:json.dump(result,stream)
print('FOUNDATIONS_COMPLETE',result['passed'],len(rows),len(checks),directory)