"""Конический втяжной электромагнит: встроенный Python FEMM 5.

Собственная учебная геометрия по методическим идеям книги О. Б. Буля,
не копия П5. Запустите Run Python Script в FEMM 5: создаётся новый
calculated-NN рядом с программой. 12 отдельных DC-решений, без интерфейсных
снимков. CONICAL_OUTPUT_ROOT задаёт существующий каталог результатов.
"""
import os
import sys
import math
import json

script_path = globals().get('__file__', '') or (sys.argv[0] if sys.argv else '')
BASE = globals().get('CONICAL_OUTPUT_ROOT', os.path.dirname(os.path.abspath(script_path)) if script_path else os.getcwd())
counter = 1
while os.path.exists(BASE+'/calculated-%02d' % counter):
    counter += 1
OUT = BASE+'/calculated-%02d' % counter
os.mkdir(OUT)
MU0 = 4e-7*math.pi
N = 300
AIR = 'Воздух'
COIL = 'Обмотка возбуждения — немагнитная'
LINEAR = 'Сталь учебная — мю 2000'
NONLINEAR = 'Сталь учебная — условная B(H)'
CIRCUIT = 'Катушка втяжного привода'
BOUNDARY = 'Внешняя граница A=0'
GROUPS = {'1':'Неподвижная чашка и стоп','2':'Конический якорь','3':'Обмотка возбуждения','4':'Воздушная область'}
CORE = [[0,0],[30,0],[30,55],[22,55],[22,8],[8,8],[8,35],[4,43],[0,43]]
BH = [[0,0],[.4,.4/(MU0*2000)],[.8,.8/(MU0*2000)],[1.2,600],
      [1.45,1200],[1.6,2600],[1.75,6000],[1.9,18000],[2.1,60000],[2.3,180000]]
CONTRACT = {'revision':'bul-independent-conical-v1','problem_type':'axisymmetric',
    'coordinates':'r,z in mm; r >= 0; z upward; full 360-degree body',
    'method_source':'О. Б. Буль, 2005, П5, с.297–318: методика; размеры и B(H) самостоятельные.',
    'core_polygon_mm':CORE,'armature_polygon_mm_at_gap_d':[[0,'43+d'],[4,'43+d'],[8,'35+d'],[20.5,'35+d'],[20.5,'53+d'],[0,'53+d']],
    'winding_rectangle_mm':[10,10,19,32],'turns':N,'groups':GROUPS,
    'linear_mu_r':2000,'BH_T_A_per_m':BH,
    'BH_origin':'Условная монотонная безгистерезисная кривая; начальная мю=2000; не паспорт марки стали.',
    'working_geometry':{'d_mm':'Осевое открытие между центральными плоскими площадками',
        'fixed_cone':'z+2r=51, r4..8','moving_cone':'z+2r=51+d, r4..8',
        'normal_cone_gap_mm':'d/sqrt(5), в области перекрытия наклонных поверхностей',
        'radial_return_gap_mm':1.5,'radial_return_radii_mm':[20.5,22],
        'radial_overlap_mm':'min(55,53+d)-(35+d); 18mm for d<=2, 20-d for d>=2'},
    'domain_mm':[0,100,-50,120],
    'force_sign':'Fz<0 притягивает якорь вниз; d увеличивается при движении вверх; Fz=dWprime/dd при фиксированном I.',
    'omissions':['механические направляющие и шток','трение и пружина','нагрев','гистерезис','eddy currents','переходный режим','промышленная проверка']}
with open(OUT+'/geometry-contract.json','w') as f:
    json.dump(CONTRACT,f)
RECORDS = []
CHECKS = []
NODES = {}
SEGMENTS = {}


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


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


def armature(d):
    return [[0,43+d],[4,43+d],[8,35+d],[20.5,35+d],[20.5,53+d],[0,53+d]]


def revolved_volume(points):
    # 2*pi*integral(r dA): first area moment of the rz polygon, mm^3 -> m^3.
    moment = 0
    for k in range(len(points)):
        a,b = points[k],points[(k+1)%len(points)]
        moment += (a[0]+b[0])*(a[0]*b[1]-b[0]*a[1])/6
    return abs(2*math.pi*moment)*1e-9


