import logging
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')  # headless backend
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.signal import savgol_filter as _savgol_filter

logger = logging.getLogger(__name__)


def load_data(file_path_or_buffer):
    """加载CSV文件中的单摆数据 (path or file-like)"""
    try:
        data = pd.read_csv(file_path_or_buffer)
        return data
    except Exception as e:
        logger.exception(f"读取数据文件时出错: {e}")
        return None


def fit_circle(x, y):
    """
    通过轨迹点拟合圆，找出圆心和半径
    返回: x_c, y_c, r
    """
    A = np.column_stack([x, y, np.ones(len(x))])
    b = x**2 + y**2
    try:
        c = np.linalg.lstsq(A, b, rcond=None)[0]
        x_c = c[0] / 2
        y_c = c[1] / 2
        r = np.sqrt(c[2] + x_c**2 + y_c**2)
        return x_c, y_c, r
    except np.linalg.LinAlgError:
        logger.warning("拟合圆时出现线性代数错误，使用备用方法")
        x_c = np.mean(x)
        y_c = np.min(y) - np.sqrt(np.var(x))
        distances = np.sqrt((x - x_c)**2 + (y - y_c)**2)
        r = np.mean(distances)
        return x_c, y_c, r


def calculate_angle_from_pendulum_length(data, pendulum_length_m):
    """
    从位置数据和给定的摆长计算单摆角度
    返回: t, theta, L(米), pivot(x_c,y_c)
    """
    t = data['time'].values
    x = data['x'].values  # 像素坐标
    y = data['y'].values  # 像素坐标
    
    # 拟合圆找到支点和半径（像素单位）
    x_c, y_c, r_pixels = fit_circle(x, y)
    
    # 根据给定的摆长计算像素到米的转换系数
    mm_per_pixel = (pendulum_length_m * 1000) / r_pixels  # 毫米每像素
    
    # 转换为毫米坐标
    x_mm = x * mm_per_pixel
    y_mm = y * mm_per_pixel
    x_c_mm = x_c * mm_per_pixel 
    y_c_mm = y_c * mm_per_pixel
    
    # 计算角度
    theta = np.arctan2(x_mm - x_c_mm, y_mm - y_c_mm)
    t = t - t[0]  # 时间归零

    logger.info(f"给定摆长: {pendulum_length_m:.4f} 米")
    logger.info(f"拟合得到的支点坐标: ({x_c:.2f}, {y_c:.2f}) 像素")
    logger.info(f"转换后的支点坐标: ({x_c_mm:.2f}, {y_c_mm:.2f}) 毫米")
    logger.info(f"计算得到的像素转换系数: {mm_per_pixel:.4f} mm/pixel")
    
    return t, theta, pendulum_length_m, (x_c_mm, y_c_mm)


def damped_oscillator(t, theta0, gamma, omega, phi):
    """阻尼简谐振动模型"""
    return theta0 * np.exp(-gamma * t) * np.cos(omega * t + phi)


def fit_pendulum_model(t, theta):
    """拟合阻尼简谐振动模型到实验数据 -> (popt, pcov)"""
    theta0_guess = np.max(np.abs(theta))
    gamma_guess = 0.1
    indices = np.where(np.diff(np.signbit(theta)))[0]
    if len(indices) >= 2:
        T = 2 * (t[indices[1]] - t[indices[0]])
    else:
        T = (t[-1] - t[0]) / 2
    omega_guess = 2 * np.pi / T
    phi_guess = 0

    p0 = [theta0_guess, gamma_guess, omega_guess, phi_guess]
    bounds = ([0, 0, 0, -np.pi], [np.inf, np.inf, np.inf, np.pi])

    try:
        popt, pcov = curve_fit(
            damped_oscillator, t, theta, p0=p0, bounds=bounds, maxfev=10000
        )
        return popt, pcov
    except Exception as e:
        logger.exception(f"拟合过程中出错: {e}")
        return None, None


