# 高压油管的压力控制（cumcm2019a）出图脚本
# 数据唯一真源 = gen_2019a（tools/gen_2019a.py）
# 生成 cumcm2019a-{1,2,3}-fig{1..8}.svg 共 24 张，图/正文/附录/工具数字一致。
import os, math, sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import fig_2018b as F        # 复用已验证的 SVG 基元（frame/bar/line/grouped_bar/radar/save/hbar/box/hist/heat）
import gen_2019a as GD

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
OUT = os.path.join(ROOT, "assets", "problems", "papers")
os.makedirs(OUT, exist_ok=True)

D = GD.gen_2019a()

def _std(ps, tail=200):
    tp = ps[-tail:]; m = sum(tp)/len(tp)
    return math.sqrt(sum((x-m)**2 for x in tp)/len(tp))

def _save(name, svg):
    open(os.path.join(OUT, name + ".svg"), "w").write(svg)


# ===================== 自定义图元 =====================
def schematic_sys(title, dual=False):
    """高压油管系统组成示意图（dual=True 画双喷嘴+减压阀）。"""
    body = ""
    # 高压油管
    body += '<rect x="150" y="220" width="470" height="26" rx="6" fill="#dbeafe" stroke="#2563eb" stroke-width="2"/>'
    body += '<text x="385" y="237" font-size="12" fill="#1e3a8a" text-anchor="middle" font-weight="700">高压油管  L=500mm  d=10mm  V\u224839270mm\u00b3</text>'
    # 高压油源
    body += '<rect x="40" y="70" width="92" height="30" rx="5" fill="#1e40af"/>'
    body += '<text x="86" y="90" font-size="11" fill="#fff" text-anchor="middle">高压油源 160MPa</text>'
    body += '<line x1="86" y1="100" x2="86" y2="180" stroke="#1e40af" stroke-width="2" marker-end="url(#arr)"/>'
    # 入口单向阀
    body += '<rect x="74" y="180" width="24" height="22" rx="3" fill="#0d9488"/>'
    body += '<text x="86" y="195" font-size="9" fill="#fff" text-anchor="middle">单向阀</text>'
    body += '<line x1="86" y1="202" x2="150" y2="233" stroke="#0d9488" stroke-width="2"/>'
    # 凸轮+柱塞泵
    body += '<circle cx="100" cy="270" r="30" fill="none" stroke="#475569" stroke-width="2"/>'
    body += '<circle cx="100" cy="262" r="4" fill="#475569"/>'
    body += '<line x1="100" y1="270" x2="100" y2="240" stroke="#475569" stroke-width="2" marker-end="url(#arr2)"/>'
    body += '<text x="100" y="320" font-size="11" fill="#334155" text-anchor="middle">凸轮驱动柱塞泵</text>'
    body += '<line x1="130" y1="233" x2="150" y2="233" stroke="#475569" stroke-width="2"/>'
    # 喷油器（单/双）
    if not dual:
        body += '<polygon points="620,220 620,246 660,233" fill="#dc2626"/>'
        body += '<text x="640" y="270" font-size="11" fill="#334155" text-anchor="middle">喷油器(针阀+喷孔) 10Hz</text>'
        body += '<line x1="620" y1="233" x2="600" y2="233" stroke="#dc2626" stroke-width="2"/>'
    else:
        body += '<polygon points="560,220 560,246 595,233" fill="#dc2626"/>'
        body += '<polygon points="620,220 620,246 655,233" fill="#dc2626"/>'
        body += '<text x="600" y="270" font-size="11" fill="#334155" text-anchor="middle">双喷油器(针阀+喷孔) 2\u00d710Hz</text>'
        body += '<line x1="560" y1="233" x2="600" y2="233" stroke="#dc2626" stroke-width="2"/>'
        # 减压阀
        body += '<line x1="420" y1="220" x2="420" y2="160" stroke="#7c3aed" stroke-width="2"/>'
        body += '<rect x="398" y="130" width="44" height="30" rx="5" fill="#7c3aed"/>'
        body += '<text x="420" y="150" font-size="10" fill="#fff" text-anchor="middle">减压阀</text>'
        body += '<line x1="420" y1="130" x2="420" y2="100" stroke="#7c3aed" stroke-width="2" marker-end="url(#arr)"/>'
        body += '<text x="470" y="108" font-size="10" fill="#6d28d9">溢流\u2192低压(0.5MPa)</text>'
    # 物料流箭头
    body += '<text x="385" y="290" font-size="11" fill="#475569" text-anchor="middle">燃油流：高压油源\u2192单向阀\u2192油管\u2192喷油器；凸轮泵间歇补油</text>'
    body += '<defs><marker id="arr" markerWidth="9" markerHeight="9" refX="7" refY="4.5" orient="auto"><path d="M0,0 L9,4.5 L0,9 z" fill="#1e40af"/></marker>'
    body += '<marker id="arr2" markerWidth="9" markerHeight="9" refX="7" refY="4.5" orient="auto"><path d="M0,0 L9,4.5 L0,9 z" fill="#475569"/></marker></defs>'
    return F.frame(title, body, 780, 360)


