# FEMM 5, встроенный Python. Откройте учебную PM FEM перед запуском.
# Новый вывод: PM-4А-результаты.json и отдельная рабочая FEM/ANS.
# Статика ХХ; ЭДС восстанавливается из угловой серии потокосцепления.
import json, math, os
COILS=[[-0.1, 65.3, 'A', 1], [12.7, 62.5, 'A', 1], [22.4, 60.0, 'C', -1], [29.7, 55.9, 'C', -1], [48.3, 41.4, 'B', 1], [52.8, 33.3, 'B', 1], [59.5, 22.4, 'B', 1], [60.9, 11.6, 'A', -1], [61.2, 0.4, 'A', -1], [62.0, -9.9, 'A', -1], [59.8, -21.0, 'C', 1], [55.3, -31.9, 'C', 1], [48.9, -40.8, 'C', 1], [40.3, -48.6, 'B', -1], [32.2, -54.7, 'B', -1], [20.5, -58.9, 'B', -1], [11.3, -61.2, 'A', 1], [-1.8, -62.5, 'A', 1], [-11.8, -62.8, 'A', 1], [-21.9, -58.9, 'C', -1], [-33.3, -56.1, 'C', -1], [-40.8, -49.7, 'C', -1], [-47.8, -41.1, 'B', 1], [-55.3, -32.7, 'B', 1], [-59.5, -22.7, 'B', 1], [-63.1, -10.7, 'A', -1], [-62.8, -0.1, 'A', -1], [-65.1, 11.6, 'A', -1], [-62.8, 22.4, 'C', 1], [-56.7, 34.1, 'C', 1], [-51.7, 43.3, 'C', 1], [-42.2, 51.1, 'B', -1], [-32.7, 56.1, 'B', -1], [-21.6, 60.3, 'B', -1], [-9.3, 61.7, 'A', 1], [41.7, 48.9, 'C', -1]]
PHASES=["Фаза A","Фаза B","Фаза C"]
NCOILS_PHASE=6
TARGET_PHASE_RMS_V=220.0
FREQUENCY_ELECTRICAL_HZ=50.0
POLE_PAIRS=2
PERIOD_MECHANICAL_DEG=180.0
MAX_HARMONIC=23
OUTPUT="PM-4А-результаты.json"
def scalar(v):return complex(v).real
def set_turns(turns):
    for x,y,phase,sign in COILS:
        mi_selectlabel(x,y)
        mi_setblockprop("Медь",1,0,"Фаза "+phase,0,10,sign*turns)
        mi_clearselected()
def rotate(delta):
    mi_selectgroup(20)
    mi_moverotate(0,0,delta)
    mi_clearselected()
def solve(name):
    mi_saveas(name+".fem")
    mi_analyze(1)
    mi_loadsolution()
    flux=[scalar(mo_getcircuitproperties(p)[2]) for p in PHASES]
    current=[scalar(mo_getcircuitproperties(p)[0]) for p in PHASES]
    mo_clearblock();mo_groupselectblock(20)
    torque=scalar(mo_blockintegral(22));mo_clearblock()
    row={"flux_Wb":flux,"currents_A":current,"cogging_torque_Nm":torque,"nodes":mo_numnodes(),"elements":mo_numelements()}
    mo_close()
    return row
def spectrum(rows):
    count=len(rows);answer=[]
    for phase in range(3):
        harmonics=[]
        for h in range(1,MAX_HARMONIC+1):
            cr=sum(row["flux_Wb"][phase]*math.cos(2*math.pi*h*k/count) for k,row in enumerate(rows))
            si=-sum(row["flux_Wb"][phase]*math.sin(2*math.pi*h*k/count) for k,row in enumerate(rows))
            amp=2*math.sqrt(cr*cr+si*si)/count
            harmonics.append({"electrical_order":h,"psi_peak_Wb":amp,"emf_rms_V":2*math.pi*FREQUENCY_ELECTRICAL_HZ*h*amp/math.sqrt(2),"flux_phase_rad":math.atan2(si,cr)})
        e1=harmonics[0]["emf_rms_V"]
        answer.append({"phase":PHASES[phase],"fundamental_rms_V":e1,"band_rms_V":math.sqrt(sum(h["emf_rms_V"]**2 for h in harmonics)),"THD_2_to_23":math.sqrt(sum(h["emf_rms_V"]**2 for h in harmonics[1:]))/e1,"harmonics":harmonics})
    return answer