def node(p):
    key = (round(p[0],9),round(p[1],9))
    if key not in NODES:
        mi_addnode(key[0],key[1])
        NODES[key] = True
    return key


def segment(a,b,mesh,group=0,boundary='<None>'):
    a,b = node(a),node(b)
    key = tuple(sorted((a,b)))
    if key not in SEGMENTS:
        mi_addsegment(a[0],a[1],b[0],b[1])
        SEGMENTS[key] = True
    mi_selectsegment((a[0]+b[0])/2,(a[1]+b[1])/2)
    mi_setsegmentprop(boundary,mesh,0,0,group)
    mi_clearselected()


def polygon(points,mesh,group):
    for k in range(len(points)):
        segment(points[k],points[(k+1)%len(points)],mesh,group)


def label(r,z,material,mesh,group,circuit='<None>',turns=0):
    mi_addblocklabel(r,z)
    mi_selectlabel(r,z)
    mi_setblockprop(material,0,mesh,circuit,0,group,turns)
    mi_clearselected()


def build(d,current,nonlinear,mesh_factor,domain_factor):
    global NODES,SEGMENTS
    NODES,SEGMENTS = {},{}
    newdocument(0)
    mi_probdef(0,'millimeters','axi',1e-9,1,25,0)
    mi_setcomment('Самостоятельный конический втяжной привод для курса FEMM 5. Не копия П5 книги Буля. '+
        'Группы: 1 — неподвижная чашка и стоп; 2 — конический якорь; 3 — катушка; 4 — воздух. '+
        'Осевое открытие d=%g мм; нормальный зазор конуса d/sqrt(5); радиальный возвратный промежуток 1.5 мм. '%d+
        'Перекрытие возврата меняется с d. Ток %g А; N=300. '%current+
        ('Условная B(H), без гистерезиса.' if nonlinear else 'Учебная линейная сталь, мю=2000.'))
    for name,mu in [(AIR,1),(COIL,1),(LINEAR,2000),(NONLINEAR,2000)]:
        mi_addmaterial(name,mu,mu,0,0,0)
    for b,h in BH:
        mi_addbhpoint(NONLINEAR,b,h)
    mi_addcircprop(CIRCUIT,current,1)
    mi_addboundprop(BOUNDARY,0,0,0,0,0,0,0,0,0)
    hs = .65*mesh_factor
    hg = min(.14,d/(math.sqrt(5)*6))*mesh_factor
    ha = 2.0*mesh_factor
    polygon(CORE,hs,1)
    moving = armature(d)
    polygon(moving,hs,2)
    polygon([[10,10],[19,10],[19,32],[10,32]],hs,3)
    steel = NONLINEAR if nonlinear else LINEAR
    label(26,25,steel,hs,1)
    label(10,48+d,steel,hs,2)
    label(14,22,COIL,hs,3,CIRCUIT,N)
    label(55,20,AIR,ha,4)
    # All working surfaces are refined using their shortest normal spacing.
    for r,z,group in [(2,43,1),(6,39,1),(2,43+d,2),(6,39+d,2)]:
        mi_selectsegment(r,z)
        mi_setsegmentprop('<None>',hg,0,0,group)
        mi_clearselected()
    for r,z,group in [(22,45,1),(20.5,44+d,2)]:
        mi_selectsegment(r,z)
        mi_setsegmentprop('<None>',min(.25,hs),0,0,group)
        mi_clearselected()
    rmax,zmin,zmax = 100*domain_factor,-50*domain_factor,120*domain_factor
    cuts = [zmin,0,43,43+d,53+d,zmax]
    for k in range(len(cuts)-1):
        segment((0,cuts[k]),(0,cuts[k+1]),hg if k==2 else hs,4,BOUNDARY)
    for a,b in [((0,zmin),(rmax,zmin)),((rmax,zmin),(rmax,zmax)),((rmax,zmax),(0,zmax))]:
        segment(a,b,ha,4,BOUNDARY)
    return {'steel_mm':hs,'normal_gap_edge_mm':hg,'air_mm':ha}


