"""Геометрическая оптимизация SPM: программа встроенного Python FEMM 5.

Запуск: File → Run Script. Полная серия продолжительна; результаты пишутся
рядом с программой. Для ознакомления аргумент smoke рассчитывает один угол.
Все конструкции строятся заново; внешние пакеты и исходные FEM не нужны.
"""
import math, json, os, sys, time
assert sys.implementation.name == 'micropython', 'Нужен встроенный Python FEMM 5'
_script = globals().get('__file__') or (sys.argv[0] if sys.argv else '')
if not _script or not os.path.isfile(_script):
    raise RuntimeError('Запустите сохранённый файл из полного комплекта')
_program_dir=os.path.dirname(os.path.abspath(_script))
os.chdir(_program_dir)
_output_dir=globals().get('FEMM_OPT_OUTPUT_DIR','')
if not _output_dir:
    try: os.mkdir('runs')
    except OSError: pass
    counter=1
    while os.path.exists('runs/расчёт-'+str(counter)):
        counter+=1
    _output_dir=os.path.join(_program_dir,'runs','расчёт-'+str(counter))
    os.mkdir(_output_dir)
os.chdir(_output_dir)
print('Новый каталог результатов',os.getcwd())
for folder in ('calculated', 'models'):
    try: os.mkdir(folder)
    except OSError: pass
MU0 = 4e-7 * math.pi
P = 2
DEPTH = 80.0
RO = 57.0
CURRENT = 10.0
N_SEARCH = 12
MIN_T = 8.0
MAX_B = 1.95
CONSTRAINT_VERSION = 'учебная-нагрузка-8-Нм-v2'
nodes = {}; edges = set(); arcs = set(); angles = []
field_count = 0

def dump(path, value):
    with open(path, 'w') as handle: json.dump(value, handle)

def polar(r, angle):
    return r*math.cos(math.radians(angle)), r*math.sin(math.radians(angle))

def node(x, y):
    key = (round(x, 7), round(y, 7))
    if key not in nodes: mi_addnode(x, y); nodes[key] = 1

def line(x1, y1, x2, y2):
    node(x1, y1); node(x2, y2)
    key = tuple(sorted(((round(x1, 7), round(y1, 7)), (round(x2, 7), round(y2, 7)))))
    if key not in edges: mi_addsegment(x1, y1, x2, y2); edges.add(key)

def arc(r, a, b, fine):
    a = a % 360; span = (b-a) % 360
    if span < 1e-7: return
    cuts = [a] + sorted([a+(t-a)%360 for t in angles if 1e-6 < (t-a)%360 < span-1e-6]) + [a+span]
    for left, right in zip(cuts[:-1], cuts[1:]):
        x1, y1 = polar(r, left); x2, y2 = polar(r, right)
        node(x1, y1); node(x2, y2)
        key = (round(r, 7), round(left%360, 7), round(right%360, 7))
        if key not in arcs: mi_addarc(x1, y1, x2, y2, right-left, 2.5 if fine else 4); arcs.add(key)

def circle(r, fine):
    cuts = sorted(set([0, 90, 180, 270] + [round(a%360, 7) for a in angles]))
    for a, b in zip(cuts, cuts[1:]+[cuts[0]+360]): arc(r, a, b, fine)

def sector(r1, r2, a, b, fine):
    arc(r1, a, b, fine); arc(r2, a, b, fine)
    x1, y1 = polar(r1, a); x2, y2 = polar(r2, a); line(x1, y1, x2, y2)
    x1, y1 = polar(r1, b); x2, y2 = polar(r2, b); line(x1, y1, x2, y2)

def label(x, y, material, circuit='', turns=0, mesh=2, group=0, direction=0):
    mi_addblocklabel(x, y); mi_selectlabel(x, y)
    mi_setblockprop(material, 0, mesh, circuit, direction, group, turns)
    mi_clearselected()