out={"scope":"PM-адаптация статора4А; нулевые токи; серия статических полей, ЭДС при1500об/мин. Нет внешней нагрузки, нагрева, прочности и полного КПД.","target_phase_fundamental_rms_V":TARGET_PHASE_RMS_V,"coil_turns_pilot":1,"phase_turns_pilot":6,"frequency_electrical_Hz":FREQUENCY_ELECTRICAL_HZ,"speed_rpm":1500,"series":[]}
set_turns(1)
for count in [48,96]:
    rows=[];current_angle=0.0
    for k in range(count):
        angle=PERIOD_MECHANICAL_DEG*k/count
        if k:rotate(angle-current_angle)
        current_angle=angle
        row=solve("PM-4А-временное-поле");row["mechanical_angle_deg"]=angle
        rows.append(row)
        out["progress"]={"samples":count,"completed":k+1}
        if (k+1)%12==0:
            with open(OUTPUT,"w") as f:json.dump(out,f)
            print("PM_SERIES",count,k+1)
    rotate(-current_angle)
    item={"samples":count,"rows":rows,"spectrum":spectrum(rows)}
    out["series"].append(item)
    with open(OUTPUT,"w") as f:json.dump(out,f)
fine=out["series"][-1]["spectrum"]
e1_mean=sum(p["fundamental_rms_V"] for p in fine)/3
turns=int(TARGET_PHASE_RMS_V/e1_mean+0.5)
set_turns(turns)
out["turns_per_coil"]=turns;out["phase_turns"]=turns*NCOILS_PHASE
out["selected_spectrum"]=[]
for p in fine:
    out["selected_spectrum"].append({"phase":p["phase"],"fundamental_rms_V":turns*p["fundamental_rms_V"],"band_rms_V":turns*p["band_rms_V"],"THD_2_to_23":p["THD_2_to_23"]})
out["angular_refinement_relative"]=[abs(out["series"][0]["spectrum"][i]["fundamental_rms_V"]-fine[i]["fundamental_rms_V"])/fine[i]["fundamental_rms_V"] for i in range(3)]
out["scaling_verification"]=[]
for k in [0,8,16]:
    angle=PERIOD_MECHANICAL_DEG*k/96
    if angle:rotate(angle)
    row=solve("PM-4А-проверка-витков")
    if angle:rotate(-angle)
    reference=out["series"][-1]["rows"][k]["flux_Wb"]
    scale=max(abs(turns*v) for v in reference)
    error=max(abs(row["flux_Wb"][j]-turns*reference[j]) for j in range(3))/scale
    out["scaling_verification"].append({"mechanical_angle_deg":angle,"relative_error":error,"native":row,"reference_pilot_flux_Wb":reference})
# Финальная сохранённая FEM/ANS относится к выбранному целому числу витков.
final=solve("PM-генератор-4А-Выбранная-обмотка")
mi_close()
opendocument("PM-генератор-4А-Выбранная-обмотка.fem")
mi_loadsolution()
reopened=[scalar(mo_getcircuitproperties(p)[2]) for p in PHASES]
mo_close();mi_close()
out["final_static"]=final
out["final_solution_reopened_flux_Wb"]=reopened
out["passed"]=max(out["angular_refinement_relative"])<.005 and max(r["relative_error"] for r in out["scaling_verification"])<.001 and max(abs(reopened[j]-final["flux_Wb"][j]) for j in range(3))<1e-10
with open(OUTPUT,"w") as f:json.dump(out,f)
print("PM_COMPLETE",out["passed"],turns,out["selected_spectrum"],out["angular_refinement_relative"])