def identify_periods(t, theta):
    """识别单摆的每个周期并计算每个周期的最大摆角"""
    zero_crossings = []
    for i in range(1, len(theta)):
        if theta[i-1] > 0 and theta[i] <= 0:
            t_interp = t[i-1] + (t[i] - t[i-1]) * (0 - theta[i-1]) / (theta[i] - theta[i-1])
            zero_crossings.append((i, t_interp))

    periods = []
    for i in range(len(zero_crossings) - 1):
        start_idx, start_t = zero_crossings[i]
        end_idx, end_t = zero_crossings[i+1]
        period_indices = range(start_idx, end_idx)
        period_t = t[period_indices]
        period_theta = theta[period_indices]
        period_duration = end_t - start_t
        max_amp = np.max(np.abs(period_theta))
        periods.append({
            'start_t': start_t,
            'end_t': end_t,
            'duration': period_duration,
            'max_amplitude': max_amp,
            'indices': period_indices
        })
    return periods


def calculate_large_angle_g(L, periods, gamma):
    """计算同时考虑大摆角效应和阻尼效应的重力加速度"""
    g_values, weights = [], []
    for period in periods:
        T = period['duration']
        theta_max = period['max_amplitude']
        g_approx = (4 * np.pi**2 * L) / T**2
        omega0_approx = np.sqrt(g_approx / L)
        angle_correction = 1 + theta_max**2/16 + 11*theta_max**4/3072
        damping_factor = 1 + gamma**2/(8 * omega0_approx**2)
        combined_correction = angle_correction * damping_factor
        g = (4 * np.pi**2 * L) / (T**2 / combined_correction**2)
        weight = 1.0 / max(0.01, theta_max)
        g_values.append(g)
        weights.append(weight)

    if g_values:
        g_corrected = np.average(g_values, weights=weights)
        deviations = [abs(g - g_corrected) / g_corrected * 100 for g in g_values]
        max_dev = max(deviations) if deviations else 0
        logger.info("\n同时考虑大摆角和阻尼效应的重力加速度计算:")
        logger.info(f"阻尼系数 γ = {gamma:.4f} s⁻¹; 周期数量: {len(periods)}; 最大偏差: {max_dev:.2f}%")
        return g_corrected
    return None


def calculate_physical_parameters(L, popt):
    """
    从拟合参数计算物理参数
    返回: g, b_over_m, g_no_damping
    """
    theta0, gamma, omega, phi = popt
    g = L * (omega**2 + gamma**2)
    b_over_m = 2 * gamma
    omega0 = np.sqrt(omega**2 + gamma**2)
    g_no_damping = L * omega0**2
    return g, b_over_m, g_no_damping


plt.rcParams['font.sans-serif'] = ['PingFang HK']


def visualize_results(t, theta, theta_smooth, theta_fit, fitted_params, data=None, pivot=None, mm_per_pixel=None, output_path='pendulum_fit_analysis.png'):
    """可视化实验数据和拟合结果"""
    fig = plt.figure(figsize=(12, 18))
    try:
        ax1 = fig.add_subplot(3, 1, 1)
        ax1.plot(t, theta, 'bo', alpha=0.3, markersize=2, label='原始数据')
        ax1.plot(t, theta_smooth, 'g-', label='平滑后数据')
        ax1.plot(t, theta_fit, 'r-', label='拟合曲线')
        ax1.set_xlabel('时间 t (s)')
        ax1.set_ylabel('摆角 θ (rad)')
        ax1.set_title('单摆摆角随时间的变化')
        ax1.legend(); ax1.grid(True)

        ax2 = fig.add_subplot(3, 1, 2)
        residuals = theta_smooth - theta_fit
        ax2.plot(t, residuals, 'k.')
        ax2.axhline(y=0, color='r', linestyle='-')
        ax2.set_xlabel('时间 t (s)'); ax2.set_ylabel('残差 (rad)'); ax2.set_title('拟合残差'); ax2.grid(True)

        if data is not None and pivot is not None and mm_per_pixel is not None:
            ax3 = fig.add_subplot(3, 1, 3)
            # 原始像素坐标
            x_pixels = data['x'].values; y_pixels = data['y'].values
            # 转换为毫米
            x_mm = x_pixels * mm_per_pixel; y_mm = y_pixels * mm_per_pixel
            x_c_mm, y_c_mm = pivot
            ax3.plot(x_mm, y_mm, 'b.', markersize=2, label='轨迹点')
            ax3.plot(x_c_mm, y_c_mm, 'ro', markersize=8, label='拟合支点')
            theta_circle = np.linspace(0, 2*np.pi, 100)
            r_mm = np.sqrt(np.mean((x_mm - x_c_mm)**2 + (y_mm - y_c_mm)**2))
            x_circle = x_c_mm + r_mm * np.cos(theta_circle)
            y_circle = y_c_mm + r_mm * np.sin(theta_circle)
            ax3.plot(x_circle, y_circle, 'r-', label=f'拟合圆 (r={r_mm:.2f}mm)')
            ax3.set_xlabel('x (mm)'); ax3.set_ylabel('y (mm)'); ax3.set_title('单摆轨迹和拟合圆')
            ax3.legend(); ax3.axis('equal'); ax3.grid(True)
        fig.savefig(output_path, dpi=300)
    finally:
        plt.close(fig)