def close():
    for name in ('mo_close', 'mi_close'):
        try: globals()[name]()
        except Exception: pass

def build(full_arc, thickness, theta, current, gamma, fine=False):
    global nodes, edges, arcs, angles
    close(); nodes = {}; edges = set(); arcs = set()
    half = full_arc/2; ri = RO-thickness
    angles = [0, 90, 180, 270]
    for k in range(12): angles.extend([k*30-7, k*30+7])
    for k in range(4): angles.extend([theta+k*90-half, theta+k*90+half])
    newdocument(0); mi_probdef(0, 'millimeters', 'planar', 1e-8, DEPTH, 25)
    mi_addmaterial('Воздух', 1, 1, 0, 0, 0)
    mi_addmaterial('Обмотка — идеальный ток', 1, 1, 0, 0, 0)
    mi_addmaterial('Электротехническая сталь — учебная B–H', 1, 1, 0, 0, 0)
    for h, b in [(0,0),(100,.25),(200,.5),(400,.9),(800,1.2),(2000,1.45),(5000,1.6),(10000,1.7),(20000,1.8),(50000,1.95),(100000,2.1),(400000,2.5)]:
        mi_addbhpoint('Электротехническая сталь — учебная B–H', b, h)
    mi_addmaterial('Постоянный магнит — Br 1,2 Тл', 1.05, 1.05, 1.2/(MU0*1.05), 0, 0)
    mi_addboundprop('Внешняя граница A=0', 0,0,0,0,0,0,0,0,0)
    side=170
    for a,b,c,d in [(-side,-side,side,-side),(side,-side,side,side),(side,side,-side,side),(-side,side,-side,-side)]: line(a,b,c,d)
    for x,y in [(0,-side),(side,0),(0,side),(-side,0)]:
        mi_selectsegment(x,y); mi_setsegmentprop('Внешняя граница A=0',0,1,0,0); mi_clearselected()
    label(144,144,'Воздух',mesh=8)
    mapping=[('Фаза A',1),('Фаза C',-1),('Фаза B',1),('Фаза A',-1),('Фаза C',1),('Фаза B',-1)]
    for c,a in [('Фаза A',0),('Фаза B',-120),('Фаза C',120)]:
        mi_addcircprop(c,current*math.cos(math.radians(gamma+a)),1)
    circle(60,fine); circle(100,fine)
    label(0,94,'Электротехническая сталь — учебная B–H',mesh=1.8 if fine else 3,group=1)
    for k in range(12):
        center=k*30; sector(60,85,center-7,center+7,fine)
        x,y=polar(73,center); c,sign=mapping[k%6]
        label(x,y,'Обмотка — идеальный ток',c,sign*40,mesh=1.5 if fine else 2.5,group=1)
    circle(18,fine); circle(ri,fine); circle(RO,fine)
    label(*polar((18+ri)/2,theta+45),'Электротехническая сталь — учебная B–H',mesh=1.5 if fine else 2.5,group=2)
    for k in range(4):
        center=theta+k*90; sector(ri,RO,center-half,center+half,fine)
        label(*polar((ri+RO)/2,center),'Постоянный магнит — Br 1,2 Тл',mesh=1 if fine else 1.8,group=2,direction='theta' if k%2==0 else 'theta+180')
        label(*polar((ri+RO)/2,center+45),'Воздух',mesh=.7 if fine else 1.3)
    label(0,0,'Воздух',mesh=1.5 if fine else 2.5)
    label(*polar(58.5,17),'Воздух',mesh=.3 if fine else .6)
    mi_selectcircle(0,0,58.01,4); mi_setgroup(2); mi_clearselected()