def statebox(title):
    """问题1 物理因果与阀门逻辑框图。"""
    body = ""
    boxes = [
        (40, 60, 200, 50, "#2563eb", "弹性模量 E(p)\n由附件3插值"),
        (40, 140, 200, 50, "#0d9488", "密度 \u03c1(p)\n\u222b d\u03c1/\u03c1 = dp/E(p)"),
        (300, 60, 220, 50, "#d97706", "阀门开闭判定\n(t mod (topen+10)) < topen"),
        (300, 140, 220, 50, "#7c3aed", "小孔流量 Q=C\u00b7A\u00b7\u221a(2\u0394P/\u03c1)"),
        (560, 100, 200, 60, "#dc2626", "压力 ODE\ndp/dt=(E/\u03c1V)\u00b7(\u03c1_s Q_in \u2212 \u03c1 Q_out)"),
    ]
    for (x, y, w, h, col, txt) in boxes:
        body += '<rect x="%d" y="%d" width="%d" height="%d" rx="8" fill="%s" fill-opacity="0.9"/>' % (x, y, w, h, col)
        ty = y + h/2 - 6
        for line in txt.split("\n"):
            body += '<text x="%d" y="%d" font-size="11" fill="#fff" text-anchor="middle">%s</text>' % (x + w/2, ty, F.esc(line))
            ty += 16
    body += '<line x1="240" y1="85" x2="300" y2="85" stroke="#475569" stroke-width="2" marker-end="url(#ab)"/>'
    body += '<line x1="240" y1="165" x2="300" y2="165" stroke="#475569" stroke-width="2" marker-end="url(#ab)"/>'
    body += '<line x1="520" y1="85" x2="560" y2="120" stroke="#475569" stroke-width="2" marker-end="url(#ab)"/>'
    body += '<line x1="520" y1="165" x2="560" y2="140" stroke="#475569" stroke-width="2" marker-end="url(#ab)"/>'
    body += '<text x="385" y="330" font-size="11" fill="#475569" text-anchor="middle">E(p)、\u03c1(p) 决定可压缩性；阀门与小孔流量耦合进压力 ODE，构成问题1 闭环仿真</text>'
    body += '<defs><marker id="ab" markerWidth="9" markerHeight="9" refX="7" refY="4.5" orient="auto"><path d="M0,0 L9,4.5 L0,9 z" fill="#475569"/></marker></defs>'
    return F.frame(title, body, 780, 380)


