#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""PD132T vs PD136HA 响应对比图（纯 stdlib 手写 SVG）

面板 A：同一判据下的 CAPVD 响应（时域；Δf/f = +0.5% 大信号 与 +0.15% 近阈值两组）
面板 B：检测延迟 vs Δf/f（对数扫描）—— 看「窗主导」与「τ 主导」的分界

转 PNG：chromium --headless=new --disable-gpu --no-sandbox \
        --screenshot=docs/img/response-compare.svg.png --window-size=1520,980 \
        file://$PWD/docs/img/response-compare-pd132t-pd136ha.svg
"""
import math

FONT = "'Noto Serif CJK SC','Source Han Serif SC','WenQuanYi Zen Hei',sans-serif"
W, H = 1520.0, 980.0
o = []
def add(s): o.append(s)
def esc(s): return s.replace('&', '&amp;').replace('<', '&lt;').replace('>', '&gt;')
def txt(x, y, s, fs=14, fill='#222', anchor='start', w='400'):
    add(f'<text x="{x:.1f}" y="{y:.1f}" font-family="{FONT}" font-size="{fs}" '
        f'fill="{fill}" text-anchor="{anchor}" font-weight="{w}">{esc(s)}</text>')
def line(x1, y1, x2, y2, col='#888', w=1, dash=None):
    d = f' stroke-dasharray="{dash}"' if dash else ''
    add(f'<line x1="{x1:.1f}" y1="{y1:.1f}" x2="{x2:.1f}" y2="{y2:.1f}" stroke="{col}" '
        f'stroke-width="{w}"{d}/>')
def rect(x, y, w_, h_, fill='none', stroke='#999', sw=1, rx=4, op=1.0):
    add(f'<rect x="{x:.1f}" y="{y:.1f}" width="{w_:.1f}" height="{h_:.1f}" rx="{rx}" '
        f'fill="{fill}" fill-opacity="{op}" stroke="{stroke}" stroke-width="{sw}"/>')
def poly(pts, col, w=2.2, dash=None):
    d = f' stroke-dasharray="{dash}"' if dash else ''
    p = ' '.join(f'{x:.1f},{y:.1f}' for x, y in pts)
    add(f'<polyline points="{p}" fill="none" stroke="{col}" stroke-width="{w}"{d}/>')
def dot(x, y, r=4.2, col='#c1121f'):
    add(f'<circle cx="{x:.1f}" cy="{y:.1f}" r="{r}" fill="{col}"/>')
def note(x, y, w_, h_, title, lines_, col='#0b6'):
    rect(x, y, w_, h_, fill='#f7fbff', stroke=col)
    txt(x + 12, y + 26, title, 15, '#0b4f7a', w='700')
    for i, s in enumerate(lines_):
        txt(x + 12, y + 52 + i * 21, s, 14, '#333')

# ============ 模型参数（与 measurement-window-quantization.md 同源公式） ============
def model(name, f_cap, edge, alfa, sens_in, color, f_coil=100e3, base=32768, div=32):
    xn0 = round(div * f_cap / (edge * f_coil))
    lp = base // xn0
    tw = lp * xn0 / f_cap
    return dict(name=name, color=color, alfa=alfa, thr=sens_in / 65536.0,
                tw=tw, xn0=xn0, lp=lp, tau=-tw / math.log(1 - alfa / 256.0),
                tau_pl=tw * 256.0 / alfa)

PT = model('PD132T',  3.0e6, 2, 79, 60, '#1565c0')
PH = model('PD136HA', 1.5e6, 1, 64, 81, '#e65100')

def resp(m, df, n):
    """第 n 次（0-based）CAPVD 更新后的 Δf/f（相对值）"""
    return df * (1.0 - (1.0 - m['alfa'] / 256.0) ** (n + 1))

def delay_ms(m, df, cap=20000):
    """返回 (窗数-截止更新序号 idx, 时间 ms)；idx=None 表示 cap 内不可检"""
    n = 0
    while n < cap:
        if resp(m, df, n) >= m['thr']:
            return n, (n + 1) * m['tw'] * 1e3
        n += 1
    return None, None

# ================== 面板 A：时域响应 ==================
AX0, AX1, AY0, AY1 = 110.0, 700.0, 640.0, 150.0
TMAX = 300.0          # ms
VMAX = 0.62           # %
ax = lambda t: AX0 + (t / TMAX) * (AX1 - AX0)     # t in ms
ay = lambda v: AY0 - (v / VMAX) * (AY0 - AY1)     # v in %

txt(30, 40, 'PD132T ↔ PD136HA 响应对比：同一判据口径，不同测量节拍与滤波', 21, '#111', w='700')
txt(30, 66, 'PD132T: 窗 10.880 ms / α=79 / τ=29.5 ms（产品线口径 35.3 ms）    '
            'PD136HA: 窗 21.760 ms / α=64 / τ=76.0 ms（产品线口径 87.0 ms）    站点 100 kHz', 13.5, '#555')

txt(AX0 - 10, AY1 - 18, '面板 A：CAPVD 响应（时域，一次完全占据）', 15.5, '#0b4f7a', 'start', '700')
rect(AX0, AY1, AX1 - AX0, AY0 - AY1, fill='#fcfcfc', stroke='#ccc')
base_l = AY0 + 38
line(AX0, AY0, AX1, AY0, '#888', 1.5)
line(AX0, AY1, AX0, AY0, '#888', 1.5)
for i in range(7):
    v = i * 0.1
    line(AX0, ay(v), AX1, ay(v), '#eee', 1)
    txt(AX0 - 8, ay(v) + 5, '%.1f%%' % v, 12.5, '#666', 'end')
txt(AX0 + 10, AY1 + 22, 'Δf/f', 13, '#444', w='700')
for t in range(0, 301, 50):
    line(ax(t), AY0, ax(t), AY0 + 6, '#888', 1)
    txt(ax(t), base_l, str(t), 12.5, '#666', 'middle')
txt((AX0 + AX1) / 2, base_l + 24, '时间（ms）', 13.5, '#444', 'middle', '700')

# 阈值线（各型号档1 进入）
line(AX0, ay(PT['thr'] * 100), AX1, ay(PT['thr'] * 100), '#1565c0', 1.4, '6,4')
line(AX0, ay(PH['thr'] * 100), AX1, ay(PH['thr'] * 100), '#e65100', 1.4, '6,4')
txt(ax(205), ay(PT['thr'] * 100) + 20, 'PD132T 档1 0.0916%', 12.5, '#1565c0')
txt(ax(205), ay(PH['thr'] * 100) + 20, 'PD136HA 档1 0.1236%', 12.5, '#e65100')

# 两条响应（阶梯）
def stair_ms(m, df, tmax=TMAX):
    pts = [(ax(0.0), ay(0.0))]
    t = 0.0
    n = 0
    while t <= tmax:
        pts.append((ax(t), ay(resp(m, df * 100, n) if False else resp(m, df, n) * 100)))
        t += m['tw'] * 1e3
        n += 1
    return pts

for m in (PT, PH):
    for df, dash, lw in ((0.005, None, 2.6), (0.0015, '7,5', 1.8)):
        pts = [(ax(0.0), ay(0.0))]
        t = 0.0
        n = 0
        while t <= TMAX:
            pts.append((ax(t), ay(resp(m, df, n) * 100)))
            t += m['tw'] * 1e3
            n += 1
        poly(pts, m['color'], lw, dash)
# 越线标记
for m, df in ((PT, 0.005), (PH, 0.005), (PT, 0.0015), (PH, 0.0015)):
    i, ms = delay_ms(m, df)
    if i is None:
        continue
    dot(ax(ms), ay(resp(m, df, i) * 100), 4.6 if df == 0.005 else 3.4, m['color'])

iT5, msT5 = delay_ms(PT, 0.005); iH5, msH5 = delay_ms(PH, 0.005)
iT15, msT15 = delay_ms(PT, 0.0015); iH15, msH15 = delay_ms(PH, 0.0015)
note(110, 716, 580, 158, '越线时刻与读图要点（同一判据口径）', [
    '大信号 +0.5%%：PD132T 第 %d 窗 ⇒ %.1f ms ｜ PD136HA 第 %d 窗 ⇒ %.1f ms' % (iT5 + 1, msT5, iH5 + 1, msH5),
    '近阈值 +0.15%%：PD132T 第 %d 窗 ⇒ %.1f ms ｜ PD136HA 第 %d 窗 ⇒ %.1f ms' % (iT15 + 1, msT15, iH15 + 1, msH15),
    '⇒ 大信号由「窗长」决定（2×）；近阈值由「τ」决定（≈4.7×）',
    '面板 B：右侧平坦段 = 1 窗下限；左侧陡升 = IIR(τ) 主导；阈值以下永不检到',
], '#0b6')
txt(110, 890, '蓝 = PD132T，橙 = PD136HA；实线 = 大信号（+0.5%），虚线 = 近阈值（+0.15%）', 13, '#555')

# ================== 面板 B：检测延迟 vs Δf/f（对数扫描） ==================
BX0, BX1, BYB, BYT = 830.0, 1430.0, 830.0, 170.0
DLO, DHI = 4.0, 3000.0      # ms
RLO, RHI = 0.02, 5.0        # %
bx = lambda v: BX0 + (math.log10(v / RLO) / math.log10(RHI / RLO)) * (BX1 - BX0)
by = lambda d: BYB - (math.log10(d / DLO) / math.log10(DHI / DLO)) * (BYB - BYT)

txt(BX0 - 10, BYT - 18, '面板 B：检测延迟随信号强度（对数扫描）', 15.5, '#0b4f7a', 'start', '700')
rect(BX0, BYT, BX1 - BX0, BYB - BYT, fill='#fcfcfc', stroke='#ccc')
line(BX0, BYB, BX1, BYB, '#888', 1.5)
line(BX0, BYT, BX0, BYB, '#888', 1.5)
for d in (4, 10, 30, 100, 300, 1000, 3000):
    line(BX0, by(d), BX1, by(d), '#eee', 1)
    txt(BX0 - 8, by(d) + 5, ('%g' % d) + ' ms', 12.5, '#666', 'end')
for v in (0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5):
    line(bx(v), BYB, bx(v), BYB - 6, '#888', 1)
    txt(bx(v), BYB + 22, ('%g' % v) + '%', 12.5, '#666', 'middle')
txt((BX0 + BX1) / 2, BYB + 48, '车辆引起的 Δf/f（%，对数轴）', 13.5, '#444', 'middle', '700')

# 典型车型区间
rect(bx(0.1), BYT, bx(1.0) - bx(0.1), BYB - BYT, fill='#f2c200', stroke='none', op=0.13)
txt(bx(0.1) + 6, BYT + 20, '典型车型 0.1%~1%', 12.5, '#8a6d00')

# 曲线
for m in (PT, PH):
    pts = []
    v = RLO * 1.05
    while v <= RHI:
        i, ms = delay_ms(m, v / 100.0, cap=4000)
        v *= 1.02
        if i is None:
            continue
        if ms > DHI:
            break
        pts.append((bx(v), by(ms)))
    poly(pts, m['color'], 2.6)
    # 阈值竖线
    thrp = m['thr'] * 100
    line(bx(thrp), BYB, bx(thrp), BYT, m['color'], 1.4, '6,4')
    anc = 'end' if m is PT else 'start'
    off = -6 if m is PT else 6
    txt(bx(thrp) + off, BYT - 6, '%s 档1 %.4f%%' % (m['name'], thrp), 12.5, m['color'], anc)
    # 1 窗平台标注
    line(BX0, by(m['tw'] * 1e3), bx(thrp), by(m['tw'] * 1e3), m['color'], 1.2, '3,4')
    txt(bx(thrp) + 8, by(m['tw'] * 1e3) - 7, '≈1 窗 %.2f ms' % (m['tw'] * 1e3), 12.5, m['color'])
# 不可检区域（阈值左侧）
for m in (PT, PH):
    rect(BX0, BYT, bx(m['thr'] * 100) - BX0, BYB - BYT, fill='#bbb', stroke='none', op=0.18)
txt(BX0 + 8, BYB - 12, '← 阈值以下：永远检不到（灰区）', 12.5, '#555')

dT_big = iT5 and msT5; dH_big = msH5


# ================== 页脚 ==================
note(30, 912, 700, 62, '公式（与 measurement-window-quantization.md §1~§4 同源）', [
    'CAPVD_n = Δf/f·[1 − (1−α/256)^n]（n = 第 n 次窗更新）；延迟 = 首个满足 CAPVD_n ≥ 判据 的 (n+1)·T_win',
], '#0b6')
note(750, 912, 740, 62, '注', [
    '未计输出级 50 ms 主节拍量化（实际继电器动作再加 ≤50 ms）；两型号均为单次越线判定（无进入/离开确认）。',
], '#0b6')

add('</svg>')
svg = ('<svg xmlns="http://www.w3.org/2000/svg" width="%d" height="%d" viewBox="0 0 %d %d">'
       '<rect width="%d" height="%d" fill="#fff"/>' % (W, H, W, H, W, H)) + ''.join(o)
out = 'docs/img/response-compare-pd132t-pd136ha.svg'
open(out, 'w', encoding='utf-8').write(svg)
print('SVG:', out, len(svg), 'bytes')
print('PD132T : tw=%.3f ms tau=%.1f ms(口径%.1f) 档1=%.4f%%' % (PT['tw'] * 1e3, PT['tau'] * 1e3, PT['tau_pl'] * 1e3, PT['thr'] * 100))
print('PD136HA: tw=%.3f ms tau=%.1f ms(口径%.1f) 档1=%.4f%%' % (PH['tw'] * 1e3, PH['tau'] * 1e3, PH['tau_pl'] * 1e3, PH['thr'] * 100))
print('+0.50%%: PD132T %.1f ms(第%d窗) | PD136HA %.1f ms(第%d窗)' % (msT5, iT5 + 1, msH5, iH5 + 1))
print('+0.15%%: PD132T %.1f ms(第%d窗) | PD136HA %.1f ms(第%d窗)' % (msT15, iT15 + 1, msH15, iH15 + 1))
print('+0.10%%: PD132T %.1f ms | PD136HA %.1f ms' % (delay_ms(PT, 0.001)[1] or -1, delay_ms(PH, 0.001)[1] or -1))