def solve(full_arc, thickness, theta, current, gamma, fine=False, retained=None):
    global field_count
    build(full_arc,thickness,theta,current,gamma,fine)
    target='models/'+retained if retained else 'calculated/current'
    mi_saveas(target+'.fem'); mi_analyze(1); mi_loadsolution(); mo_smooth('off')
    field_count += 1
    mo_clearblock(); mo_groupselectblock(2); torque=float(mo_blockintegral(22)); mo_clearblock()
    flux=[float(mo_getcircuitproperties(n)[2]) for n in ('Фаза A','Фаза B','Фаза C')]
    # Контроль в фиксированных внутренних точках; это не максимум в сингулярной вершине.
    bmax=0.; hparallel=1e30
    for k in range(12):
        for r in (62,68,78,90,95):
            v=mo_getpointvalues(*polar(r,k*30+15)); bmax=max(bmax,math.sqrt(float(v[1])**2+float(v[2])**2))
    for k in range(4):
        center=theta+k*90; a=math.radians(center+(180 if k%2 else 0))
        for fraction in (.2,.5,.8):
            v=mo_getpointvalues(*polar(RO-thickness+fraction*thickness,center))
            hparallel=min(hparallel,float(v[5])*math.cos(a)+float(v[6])*math.sin(a))
    assert all(math.isfinite(v) for v in [torque,bmax,hparallel]+flux)
    return {'theta_deg':theta,'torque_Nm':torque,'flux_Wb':flux,'B_control_max_T':bmax,'H_parallel_min_Am':hparallel,'elements':int(mo_numelements())}

def volume_cm3(full_arc,thickness):
    # Четыре кольцевых сектора; мм³ → см³.
    return 4*.5*math.radians(full_arc)*(RO*RO-(RO-thickness)**2)*DEPTH/1000

def evaluate(full_arc,thickness,n=N_SEARCH,fine=False,current=CURRENT,stage='coarse',retain=None,period=30.,shift=0.):
    zero=solve(full_arc,thickness,0,0,0,fine)
    psi=zero['flux_Wb']; alpha=2/3*(psi[0]-.5*psi[1]-.5*psi[2]); beta=math.sqrt(3)/3*(psi[1]-psi[2])
    offset=math.degrees(math.atan2(beta,alpha))
    rows=[]
    for k in range(n):
        theta=shift+period*k/n
        rows.append(solve(full_arc,thickness,theta,current,2*theta+offset+90,fine,retain if k==0 else None))
    values=[r['torque_Nm'] for r in rows]
    result={'id':'a'+str(round(full_arc,4))+'-t'+str(round(thickness,4)), 'arc_deg':full_arc,'thickness_mm':thickness,
        'stage':stage,'fine_mesh':fine,'angular_samples':n,'angular_period_deg':period,'angular_shift_deg':shift,'current_peak_A':current,'phase_offset_deg':offset,
        'mean_torque_Nm':sum(values)/n,'ripple_pp_Nm':max(values)-min(values),'magnet_volume_cm3':volume_cm3(full_arc,thickness),
        'B_control_max_T':max(r['B_control_max_T'] for r in rows),'H_parallel_min_Am':min(r['H_parallel_min_Am'] for r in rows),
        'torque_series':rows,'field_solutions':n+1}
    result['ripple_relative_pct']=100*result['ripple_pp_Nm']/max(abs(result['mean_torque_Nm']),1e-12)
    result['feasible']=result['mean_torque_Nm']>=MIN_T and result['B_control_max_T']<=MAX_B
    result['constraints']={'mean_torque_margin_Nm':result['mean_torque_Nm']-MIN_T,'B_control_margin_T':MAX_B-result['B_control_max_T']}
    print('Оценена конструкция', result['id'], result['mean_torque_Nm'],result['ripple_pp_Nm'],result['B_control_max_T'],result['feasible'])
    return result

def criteria(row):
    return [-row['mean_torque_Nm'],row['ripple_pp_Nm'],row['magnet_volume_cm3']]

def dominates(a,b):
    aa=criteria(a); bb=criteria(b)
    return all(x<=y for x,y in zip(aa,bb)) and any(x<y for x,y in zip(aa,bb))