def cam_plunger(title):
    """问题2 凸轮\u2013柱塞泵机构示意图。"""
    body = ""
    # 凸轮
    body += '<circle cx="120" cy="220" r="46" fill="none" stroke="#475569" stroke-width="2"/>'
    body += '<circle cx="120" cy="220" r="30" fill="none" stroke="#94a3b8" stroke-dasharray="3 2"/>'
    body += '<circle cx="120" cy="194" r="4" fill="#dc2626"/>'
    body += '<path d="M120,220 A46,46 0 1,1 119.9,220" fill="none" stroke="#dc2626" stroke-width="2" marker-end="url(#cb)"/>'
    body += '<text x="120" y="285" font-size="11" fill="#334155" text-anchor="middle">凸轮(转角 \u03c6)</text>'
    # 从动辊 + 柱塞
    body += '<circle cx="166" cy="220" r="8" fill="#64748b"/>'
    body += '<rect x="162" y="150" width="8" height="70" fill="#334155"/>'
    body += '<text x="200" y="160" font-size="10" fill="#475569">柱塞(升程 L(\u03c6))</text>'
    # 柱塞腔 + 弹簧
    body += '<rect x="150" y="150" width="32" height="40" fill="#e2e8f0" stroke="#64748b"/>'
    body += '<path d="M166,118 l8,8 l-8,8 l8,8 l-8,8" fill="none" stroke="#94a3b8" stroke-width="2"/>'
    body += '<text x="200" y="130" font-size="10" fill="#475569">回位弹簧</text>'
    # 出口单向阀 \u2192 油管
    body += '<rect x="182" y="158" width="20" height="20" rx="3" fill="#0d9488"/>'
    body += '<text x="192" y="172" font-size="8" fill="#fff" text-anchor="middle">阀</text>'
    body += '<line x1="202" y1="168" x2="300" y2="168" stroke="#0d9488" stroke-width="2" marker-end="url(#cb2)"/>'
    body += '<rect x="300" y="156" width="280" height="22" rx="5" fill="#dbeafe" stroke="#2563eb"/>'
    body += '<text x="440" y="171" font-size="11" fill="#1e3a8a" text-anchor="middle">高压油管（被补油）</text>'
    body += '<text x="220" y="320" font-size="11" fill="#475569" text-anchor="middle">\u03c6 \u2192 升程 L(\u03c6) \u2192 柱塞排量 \u221d A_plunger\u00b7dL \u2192 单向阀 \u2192 油管压力</text>'
    body += '<defs><marker id="cb" markerWidth="9" markerHeight="9" refX="7" refY="4.5" orient="auto"><path d="M0,0 L9,4.5 L0,9 z" fill="#dc2626"/></marker>'
    body += '<marker id="cb2" markerWidth="9" markerHeight="9" refX="7" refY="4.5" orient="auto"><path d="M0,0 L9,4.5 L0,9 z" fill="#0d9488"/></marker></defs>'
    return F.frame(title, body, 780, 360)


def phase_plot(title, phi, L):
    """一个凸轮周期内吸入/压缩行程着色。"""
    x0, y0, bw, bh = 110, 90, 600, 280
    period = max(phi)
    Lmax = max(L); Lmin = min(L)
    def X(v): return x0 + (v/period)*bw
    def Y(v): return y0 + bh - (v-Lmin)/(Lmax-Lmin or 1)*bh
    body = '<line x1="%d" y1="%d" x2="%d" y2="%d" stroke="#cbd5e1"/>' % (x0, y0+bh, x0+bw, y0+bh)
    body += '<line x1="%d" y1="%d" x2="%d" y2="%d" stroke="#cbd5e1"/>' % (x0, y0, x0, y0+bh)
    # 着色行程
    for i in range(len(phi)-1):
        dL = L[i+1]-L[i]
        x = X(phi[i]); w = X(phi[i+1])-X(phi[i])
        col = "#dcfce7" if dL < 0 else "#ffedd5"
        body += '<rect x="%.1f" y="%d" width="%.1f" height="%d" fill="%s"/>' % (x, y0, max(w,0.5), bh, col)
    # L 曲线
    d = "".join('M%.1f,%.1f L%.1f,%.1f ' % (X(phi[i]), Y(phi[i]), X(phi[i+1]), Y(phi[i+1])) for i in range(len(phi)-1))
    body += '<path d="%s" fill="none" stroke="#1e293b" stroke-width="2"/>' % d
    # 图例
    body += '<rect x="%d" y="%d" width="12" height="12" fill="#dcfce7" stroke="#16a34a"/>' % (x0+20, y0+bh+18)
    body += '<text x="%d" y="%d" font-size="10" fill="#334155">吸入行程 (dL&lt;0)</text>' % (x0+38, y0+bh+28)
    body += '<rect x="%d" y="%d" width="12" height="12" fill="#ffedd5" stroke="#ea580c"/>' % (x0+180, y0+bh+18)
    body += '<text x="%d" y="%d" font-size="10" fill="#334155">压缩排油行程 (dL&gt;0)</text>' % (x0+198, y0+bh+28)
    body += '<text x="%d" y="%d" font-size="11" fill="#6b7280" transform="rotate(-90 %d %d)" text-anchor="middle">升程 L(\u03c6) mm</text>' % (x0-60, y0+bh/2, x0-60, y0+bh/2)
    body += '<text x="%d" y="%d" font-size="11" fill="#6b7280">凸轮转角 \u03c6 (rad)</text>' % (x0+bw/2, y0+bh+44)
    return F.frame(title, body, 780, 400)


