SOR
161 lines
5.9 KiB
Python
161 lines
5.9 KiB
Python
import math
|
|
from decimal import Decimal, ROUND_HALF_EVEN
|
|
|
|
"""
|
|
Speed of Sound Calculator (Successive Difference Method)
|
|
声速测量计算器 (逐差法)
|
|
|
|
参考了 Newton's rings.py 的高精度与不确定度修约逻辑。
|
|
|
|
计算公式说明:
|
|
1. 理论声速公式:
|
|
v_t = 331.45 * sqrt(1 + t / 273.15)
|
|
其中 t 为摄氏温度。
|
|
|
|
2. 逐差法计算波长 λ:
|
|
本实验两组数据均为 12 个点，分成两组进行逐差：Δx_i = x_{i+6} - x_i (i=1..6)
|
|
- 干涉法 (驻波法) 找振幅最大处：相邻两次最大处间距为 λ/2
|
|
因此 Δx = 6 * (λ/2) = 3λ => λ = Δx_avg / 3
|
|
- 相位法 (李萨如图像) 变为相同方向直线：相邻两次相位改变 2π，间距为 λ
|
|
因此 Δx = 6 * λ => λ = Δx_avg / 6
|
|
|
|
3. 声速公式:
|
|
v = f * λ
|
|
|
|
4. 不确定度传递公式:
|
|
- Δx_avg 的 A 类不确定度: u_A(Δx_avg) = S_Δ / sqrt(M) (其中 S_Δ 为逐差后 Δx 的标准差, M=6)
|
|
- 坐标 x 的 B 类不确定度: u_B(x) = Δ_inst / sqrt(3) (假设均匀分布)
|
|
- 逐差平均值的 B 类不确定度传播: u_B(Δx_avg) = sqrt(2/M) * u_B(x)
|
|
- Δx_avg 的合成不确定度: u_c(Δx_avg) = sqrt(u_A^2 + u_B^2)
|
|
- 波长 λ 的不确定度:
|
|
干涉法 u(λ) = u_c(Δx_avg) / 3
|
|
相位法 u(λ) = u_c(Δx_avg) / 6
|
|
- 声速 v 的不确定度 (忽略 f 的极小不确定度): u(v) = f * u(λ)
|
|
"""
|
|
|
|
def scientific_round(value, uncertainty):
|
|
"""
|
|
按照“四舍六入五凑偶”修约：
|
|
1. 不确定度 u 保留一位有效数字
|
|
2. 平均值末位与 u 对齐
|
|
"""
|
|
if uncertainty == 0:
|
|
return str(value), "0"
|
|
|
|
u_float = float(uncertainty)
|
|
first_digit_pos = math.floor(math.log10(u_float))
|
|
prec = Decimal('1e' + str(first_digit_pos))
|
|
|
|
u_final = uncertainty.quantize(prec, rounding=ROUND_HALF_EVEN)
|
|
val_final = value.quantize(prec, rounding=ROUND_HALF_EVEN)
|
|
|
|
return val_final, u_final
|
|
|
|
def calculate_method(method_name, f, k, delta_inst):
|
|
print(f"\n========================================")
|
|
print(f" {method_name} ")
|
|
print(f"========================================")
|
|
print("请输入12个坐标数据 (单位: mm)。")
|
|
|
|
x_vals = []
|
|
for i in range(1, 13):
|
|
while True:
|
|
try:
|
|
val_str = input(f"x_{i:02d} = ").strip()
|
|
if not val_str:
|
|
continue
|
|
x_vals.append(Decimal(val_str))
|
|
break
|
|
except Exception:
|
|
print("无效输入！请输入正确的数值。")
|
|
|
|
# 1. 逐差计算
|
|
M = 6
|
|
deltas = []
|
|
print("\n--- 逐差过程 ---")
|
|
for i in range(M):
|
|
diff = x_vals[i + M] - x_vals[i]
|
|
deltas.append(diff)
|
|
print(f"Δx_{i+1} = x_{i+M+1} - x_{i+1} = {x_vals[i+M]} - {x_vals[i]} = {diff} mm")
|
|
|
|
mean_delta = sum(deltas) / Decimal(str(M))
|
|
|
|
# 2. 波长与声速计算
|
|
# 将毫米转换为米计算速度
|
|
wavelength_mm = mean_delta / Decimal(str(k))
|
|
wavelength_m = wavelength_mm / Decimal("1000")
|
|
velocity = Decimal(str(f)) * wavelength_m
|
|
|
|
# 3. 不确定度计算
|
|
# u_A
|
|
variance_delta = sum((d - mean_delta)**2 for d in deltas) / Decimal(str(M - 1))
|
|
s_delta = Decimal(str(math.sqrt(float(variance_delta))))
|
|
u_a_delta = s_delta / Decimal(str(math.sqrt(M)))
|
|
|
|
# u_B
|
|
u_x = delta_inst / Decimal(str(math.sqrt(3)))
|
|
u_b_delta = Decimal(str(math.sqrt(2.0 / M))) * u_x
|
|
|
|
# u_c(Δx_avg)
|
|
u_delta_combined = Decimal(str(math.sqrt(float(u_a_delta**2 + u_b_delta**2))))
|
|
|
|
# 传播到波长和声速
|
|
u_wavelength_mm = u_delta_combined / Decimal(str(k))
|
|
u_wavelength_m = u_wavelength_mm / Decimal("1000")
|
|
u_velocity = Decimal(str(f)) * u_wavelength_m
|
|
|
|
# 4. 修约与展示
|
|
v_final, uv_final = scientific_round(velocity, u_velocity)
|
|
|
|
print("\n--- 中间统计结果 (高精度) ---")
|
|
print(f"平均差值 Δx_avg: {mean_delta:.6f} mm")
|
|
print(f"Δx 的标准差 S_Δ: {s_delta:.6f} mm")
|
|
print(f"Δx_avg 的 A 类不确定度 u_A: {u_a_delta:.6f} mm")
|
|
print(f"Δx_avg 的 B 类不确定度 u_B: {u_b_delta:.6f} mm")
|
|
print(f"波长 λ: {wavelength_mm:.6f} mm")
|
|
print(f"波长的不确定度 u(λ): {u_wavelength_mm:.6f} mm")
|
|
|
|
print("\n========================================")
|
|
print(f" {method_name} 最终结果 (四舍六入五凑偶)")
|
|
print(f" v = {v_final} ± {uv_final} m/s")
|
|
print("========================================")
|
|
|
|
def main():
|
|
print("========================================")
|
|
print(" 声速测量计算器 ")
|
|
print("========================================")
|
|
|
|
t_input = input("请输入环境温度 t (℃) [默认 24]: ").strip() or "24"
|
|
f_input = input("请输入谐振波频率 f (Hz) [默认 37000]: ").strip() or "37000"
|
|
delta_inst_input = input("请输入游标卡尺/测微鼓轮分度值误差 Δ_inst (mm) [默认 0.02]: ").strip() or "0.02"
|
|
|
|
t = float(t_input)
|
|
f = float(f_input)
|
|
delta_inst = Decimal(delta_inst_input)
|
|
|
|
# 计算理论声速
|
|
v_t = 331.45 * math.sqrt(1 + t / 273.15)
|
|
print(f"\n理论声速 v_t = {v_t:.2f} m/s (温度 {t} ℃)")
|
|
|
|
while True:
|
|
print("\n请选择计算方法:")
|
|
print("1. 干涉法 (驻波法) - 振幅最大处")
|
|
print("2. 相位法 (李萨如法) - 相同方向直线")
|
|
print("3. 退出")
|
|
choice = input("请输入选项 (1/2/3): ").strip()
|
|
|
|
if choice == '1':
|
|
# 干涉法 k = 3 (Δx_avg = 3λ)
|
|
calculate_method("干涉法 (振幅最大处)", f, 3, delta_inst)
|
|
elif choice == '2':
|
|
# 相位法 k = 6 (Δx_avg = 6λ)
|
|
calculate_method("相位法 (相同方向直线)", f, 6, delta_inst)
|
|
elif choice == '3':
|
|
print("退出程序。")
|
|
break
|
|
else:
|
|
print("无效选项，请重新输入。")
|
|
|
|
if __name__ == "__main__":
|
|
main()
|