292 lines
11 KiB
Python
292 lines
11 KiB
Python
import logging
|
||
import numpy as np
|
||
import pandas as pd
|
||
import matplotlib
|
||
matplotlib.use('Agg') # headless backend
|
||
import matplotlib.pyplot as plt
|
||
from matplotlib import font_manager
|
||
from scipy.optimize import curve_fit
|
||
from scipy.signal import savgol_filter as _savgol_filter
|
||
import os
|
||
|
||
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
|
||
|
||
|
||
# 设置中文字体
|
||
font_path = './fonts/PingFangSC-Medium.ttf'
|
||
if os.path.exists(font_path):
|
||
prop = font_manager.FontProperties(fname=font_path)
|
||
plt.rcParams['font.family'] = prop.get_name()
|
||
plt.rcParams['axes.unicode_minus'] = False
|
||
logger.info(f"成功加载字体文件: {font_path}")
|
||
else:
|
||
# 备选方案
|
||
plt.rcParams['font.sans-serif'] = ['PingFang SC', 'SimHei', 'DejaVu Sans']
|
||
plt.rcParams['axes.unicode_minus'] = False
|
||
logger.warning(f"字体文件不存在: {font_path}, 使用系统字体")
|
||
|
||
|
||
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
|
||
|
||
|
||
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 |