# -*- coding: utf-8 -*-
"""
对标 V3.0 解析报告的 FEA 验证扫描
==================================
参照《轴向磁通电机轴向磁拉力计算与轴承选型校核报告V3.0-20260826.pdf》:
  - 其基准: 磁钢 20°C, 空载 Fz=483 N (中值口径 500 N)
  - 其表3-3: g=0.6/1.0/1.5 mm -> 601/483/378 N, kneg(1mm)=250 N/mm
本模型 .mot 磁钢温度默认 100°C (Br -0.12%/K), 首先归一到 20°C 再扫气隙。

工况: 磁钢 20°C x 气隙 {0.6, 1.0, 1.5} mm, 各求解一次, 读 OC/OL 净轴向力
及气隙磁密; 有限差分求磁负刚度 kneg。

用法: python axial_compare.py [--quit]
"""
import json
import os
import sys
import time

BASE = os.path.dirname(os.path.abspath(__file__))
MOT_SRC = os.path.join(BASE, "MARS-12S10P_SSSR_D76-C150_V5.0-0819.mot")
OUT_DIR = os.path.join(BASE, "output_motorcad")

GAPS = [0.6, 1.0, 1.5]
MAGNET_TEMP_C = 20.0
MAX_TSTEPS = 64
MAX_NODES = 40


def stats(ys):
    if not ys:
        return None
    return {"mean": sum(ys) / len(ys), "min": min(ys), "max": max(ys),
            "pk2pk": max(ys) - min(ys), "n": len(ys)}


def read_nodes(mc, graph, sec, tstep):
    xs, ys = [], []
    for i in range(MAX_NODES):
        try:
            x, y = mc.get_magnetic_3d_graph_point(graph, sec, i, tstep)
        except Exception:
            break
        xs.append(x)
        ys.append(y)
    return xs, ys


def net_force_series(mc, graph):
    series = []
    for tstep in range(MAX_TSTEPS):
        total, got = 0.0, False
        for sec in (1, 2):
            xs, ys = read_nodes(mc, graph, sec, tstep)
            if not ys:
                continue
            got = True
            nu = len(ys) - 1 if (len(xs) > 1 and
                                 abs(xs[-1] - xs[0] - 360.0) < 1e-6) else len(ys)
            total += sum(ys[:nu])
        if not got:
            break
        series.append(total)
    return series


def read_2d(mc, graph, maxpts=64):
    ys = []
    for i in range(maxpts):
        try:
            _, y = mc.get_magnetic_graph_point(graph, i)
        except Exception:
            break
        ys.append(y)
    return ys


def main(argv):
    quit_after = "--quit" in argv
    os.makedirs(OUT_DIR, exist_ok=True)
    ts = time.strftime("%m%d_%H%M%S")

    from ansys.motorcad.core import MotorCAD
    print("启动 Motor-CAD (前台) ...")
    mc = MotorCAD()
    try:
        mc.set_visible(True)   # /SCRIPTING 模式部分机器窗口不显示, 强制可见
    except Exception:
        pass
    results = {"when": ts, "magnet_temp_C": MAGNET_TEMP_C,
               "reference": "V3.0 报告表3-1/3-3: 20°C, g=0.6/1.0/1.5 -> "
                            "601/483/378 N, kneg(1mm)=250 N/mm",
               "cases": []}
    try:
        mc.load_from_file(MOT_SRC)
        out_mot = os.path.join(OUT_DIR, "MARS_SSSR_compare_%s.mot" % ts)
        mc.save_to_file(out_mot)
        results["work_mot"] = out_mot

        t_before = mc.get_variable("Magnet_Temperature")
        results["magnet_temp_before_C"] = t_before
        mc.set_variable("Magnet_Temperature", MAGNET_TEMP_C)
        print("磁钢温度: %s -> %s °C" % (t_before, MAGNET_TEMP_C))
        for var in ["ElectromagneticForcesCalc_Load",
                    "ElectromagneticForcesCalc_OC"]:
            mc.set_variable(var, True)

        for g in GAPS:
            mc.set_variable("Airgap", g)
            back = mc.get_variable("Airgap")
            print("== 气隙 %.1f mm (回读 %s), 求解 ..." % (g, back))
            t0 = time.time()
            mc.do_magnetic_calculation()
            dt = time.time() - t0
            case = {"airgap_mm": back, "solve_seconds": dt}
            for graph, key in [("Fr_Rotor_OC_Lumped", "F_OC"),
                               ("Fr_Rotor_OL_Lumped", "F_OL")]:
                s = net_force_series(mc, graph)
                case[key] = stats(s)
            bys = read_2d(mc, "FluxDensityAirgap")
            if bys:
                case["B2_mean_T2"] = sum(b * b for b in bys) / len(bys)
            print("   F_OC=%.1f N, F_OL=%.1f N, mean(B²)=%.3f (耗时 %.0fs)"
                  % (case["F_OC"]["mean"], case["F_OL"]["mean"],
                     case.get("B2_mean_T2", -1), dt))
            results["cases"].append(case)

        # ---- 磁负刚度 (有限差分, OC 口径) ----
        cs = results["cases"]
        if len(cs) == 3:
            f = [c["F_OC"]["mean"] for c in cs]
            g0, g1, g2 = [c["airgap_mm"] for c in cs]
            k_low = -(f[1] - f[0]) / (g1 - g0)     # 0.6~1.0 段
            k_high = -(f[2] - f[1]) / (g2 - g1)    # 1.0~1.5 段
            k_mid = -(f[2] - f[0]) / (g2 - g0)     # 全段中心差分
            results["kneg_N_per_mm"] = {"seg_0.6_1.0": k_low,
                                        "seg_1.0_1.5": k_high,
                                        "central_at_1.0": k_mid}
            print("kneg: 0.6~1.0段 %.0f, 1.0~1.5段 %.0f, 中心差分 %.0f N/mm "
                  "(报告解析值 250)" % (k_low, k_high, k_mid))

        res_path = os.path.join(OUT_DIR, "compare_results_%s.json" % ts)
        with open(res_path, "w", encoding="utf-8") as fjson:
            json.dump(results, fjson, ensure_ascii=False, indent=2)
        print("RESULTS: %s" % res_path)
        return 0
    finally:
        if quit_after:
            try:
                mc.quit()
            except Exception:
                pass
        else:
            print("[提示] Motor-CAD 保持前台打开供检查。")


if __name__ == "__main__":
    sys.exit(main(sys.argv[1:]))