#!/usr/bin/env python3
"""
Clarke 变换验证脚本

验证目标：
1. 对称三相电流 → Iα、Iβ 正交
2. 正变换 → 反变换 → 还原
3. Q15 定点数精度对比
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')  # 无头模式
import matplotlib.pyplot as plt
import os

SQRT3_INV = 1.0 / np.sqrt(3)
SQRT3_OVER2 = np.sqrt(3) / 2.0

def clarke_transform(Ia, Ib):
    """Clarke 变换浮点版"""
    Ialpha = Ia
    Ibeta = (Ia + 2.0 * Ib) * SQRT3_INV
    return Ialpha, Ibeta

def clarke_inv_transform(Ialpha, Ibeta):
    """Clarke 反变换浮点版"""
    Ia = Ialpha
    Ib = -0.5 * Ialpha + SQRT3_OVER2 * Ibeta
    Ic = -(Ia + Ib)
    return Ia, Ib, Ic

# ───────── 测试 1：对称三相电流验证 ─────────
print("=" * 60)
print("测试 1：对称三相电流（Ia=1, Ib=-0.5, Ic=-0.5）")
print("=" * 60)

Ia, Ib = 1.0, -0.5
Ialpha, Ibeta = clarke_transform(Ia, Ib)
print(f"  输入: Ia={Ia:.3f}, Ib={Ib:.3f}, Ic={-(Ia+Ib):.3f}")
print(f"  输出: Iα={Ialpha:.6f}, Iβ={Ibeta:.6f}")
print(f"  预期: Iα=1.000000, Iβ=0.000000")
print(f"  结论: {'✓ 正确' if abs(Ialpha-1.0)<0.001 and abs(Ibeta)<0.001 else '✗ 错误'}")

# ───────── 测试 2：正反变换一致性 ─────────
print()
print("=" * 60)
print("测试 2：正反变换一致性")
print("=" * 60)

np.random.seed(42)
max_err = 0
for _ in range(1000):
    Ia_orig = np.random.uniform(-1, 1)
    Ib_orig = np.random.uniform(-1, 1)
    # 确保 Ic 在合理范围内 (Ic = -Ia - Ib)
    Ic_orig = -(Ia_orig + Ib_orig)
    
    Ialpha, Ibeta = clarke_transform(Ia_orig, Ib_orig)
    Ia_back, Ib_back, Ic_back = clarke_inv_transform(Ialpha, Ibeta)
    
    err = max(abs(Ia_orig - Ia_back), abs(Ib_orig - Ib_back))
    max_err = max(max_err, err)

print(f"  最大误差: {max_err:.2e}")
print(f"  结论: {'✓ 通过' if max_err < 1e-6 else '✗ 异常'}")

# ───────── 测试 3：三相正弦信号 → 正交输出 ─────────
print()
print("=" * 60)
print("测试 3：三相正弦 → Clarke → 正交 Iα/Iβ")
print("=" * 60)

theta = np.linspace(0, 2 * np.pi, 200)
Ia_sin = np.sin(theta)
Ib_sin = np.sin(theta - 2 * np.pi / 3)
Ic_sin = np.sin(theta + 2 * np.pi / 3)

Ia_c, Ib_c = Ia_sin, Ib_sin  # Ia=sin(θ), Ib=sin(θ-120°)
Ialpha_arr, Ibeta_arr = clarke_transform(Ia_c, Ib_c)

# 验证正交性：Iα, Iβ 的相关系数应接近 0
corr = np.corrcoef(Ialpha_arr, Ibeta_arr)[0, 1]
print(f"  Iα 与 Iβ 的相关系数: {corr:.6f} (应为 0)")
print(f"  结论: {'✓ 正交' if abs(corr) < 0.01 else '✗ 不正交'}")

# 验证幅值相等
amp_alpha = np.max(np.abs(Ialpha_arr))
amp_beta  = np.max(np.abs(Ibeta_arr))
print(f"  Iα 幅值: {amp_alpha:.6f}")
print(f"  Iβ 幅值: {amp_beta:.6f}")
print(f"  结论: {'✓ 幅值一致' if abs(amp_alpha - amp_beta) < 0.001 else '✗ 幅值不等'}")

# 生成验证图
fig, axes = plt.subplots(2, 2, figsize=(12, 8))

# 图1：三相电流
axes[0, 0].plot(theta, Ia_sin, label='Ia')
axes[0, 0].plot(theta, Ib_sin, label='Ib')
axes[0, 0].plot(theta, Ic_sin, label='Ic')
axes[0, 0].set_title('三相电流 (ABC)')
axes[0, 0].set_xlabel('电角度 (rad)')
axes[0, 0].legend()
axes[0, 0].grid(True)

# 图2：αβ 电流
axes[0, 1].plot(theta, Ialpha_arr, label='Iα')
axes[0, 1].plot(theta, Ibeta_arr, label='Iβ')
axes[0, 1].set_title('Clarke 变换后 (αβ)')
axes[0, 1].set_xlabel('电角度 (rad)')
axes[0, 1].legend()
axes[0, 1].grid(True)

# 图3：相量图
axes[1, 0].plot(Ialpha_arr, Ibeta_arr)
axes[1, 0].set_title('αβ 相量图')
axes[1, 0].set_xlabel('Iα')
axes[1, 0].set_ylabel('Iβ')
axes[1, 0].axis('equal')
axes[1, 0].grid(True)

# 图4：Ia + Ib + Ic = 0 验证
sum_abc = Ia_sin + Ib_sin + Ic_sin
axes[1, 1].plot(theta, sum_abc)
axes[1, 1].set_title('Ia + Ib + Ic (应为 0)')
axes[1, 1].set_xlabel('电角度 (rad)')
axes[1, 1].set_ylabel('和')
axes[1, 1].grid(True)
axes[1, 1].set_ylim(-1e-14, 1e-14)  # 浮点误差级别

plt.tight_layout()

# 保存图
out_dir = os.path.dirname(os.path.abspath(__file__))
plt.savefig(os.path.join(out_dir, 'clarke_verify.png'), dpi=150)
print(f"\n  验证图已保存: {os.path.join(out_dir, 'clarke_verify.png')}")

# ───────── 测试 4：Q15 定点数精度对比 ─────────
print()
print("=" * 60)
print("测试 4：float vs Q15 定点数精度对比")
print("=" * 60)

SQRT3_INV_Q15 = 18919
Q15_FACTOR = 32768.0

theta_test = np.linspace(0, 2 * np.pi, 100)
q15_errors = []

for t in theta_test:
    Ia_f = np.sin(t)
    Ib_f = np.sin(t - 2 * np.pi / 3)
    
    # float 结果
    Ia_f_result, Ibeta_f_result = clarke_transform(Ia_f, Ib_f)
    
    # Q15 模拟
    Ia_q15 = int(round(Ia_f * Q15_FACTOR))
    Ib_q15 = int(round(Ib_f * Q15_FACTOR))
    
    tmp = Ia_q15 + 2 * Ib_q15
    tmp = tmp * SQRT3_INV_Q15
    Ibeta_q15 = tmp >> 15
    
    # Q15 转回 float 比较
    Ibeta_q15_f = Ibeta_q15 / Q15_FACTOR
    
    q15_errors.append(abs(Ibeta_q15_f - Ibeta_f_result))

print(f"  Q15 Iβ 最大误差: {max(q15_errors):.2e}")
print(f"  Q15 Iβ 平均误差: {np.mean(q15_errors):.2e}")
print(f"  Q15 Iβ 相对精度: ~{1/np.mean(q15_errors):.0f} 分之一")

print()
print("=" * 60)
print("所有测试完成 ✓")
print("=" * 60)