def visualize_physical_parameters(L, popt, g, b_over_m, g_no_damping, g_corrected=None, output_path='pendulum_parameters_summary.png'):
    """创建一个展示所有物理参数的可视化页面"""
    theta0, gamma, omega, phi = popt
    T = 2 * np.pi / omega
    f = 1 / T
    decay_time = 1 / gamma
    Q = omega / (2 * gamma)

    fig = plt.figure(figsize=(10, 8))
    try:
        gs = plt.GridSpec(2, 2, height_ratios=[2, 1])
        ax1 = plt.subplot(gs[0, :]); ax1.axis('off')
        param_names = [
            "摆长 L\n(拟合得到的单摆长度)",
            "初始振幅 θ₀\n(摆角的最大值)",
            "阻尼系数 γ\n(振幅衰减速率)",
            "相对阻尼系数 b/m\n(摩擦力与质量的比值)",
            "角频率 ω\n(单位时间内角度变化)",
            "周期 T\n(完成一次完整摆动的时间)",
            "频率 f\n(每秒摆动次数)",
            "衰减时间 T_d\n(振幅衰减至1/e所需时间)",
            "品质因数 Q\n(振动系统的品质指标)",
        ]
        param_values = [
            f"{L*100:.1f} cm",
            f"{theta0:.4g} rad",
            f"{gamma:.4g} s^-1",
            f"{b_over_m:.4g} s^-1",
            f"{omega:.4g} rad/s",
            f"{T:.4g} s",
            f"{f:.4g} Hz",
            f"{decay_time:.4g} s",
            f"{Q:.4g}",
        ]
        table_data = list(zip(param_names, param_values))
        table = ax1.table(cellText=table_data, colLabels=["参数", "数值"], loc='center', cellLoc='center', colWidths=[0.6, 0.3])
        table.auto_set_font_size(False); table.set_fontsize(10); table.scale(1, 2.2)
        for (i, j), cell in table.get_celld().items():
            cell.set_linewidth(1)
            if i == 0:
                cell.set_facecolor('#E6E6FA'); cell.set_text_props(weight='bold')
            else:
                cell.set_facecolor('#F8F8FF')
        ax1.set_title('单摆物理参数', fontsize=16, pad=20)

        ax2 = plt.subplot(gs[1, :])
        g_methods = ["拟合模型", "无阻尼修正", "大摆角+阻尼修正" if g_corrected else "无大摆角修正"]
        g_values = [g, g_no_damping, g_corrected if g_corrected else 0]
        g_standard = 9.8
        g_errors = [(gv - g_standard) / g_standard * 100 for gv in g_values]
        bars = ax2.bar(g_methods, g_values, width=0.4)
        ax2.axhline(y=g_standard, color='k', linestyle='--', label='标准值 (9.8 m/s²)')
        for i, bar in enumerate(bars):
            height = bar.get_height(); error = g_errors[i]
            ax2.text(bar.get_x() + bar.get_width()/2, height, f'{height:.3f}\n({error:+.1f}%)', ha='center', va='bottom', fontsize=9)
        ax2.set_ylabel('重力加速度 (m/s²)'); ax2.grid(axis='y', linestyle='--', alpha=0.7); ax2.legend()
        fig.tight_layout(); fig.savefig(output_path, dpi=300)
    finally:
        plt.close(fig)


# Helper: expose the same savgol name expected by your original server code
savgol_filter = _savgol_filter