def flux(a,b):
    mo_clearcontour()
    mo_addcontour(a[0],a[1])
    mo_addcontour(b[0],b[1])
    result = real(mo_lineintegral(0)[0])
    mo_clearcontour()
    return result


def measure(d):
    current,voltage,linkage = mo_getcircuitproperties(CIRCUIT)
    mo_clearblock()
    mo_groupselectblock(2)
    force_z = real(mo_blockintegral(19))
    volume = real(mo_blockintegral(10))
    mo_clearblock()
    mo_groupselectblock()
    energy = real(mo_blockintegral(2))
    coenergy = real(mo_blockintegral(17))
    mo_clearblock()
    probes = []
    for name,r,z in [('Рабочая плоская площадка',2,43+d/2),('Рабочий конус',6+.2*d,39+.1*d),
        ('Центральный стоп',4,25),('Наружный возврат',26,25),('Якорь',10,48+d)]:
        values = mo_getpointvalues(r,z)
        probes.append({'label':name,'point_rz_mm':[r,z],'B_rz_T':[real(values[1]),real(values[2])],
            'B_abs_T':math.sqrt(real(values[1])**2+real(values[2])**2),
            'H_rz_A_per_m':[real(values[5]),real(values[6])]})
    current,linkage = real(current),real(linkage)
    exact = revolved_volume(armature(d))
    return {'current_A':current,'linkage_Wb':linkage,'static_inductance_H':linkage/current,
        'force_z_N':force_z,'attraction_N':-force_z,'energy_J':energy,'coenergy_J':coenergy,
        'energy_identity_relative_error':abs(energy+coenergy-current*linkage)/max(abs(current*linkage),1e-20),
        'armature_volume_m3':volume,'armature_volume_exact_m3':exact,'armature_volume_relative_error':relative(volume,exact),
        'central_core_section_flux_Wb':flux((.001,25),(7.999,25)),
        'outer_return_section_flux_Wb':flux((22.001,25),(29.999,25)),
        'section_scope':'Ориентированные интегралы Bn по двум названным сечениям; не весь рабочий поток и не универсальный коэффициент утечки.',
        'field_probes':probes,'elements':mo_numelements(),'nodes':mo_numnodes()}


def solve(tag,d=1,current=1,nonlinear=False,mesh_factor=1,domain_factor=1):
    mesh = build(d,current,nonlinear,mesh_factor,domain_factor)
    stem = '%02d-Конический-привод-%s'%(len(RECORDS)+1,tag)
    path = OUT+'/'+stem+'.fem'
    mi_saveas(path)
    mi_close()
    opendocument(path)
    mi_analyze(1)
    mi_loadsolution()
    row = {'tag':tag,'model':stem+'.fem','solution':stem+'.ans','axial_opening_mm':d,
        'normal_cone_gap_mm':d/math.sqrt(5),'radial_return_gap_mm':1.5,
        'radial_return_overlap_mm':min(55,53+d)-(35+d),'current_A':current,'turns':N,
        'nonlinear':nonlinear,'mesh_factor':mesh_factor,'domain_factor':domain_factor,
        'mesh':mesh,'saved_reopened_solved':True}
    row.update(measure(d))
    mo_close()
    mi_close()
    opendocument(path)
    mi_loadsolution()
    reopened = real(mo_getcircuitproperties(CIRCUIT)[2])
    row['solution_reopened_linkage_Wb'] = reopened
    row['solution_reopened_matches'] = abs(reopened-row['linkage_Wb'])<1e-12
    mo_close()
    mi_close()
    RECORDS.append(row)
    with open(OUT+'/progress.json','w') as f:
        json.dump({'completed_solves':len(RECORDS),'last':stem},f)
    print('CONICAL_SOLVED',tag,row['elements'],row['attraction_N'])
    return row


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