def mark_pareto(rows):
    valid=[r for r in rows if r['feasible']]
    for row in rows: row['pareto']=row['feasible'] and not any(dominates(other,row) for other in valid)
    return [r for r in rows if r['pareto']]

def select_compromise(rows):
    valid=[r for r in rows if r['feasible']]
    assert valid,'Нет допустимых конструкций: измените постановку, не подменяя результаты'
    values=[criteria(r) for r in valid]
    lo=[min(x[j] for x in values) for j in range(3)]; hi=[max(x[j] for x in values) for j in range(3)]
    for row in valid:
        row['normalized_criteria']=[(criteria(row)[j]-lo[j])/max(hi[j]-lo[j],1e-12) for j in range(3)]
        # Явно заданные предпочтения конструктора: момент 0,4; пульсации 0,4; PM 0,2.
        row['compromise_score']=sum(w*z for w,z in zip((.4,.4,.2),row['normalized_criteria']))
    return min(valid,key=lambda r:r['compromise_score']), {'weights':[.4,.4,.2],'ideal':lo,'nadir_sample':hi}

def write_csv(rows,path):
    keys=['id','stage','arc_deg','thickness_mm','current_peak_A','angular_samples','angular_period_deg','angular_shift_deg','fine_mesh','mean_torque_Nm','ripple_pp_Nm','ripple_relative_pct','magnet_volume_cm3','B_control_max_T','H_parallel_min_Am','feasible','pareto']
    with open(path,'w') as f:
        f.write(','.join(keys)+'\n')
        for row in rows: f.write(','.join(str(row.get(k,'')) for k in keys)+'\n')

