"""Сопоставление катушки курса с теорией конечного соленоида.

Среда: встроенный Python FEMM 5. Входная учебная копия Катушка-теория.fem
лежит рядом с программой. Каждый запуск сохраняет отдельный calculated-NN.
Радиусы 20/30 мм, длина 200 мм, N=100, I=1 А; свойства русские.
Внешняя полуокружность R=300 мм отличается от прямоугольной дальней
границы браузерного решателя. Это отдельный источник погрешности энергии.
"""
import os, sys, math, json

source = globals().get('__file__', '') or sys.argv[0]
base = os.path.dirname(os.path.abspath(source))
index = 1
while os.path.exists(base + '/calculated-%02d' % index): index += 1
out = base + '/calculated-%02d' % index
os.mkdir(out)

opendocument(base + '/Катушка-теория.fem')
# Expand the exterior boundary independently of winding geometry and sources.
mi_selectnode(0, 300)
mi_movetranslate(0, 600, 0)
mi_clearselected()
mi_selectnode(0, -300)
mi_movetranslate(0, -600, 0)
mi_clearselected()
mi_selectlabel(210, 0)
mi_setblockprop('Воздух', 0, 10, '<None>', 0, 0, 0)
mi_clearselected()
mi_selectlabel(25, 0)
mi_setblockprop('Обмоточная область', 0, 0.5, 'Катушка', 0, 1, 100)
mi_clearselected()
for x,y in ((0,-180),(15,-180),(15,180),(0,180)):
    mi_addnode(x,y)
for a,b in (((0,-180),(15,-180)),((15,-180),(15,180)),((15,180),(0,180))):
    mi_addsegment(a[0],a[1],b[0],b[1])
mi_addblocklabel(7,0)
mi_selectlabel(7,0)
mi_setblockprop('Воздух',0,.75,'<None>',0,0,0)
mi_clearselected()
path = out + '/Катушка-проверка.fem'
mi_saveas(path)
mi_close()
opendocument(path)
mi_analyze(1)
mi_loadsolution()
current, voltage, linkage = mo_getcircuitproperties('Катушка')
current = complex(current).real
linkage = complex(linkage).real
mo_clearblock()
mo_groupselectblock()
energy = complex(mo_blockintegral(2)).real
mo_clearblock()

def analytic(z):
    # Uniform rectangular winding: average finite thin-solenoid field over radius.
    a, b, length, turns = 0.02, 0.03, 0.2, 100
    count = 100
    step = (b-a)/count
    total = 0
    for k in range(count+1):
        r = a+k*step
        v = (z+length/2)/math.sqrt((z+length/2)**2+r*r) - (z-length/2)/math.sqrt((z-length/2)**2+r*r)
        total += v*(1 if k==0 or k==count else 4 if k%2 else 2)
    return 4e-7*math.pi*turns*current/(2*length*(b-a))*total*step/3

rows=[]
for z in (0,50,100,150):
    values=mo_getpointvalues(.001,z)
    field=complex(values[2]).real
    exact=analytic(z/1000)
    error=abs(field-exact)/max(abs(exact),1e-15)
    rows.append({'z_mm':z,'Bz_T':field,'analytic_Bz_T':exact,'relativeError':error})
energy_error=abs(linkage*current-2*energy)/max(abs(linkage*current),1e-15)
result={'model':'Катушка-проверка.fem','solution':'Катушка-проверка.ans',
        'current_A':current,'linkage_Wb':linkage,'energy_J':energy,
        'L_linkage_H':linkage/current,'L_energy_H':2*energy/(current*current),
        'energyIdentityError':energy_error,'axisProfile':rows,
        'scope':'Air coil; semicircular external boundary R=900 mm; local axis air strip mesh 0.75 mm; winding mesh 0.5 mm; no leads',
        'passed':all(r['relativeError']<.01 for r in rows) and energy_error<.002}
mo_close()
mi_close()
with open(out+'/coil-result.json','w') as stream:json.dump(result,stream)
print(json.dumps(result))
assert result['passed']