base = solve('ход-1мм')
gap_series = [solve('ход-0.5мм',d=.5),base,solve('ход-2мм',d=2),solve('ход-4мм',d=4)]
fine = solve('мелкая-сетка',mesh_factor=.6)
large = solve('большая-область',mesh_factor=.6,domain_factor=1.5)
check('mesh-force',relative(base['attraction_N'],fine['attraction_N'])<.025,
    relative(base['attraction_N'],fine['attraction_N']),.025,'Две локальные сетки, одинаковая геометрия и область; контрольная точка d=1мм,I=1А.')
check('domain-force',relative(large['attraction_N'],fine['attraction_N'])<.015,
    relative(large['attraction_N'],fine['attraction_N']),.015,'Внешний прямоугольник увеличен в1.5раза при той же локальной сетке; d=1мм,I=1А.')
DERIVATIVES = []
for h in [.1,.05]:
    minus = solve('коэнергия-минус-%g'%h,d=1-h,mesh_factor=.6)
    plus = solve('коэнергия-плюс-%g'%h,d=1+h,mesh_factor=.6)
    force = (minus['coenergy_J']-plus['coenergy_J'])/(2*h*.001)
    DERIVATIVES.append({'step_mm':h,'attraction_from_coenergy_N':force,'weighted_Maxwell_N':fine['attraction_N'],
        'relative_difference':relative(force,fine['attraction_N']),'minus_model':minus['model'],'plus_model':plus['model'],
        'motion_scope':'Весь якорь смещается по z; радиальная толщина возврата постоянна, перекрытие соответствует геометрии.'})
check('force-methods',DERIVATIVES[-1]['relative_difference']<.04,DERIVATIVES[-1]['relative_difference'],.04,
    'Взвешенный тензор Maxwell и -dWprime/dd при фиксированном токе; d=1мм,I=1А.')
check('derivative-step',relative(DERIVATIVES[0]['attraction_from_coenergy_N'],DERIVATIVES[1]['attraction_from_coenergy_N'])<.025,
    relative(DERIVATIVES[0]['attraction_from_coenergy_N'],DERIVATIVES[1]['attraction_from_coenergy_N']),.025,'Центральные шаги0.1и0.05мм; точка d=1мм не лежит на изломе зависимости перекрытия.')
for current in [3,6]:
    solve('условное-насыщение-%gА'%current,current=current,nonlinear=True,mesh_factor=.6)
for row in RECORDS:
    check(row['model']+'-roundtrip',row['solution_reopened_matches'],row['solution_reopened_linkage_Wb'],row['linkage_Wb'],
        'Русское имя цепи сохранено при повторном открытии точной FEM/ANS-пары.')
    check(row['model']+'-volume',row['armature_volume_relative_error']<1e-8,row['armature_volume_relative_error'],1e-8,
        'Нативный объём выбранной группы якоря и аналитический объём тела вращения.')
    check(row['model']+'-energy',row['energy_identity_relative_error']<.003,row['energy_identity_relative_error'],.003,
        'DC безгистерезисный контроль W+Wprime=I*Psi.')
    check(row['model']+'-force-sign',row['attraction_N']>0,row['attraction_N'],0,'Якорь притягивается вниз.')
check('monotone-stroke',all(gap_series[k]['attraction_N']>gap_series[k+1]['attraction_N'] for k in range(len(gap_series)-1)),
    [r['attraction_N'] for r in gap_series],None,'Только рассчитанная линейная серия d0.5/1/2/4мм; не универсальный закон любых конических приводов.')
result = {'schema':'bul-independent-conical-native-v1','implementation':repr(sys.implementation),
    'geometry_revision':CONTRACT['revision'],'records':RECORDS,'force_derivatives':DERIVATIVES,'checks':CHECKS,
    'passed':all(c['passed'] for c in CHECKS),'scope':'12 DC-решений собственного конического осесимметричного привода. Сетка/область/производная силы проверены только d1мм,I1А,линейнаямю2000. Нелинейные I3/6А демонстрационные; их сила отдельно не аттестована производной.'}
with open(OUT+'/result.json','w') as f:
    json.dump(result,f)
print('CONICAL_PASSED',result['passed'],len(RECORDS),OUT)
if not result['passed']:
    raise RuntimeError('Конический привод: критерии не выполнены; все результаты сохранены отдельно')