"""Четыре режима двух воздушных катушек: встроенный Python FEMM 5.

Рядом должен находиться multicoil-template.fem — отдельная русская
учебная копия точного готового примера I04.2. Исходник не меняется.
В модели НЕТ магнитопровода: использованы лишь Воздух и Обмоточная
область с mu=1,sigma=0. Обе катушки имеют положительные N=100/50.
Отрицательный I2 означает противодействующую МДС второй катушки.
"""
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()
TEMPLATE = globals().get('BUL_MULTICOIL_TEMPLATE',BASE+'/multicoil-template.fem')
OUTPUT = globals().get('BUL_MULTICOIL_OUTPUT_ROOT',BASE)
n=1
while os.path.exists(OUTPUT+'/calculated-%02d'%n):
    n+=1
OUT=OUTPUT+'/calculated-%02d'%n
os.mkdir(OUT)
NAMES=['Первая катушка','Вторая катушка']
RECORDS=[]
CHECKS=[]
VOLUMES=[math.pi*(.030**2-.020**2)*.200,math.pi*(.040**2-.035**2)*.200]


def real(value):
    return float(complex(value).real)


def check(name,passed,observed,limit,scope):
    CHECKS.append({'id':name,'passed':bool(passed),'observed':observed,'limit':limit,'scope':scope})


def relative(a,b):
    return abs(a-b)/max(abs(b),1e-20)


for tag,currents in [('только-первая',[1,0]),('только-вторая',[0,1]),('согласные-токи',[1,1]),('встречные-токи',[2,-.5])]:
    opendocument(TEMPLATE)
    for name,current in zip(NAMES,currents):
        mi_modifycircprop(name,1,current)
    mi_setcomment('Русская учебная копия примера I04.2: две воздушные соосные катушки. '+
        'Магнитопровода нет, mu=1,sigma=0. N1=100,N2=50, оба положительны. '+
        'I1=%gА,I2=%gА. '%(currents[0],currents[1])+
        'Группы: 1 — первая катушка; 2 — вторая катушка; 3 — воздух. '+
        'Осесимметрия: 2pi уже включено в объёмы, потокосцепления и энергию.')
    stem='%02d-Две-воздушные-катушки-%s'%(len(RECORDS)+1,tag)
    path=OUT+'/'+stem+'.fem'
    mi_saveas(path)
    mi_close()
    opendocument(path)
    mi_analyze(1)
    mi_loadsolution()
    circuit=[mo_getcircuitproperties(name)for name in NAMES]
    psi=[real(v[2])for v in circuit]
    mo_clearblock()
    mo_groupselectblock()
    energy=real(mo_blockintegral(2))
    coenergy=real(mo_blockintegral(17))
    mo_clearblock()
    volumes=[]
    for group in [1,2]:
        mo_groupselectblock(group)
        volumes.append(real(mo_blockintegral(10)))
        mo_clearblock()
    row={'tag':tag,'model':stem+'.fem','solution':stem+'.ans','currents_A':currents,
        'linkages_Wb':psi,'energy_J':energy,'coenergy_J':coenergy,
        'half_Ipsi_J':sum(currents[k]*psi[k]for k in range(2))/2,
        'coil_volumes_m3':volumes,'coil_volumes_exact_m3':VOLUMES,
        'elements':mo_numelements(),'nodes':mo_numnodes(),'saved_reopened_solved':True}
    mo_close()
    mi_close()
    opendocument(path)
    mi_loadsolution()
    reopened=[real(mo_getcircuitproperties(name)[2])for name in NAMES]
    row['solution_reopened_matches']=all(abs(reopened[k]-psi[k])<1e-12 for k in range(2))
    mo_close()
    mi_close()
    RECORDS.append(row)
    check(tag+'-roundtrip',row['solution_reopened_matches'],reopened,psi,'Русские цепи и сохранённые решения открыты повторно.')
    check(tag+'-energy',relative(energy,row['half_Ipsi_J'])<.003,relative(energy,row['half_Ipsi_J']),.003,'Линейная задача: W=I^T Psi/2, СИ.')
    check(tag+'-positive-energy',energy>0,energy,0,'Положительная энергия для ненулевого вектора токов.')
    for k in range(2):
        error=relative(volumes[k],VOLUMES[k])
        check(tag+'-volume-'+str(k+1),error<1e-8,error,1e-8,'Полный объём кольцевой обмотки; дополнительное2pi не вводится.')
    print('BUL_MULTICOIL_SOLVED',tag,energy,psi)

L=[[RECORDS[0]['linkages_Wb'][0],RECORDS[1]['linkages_Wb'][0]],
   [RECORDS[0]['linkages_Wb'][1],RECORDS[1]['linkages_Wb'][1]]]
M=(L[0][1]+L[1][0])/2
determinant=L[0][0]*L[1][1]-M*M
reciprocity=abs(L[0][1]-L[1][0])/max(abs(M),1e-20)
check('reciprocity',reciprocity<.001,reciprocity,.001,'L12 и L21 получены из разных полевых режимов.')
check('positive-definite',L[0][0]>0 and L[1][1]>0 and determinant>0,
    [L[0][0],L[1][1],determinant],0,'Симметричная часть L положительно определена: L11>0,L22>0,det>0.')
for row in RECORDS:
    i=row['currents_A']
    predicted=[sum(L[k][j]*i[j]for j in range(2))for k in range(2)]
    predicted_energy=.5*(L[0][0]*i[0]**2+2*M*i[0]*i[1]+L[1][1]*i[1]**2)
    row['matrix_linkages_Wb']=predicted
    row['matrix_energy_J']=predicted_energy
    row['mutual_energy_term_J']=M*i[0]*i[1]
    flux_error=max(abs(predicted[k]-row['linkages_Wb'][k])for k in range(2))/max(max(abs(v)for v in row['linkages_Wb']),1e-20)
    check(row['tag']+'-superposition',flux_error<1e-7,flux_error,1e-7,'Матрица из единичных токов воспроизводит отдельное поле общего режима.')
    error=relative(predicted_energy,row['energy_J'])
    check(row['tag']+'-matrix-energy',error<.003,error,.003,'Знак взаимного члена M*I1*I2 соответствует знакам токов и положительным N.')
result={'schema':'bul-multicoil-native-v1','implementation':repr(sys.implementation),
    'problem_type':'axisymmetric','materials_used':['Воздух','Обмоточная область'],
    'object':'Две воздушные соосные катушки; магнитопровода нет',
    'circuits':NAMES,'turns':[100,50],'L_matrix_H':L,'M_symmetric_H':M,
    'reciprocity_relative_difference':reciprocity,'determinant_H2':determinant,
    'coupling_factor':M/math.sqrt(L[0][0]*L[1][1]),
    'formula':'Psi=L*I; W=0.5*L11*I1^2+M*I1*I2+0.5*L22*I2^2; M=(L12+L21)/2',
    'sign_convention':'Оба N положительны; M положительна; отрицательный I2 даёт отрицательный взаимный член энергии.',
    'records':RECORDS,'checks':CHECKS,'passed':all(c['passed']for c in CHECKS),
    'scope':'Четыре магнитостатических режима русской копии I04.2; взаимность, суперпозиция и энергия одной линейной воздушной модели. Не проверка магнитопровода или электромеханики.'}
with open(OUT+'/result.json','w')as handle:
    json.dump(result,handle)
print('BUL_MULTICOIL_PASSED',result['passed'],OUT)
if not result['passed']:
    raise RuntimeError('Проверки связанных катушек не выполнены; результаты сохранены')