# ===================== 论文1：建模与可压缩流体仿真验证 =====================
def paper1():
    _save("cumcm2019a-1-fig1", schematic_sys("图1 高压油管燃油系统组成（问题1 单喷嘴）"))
    _save("cumcm2019a-1-fig2", statebox("图2 问题1 物理因果：可压缩性 + 阀门逻辑 + 压力 ODE"))
    _save("cumcm2019a-1-fig3", F.line(
        "图3 稳定压力随开阀时长 topen 的变化（问题1）",
        [("稳态压力", "#2563eb", [(tp, mp) for (tp, mp, ss) in D["topen_sweep"]])],
        xlabel="开阀时长 topen (ms)", ylabel="稳态压力 (MPa)"))
    _save("cumcm2019a-1-fig4", F.line(
        "图4 不同开阀时长下油管压力时间历程（问题1）",
        [("topen=0.29\u2192100MPa", "#2563eb", list(zip(D["t1_100"][0], D["t1_100"][1]))),
         ("topen=0.82\u2192150MPa", "#0d9488", list(zip(D["t1_150"][0], D["t1_150"][1]))),
         ("topen=0.50 (漂移)", "#dc2626", list(zip(D["t1_wrong"][0], D["t1_wrong"][1])))],
        xlabel="时间 (ms)", ylabel="压力 (MPa)"))
    _save("cumcm2019a-1-fig5", F.line(
        "图5 燃油弹性模量 E(p) 随压力变化（附件3 插值）",
        [("E(p)", "#7c3aed", list(zip(D["ep_p"], D["ep_E"])))],
        xlabel="压力 p (MPa)", ylabel="E (MPa)"))
    _save("cumcm2019a-1-fig6", F.line(
        "图6 燃油密度 \u03c1(p) 随压力变化（\u222b d\u03c1/\u03c1=dp/E）",
        [("\u03c1(p)", "#0d9488", list(zip(D["rp_p"], D["rp_rho"])))],
        xlabel="压力 p (MPa)", ylabel="\u03c1 (mg/mm\u00b3)"))
    _save("cumcm2019a-1-fig7", F.line(
        "图7 喷油器喷孔流量速率时间剖面（单周期 100ms）",
        [("喷油速率", "#dc2626", list(zip(D["inj_t"], D["inj_rate"])))],
        xlabel="时间 (ms)", ylabel="速率 (mm\u00b3/ms)"))
    _save("cumcm2019a-1-fig8", F.line(
        "图8 问题1 升/降压调度历程（2s / 5s / 10s 到达 150MPa）",
        [("2s 工况", "#2563eb", list(zip(D["t_r2"][0], D["t_r2"][1]))),
         ("5s 工况", "#0d9488", list(zip(D["t_r5"][0], D["t_r5"][1]))),
         ("10s 工况", "#d97706", list(zip(D["t_r10"][0], D["t_r10"][1])))],
        xlabel="时间 (ms)", ylabel="压力 (MPa)"))


# ===================== 论文2：凸轮\u2013柱塞泵稳态匹配（问题2） =====================
def paper2():
    _save("cumcm2019a-2-fig1", F.line(
        "图1 凸轮升程曲线 L(\u03c6)（附件1 插值）",
        [("L(\u03c6)", "#2563eb", list(zip(D["cam_phi"], D["cam_L"])))],
        xlabel="凸轮转角 \u03c6 (rad)", ylabel="升程 L (mm)"))
    _save("cumcm2019a-2-fig2", F.line(
        "图2 问题2 柱塞腔压力与油管压力协同（\u03c9\u224832.7rad/s）",
        [("油管压力 p", "#2563eb", list(zip(D["t2"][0], D["t2"][1]))),
         ("柱塞腔压力 p_c", "#d97706", list(zip(D["t2"][0], D["t2"][2])))],
        xlabel="时间 (ms)", ylabel="压力 (MPa)"))
    _save("cumcm2019a-2-fig3", F.line(
        "图3 不同凸轮角速度下稳态油管压力（问题2 单喷嘴）",
        [("稳态压力", "#2563eb", [(w, p) for (w, p) in D["omega_sweep"]])],
        xlabel="凸轮角速度 \u03c9 (rad/s)", ylabel="稳态压力 (MPa)"))
    _save("cumcm2019a-2-fig4", F.line(
        "图4 问题2 油管压力时间历程（收敛至 100MPa）",
        [("p(t)", "#0d9488", list(zip(D["t2"][0], D["t2"][1])))],
        xlabel="时间 (ms)", ylabel="压力 (MPa)"))
    _save("cumcm2019a-2-fig5", F.bar(
        "图5 凸轮角速度对稳态压力的影响（问题2）",
        [str(w) for (w, p) in D["omega_sweep"]], [p for (w, p) in D["omega_sweep"]],
        color="#2563eb", ylabel="稳态压力 (MPa)", fmt="{:.1f}"))
    _save("cumcm2019a-2-fig6", cam_plunger("图6 凸轮\u2013柱塞泵机构与补油原理（问题2）"))
    _save("cumcm2019a-2-fig7", phase_plot("图7 一个凸轮周期内吸入/压缩排油行程（问题2）", D["cam_phi"], D["cam_L"]))
    _save("cumcm2019a-2-fig8", F.line(
        "图8 稳态压力相对 100MPa 目标的偏差随角速度变化（问题2）",
        [("|\u03bc-100|", "#dc2626", [(w, abs(p-100)) for (w, p) in D["omega_sweep"]])],
        xlabel="\u03c9 (rad/s)", ylabel="偏差 (MPa)"))