def main():
    mode=sys.argv[1] if len(sys.argv)>1 else 'full'
    if mode=='smoke':
        row=solve(72,8,0,CURRENT,180,False,'проверка-русских-имён')
        dump('calculated/smoke.json',row); close(); print('SMOKE OK',row); return
    started=time.ticks_ms(); rows=[]
    # Начальный план: полный факторный набор 4 × 4.
    for a in (60.,68.,76.,84.):
        for t in (5.,6.5,8.,9.5):
            rows.append(evaluate(a,t)); dump('calculated/checkpoint.json',rows)
    mark_pareto(rows); chosen,normalization=select_compromise(rows)
    # Локальное уточнение действительно зависит от рассчитанного компромисса.
    center=(chosen['arc_deg'],chosen['thickness_mm']); seen=set((r['arc_deg'],r['thickness_mm']) for r in rows)
    for da,dt in [(-4,0),(4,0),(0,-.75),(0,.75),(-4,-.75),(-4,.75),(4,-.75),(4,.75)]:
        a=min(84.,max(60.,center[0]+da)); t=min(9.5,max(5.,center[1]+dt))
        if (a,t) not in seen:
            seen.add((a,t)); rows.append(evaluate(a,t,stage='adaptive')); dump('calculated/checkpoint.json',rows)
    pareto=mark_pareto(rows); chosen,normalization=select_compromise(rows)
    best=max([r for r in rows if r['feasible']],key=lambda r:r['mean_torque_Nm'])
    smallest=min(pareto,key=lambda r:r['magnet_volume_cm3'])
    quietest=min(pareto,key=lambda r:r['ripple_pp_Nm'])
    selection={'compromise':chosen['id'],'max_mean_torque':best['id'],'min_magnet_volume':smallest['id'],'min_ripple':quietest['id']}
    baseline=next(r for r in rows if r['arc_deg']==76 and r['thickness_mm']==8)
    verification=[]
    # Разделяем угловое и пространственное уточнение. Для сравнения используется
    # одна и та же плотная угловая сетка из 24 точек.
    unique=[]
    for row in (baseline,chosen,best):
        if row['id'] not in [r['id'] for r in unique]: unique.append(row)
    for row in unique:
        for fine in (False,True):
            verification.append(evaluate(row['arc_deg'],row['thickness_mm'],n=24,fine=fine,stage='verification',retain=('компромисс' if row['id']==chosen['id'] else 'максимальный-момент' if row['id']==best['id'] else 'базовая-конструкция')+('-уточнённая' if fine else '-основная')))
    # Полный электрический период составляет 360/p = 180 механических градусов.
    # Смещение на 1,25° проверяет также точки между исходными поисковыми углами.
    verification.append(evaluate(chosen['arc_deg'],chosen['thickness_mm'],n=72,fine=True,stage='full-electrical-period',period=180.,shift=1.25,retain='компромисс-полный-электрический-период'))
    # Проверка чувствительности компромисса к снижению тока: отдельные поля.
    verification.append(evaluate(chosen['arc_deg'],chosen['thickness_mm'],n=24,fine=True,current=7.5,stage='second-current',retain='компромисс-ток-7_5-А'))
    # Сохраняем отдельные нативные модели крайних точек фронта, не все временные поля.
    for name,row in [('Парето-минимальный-объём',smallest),('Парето-минимальные-пульсации',quietest)]:
        solve(row['arc_deg'],row['thickness_mm'],0,CURRENT, row['phase_offset_deg']+90,False,name)
    close()
    result={'schema':'femm5-optimization-machine-v1','runtime':'Встроенный MicroPython FEMM 5',
        'model':{'topology':'SPM, 12 пазов, 4 полюса','stack_length_mm':DEPTH,'stator_inner_radius_mm':60,'stator_outer_radius_mm':100,'rotor_outer_radius_mm':RO,'air_gap_mm':3,'shaft_radius_mm':18,'slot_depth_mm':25,'slot_arc_deg':14,'turns_per_slot':40,'steel':'Заданная в программе нелинейная учебная B–H','magnet_Br_T':1.2,'magnet_mu_r':1.05,'frequency_Hz':0},
        'search':{'algorithm':'16 начальных точек; локальные соседи выбранного компромисса с шагами 4° и 0,75 мм','adaptive_center_arc_deg':center[0],'adaptive_center_thickness_mm':center[1],'domain':{'arc_deg':[60,84],'thickness_mm':[5,9.5]},'current_peak_A':CURRENT,'angular_period_deg':30,'angular_samples':N_SEARCH,'constraints':{'mean_torque_min_Nm':MIN_T,'B_control_max_T':MAX_B,'version':CONSTRAINT_VERSION},'objectives':['max mean_torque_Nm','min ripple_pp_Nm','min magnet_volume_cm3'],'normalization':normalization},
        'evaluations':rows,'pareto_ids':[r['id'] for r in pareto],'selection':selection,'baseline_id':baseline['id'],'verification':verification,
        'field_solutions':field_count,'elapsed_ms':time.ticks_diff(time.ticks_ms(),started),
        'scope':['Дискретный фронт среди рассчитанных допустимых конструкций; глобальная или непрерывная оптимальность не доказана.','Квазистатические магнитные положения с заданными синусоидальными токами; нет движения во времени, преобразователя, потерь и нагрева.','Учебная характеристика стали не заменяет паспорт конкретной марки; Br и проницаемость магнита заданы при одной условной температуре.','Пиковая B ограничена в фиксированных внутренних контрольных точках стали, не в вершинах зубцов и не по всему полю.','H вдоль намагниченности — диагностическая величина; необратимое размагничивание не проверено без температурных диаграмм магнита.','Внешняя граница A=0 находится на квадрате ±170 мм; торцевые эффекты, скос, мехпрочность и технологические допуски не учтены.']}
    dump('machine-optimization-result.json',result); write_csv(rows,'machine-evaluations.csv'); write_csv(verification,'machine-verification.csv')
    print('OPTIMIZATION COMPLETE',len(rows),len(pareto),field_count,selection)

main()