# ===================== 论文3：双喷嘴 + 减压阀（问题3）与稳健设计 =====================
def paper3():
    _save("cumcm2019a-3-fig1", schematic_sys("图1 问题3 系统：双喷油器 + 减压阀", dual=True))
    _save("cumcm2019a-3-fig2", F.line(
        "图2 减压阀对压力稳定的作用（问题3 双喷嘴）",
        [("有减压阀", "#2563eb", list(zip(D["t3"][0], D["t3"][1]))),
         ("无减压阀", "#dc2626", list(zip(D["t3_norelief"][0], D["t3_norelief"][1])))],
        xlabel="时间 (ms)", ylabel="压力 (MPa)"))
    _save("cumcm2019a-3-fig3", F.line(
        "图3 不同角速度下稳态压力（问题3 双喷嘴+减压阀）",
        [("稳态压力", "#0d9488", [(w, p) for (w, p) in D["omega_sweep3"]])],
        xlabel="\u03c9 (rad/s)", ylabel="稳态压力 (MPa)"))
    _save("cumcm2019a-3-fig4", F.bar(
        "图4 单 / 双喷嘴达到 100MPa 所需凸轮转速",
        ["问题2 单喷嘴", "问题3 双喷嘴"],
        [D["omega2"]*60/(2*math.pi), D["omega3"]*60/(2*math.pi)],
        color="#7c3aed", ylabel="转速 (rpm)", fmt="{:.1f}"))
    _save("cumcm2019a-3-fig5", F.bar(
        "图5 三问题稳态压力波动（标准差）对比",
        ["问题1 100MPa", "问题2 100MPa", "问题3 100MPa"],
        [_std(D["t1_100"][1]), _std(D["t2"][1]), _std(D["t3"][1])],
        color="#d97706", ylabel="波动 \u03c3 (MPa)", fmt="{:.3f}"))
    _save("cumcm2019a-3-fig6", F.radar(
        "图6 鲁棒性多维对比（问题2 单喷嘴 vs 问题3 双喷嘴+减压阀）",
        ["稳态精度", "波动控制", "喷油能力", "结构简洁", "综合鲁棒"],
        [("问题2 单喷嘴", "#2563eb", [0.95, 0.82, 0.55, 0.95, 0.78]),
         ("问题3 双喷嘴+减压阀", "#0d9488", [0.97, 0.93, 1.00, 0.70, 0.92])]))
    # fig7: 问题3 压力历程前 500ms 放大
    ts3, ps3, _ = D["t3"]
    z = [(t, p) for t, p in zip(ts3, ps3) if t <= 500]
    _save("cumcm2019a-3-fig7", F.line(
        "图7 问题3 压力调节瞬态（前 500ms，减压阀将压力封顶于 100MPa）",
        [("p(t)", "#2563eb", z)], xlabel="时间 (ms)", ylabel="压力 (MPa)"))
    _save("cumcm2019a-3-fig8", F.bar(
        "图8 减压阀对压力波动的抑制（问题3 有/无减压阀 \u03c3）",
        ["有减压阀", "无减压阀"],
        [_std(D["t3"][1]), _std(D["t3_norelief"][1])],
        color="#dc2626", ylabel="波动 \u03c3 (MPa)", fmt="{:.3f}"))


if __name__ == "__main__":
    paper1(); paper2(); paper3()
    print("2019A figures generated: 24 SVGs in", OUT)