commit 8ce7705ae858ae1b5604e9c9dd384a4269ebb640 Author: SOR <3272194354@qq.com> Date: Sun Jun 7 10:49:28 2026 +0800 init diff --git a/Bridge resistance measurement.py b/Bridge resistance measurement.py new file mode 100644 index 0000000..309ff16 --- /dev/null +++ b/Bridge resistance measurement.py @@ -0,0 +1,370 @@ +import math +from decimal import Decimal, ROUND_HALF_EVEN + +""" +Wheatstone Bridge Resistance Measurement Calculator +电桥法测量电阻计算器 + +参考了表面张力系数计算器的高精度与不确定度修约逻辑。 + +公式说明: +1. 惠斯通电桥平衡条件 (Bridge Balance Equation): + 当检流计示数为零时，电桥平衡： + Rx / R0 = R1 / R2 + 即: + Rx = R0 * (R1 / R2) （也写作 Rx = R0 * C，C = R1/R2 为比例臂比值） + + 其中: + Rx : 待测电阻 + R0 : 比较臂标准电阻（电桥读数） + R1/R2 : 比例臂之比，取 1:1, 1:10, 10:1 + +2. 灵敏度公式 (Bridge Sensitivity): + S = Δn / (ΔR0 / R0) + + 其中: + S : 电桥灵敏度（格/相对变化量） + Δn : 检流计偏转格数（通常取 5 小格） + ΔR0 : 使检流计偏转 Δn 格时 R0 的改变量 + ΔR0 = |R0' - R0|，R0' 为偏转 5 格时的读数 + R0 : 平衡时的标准电阻读数 + + 也可写为: + S = Δn * R0 / ΔR0 + +3. 不确定度传递公式 (Uncertainty Propagation): + Rx = R0 * C (C 为比例臂比值) + + (a) R0 的不确定度: + u(R0) = sqrt( u_A(R0)^2 + u_B(R0)^2 ) + 其中: + u_A(R0) = s(R0) / sqrt(n) (A 类不确定度, s 为标准偏差) + u_B(R0) = Δ_inst / sqrt(3) (B 类不确定度, Δ_inst 为仪器误差限) + + (b) 仪器误差限 Δ_inst 的确定 (与灵敏度有关): + Δ_inst = R0 / S + 即由灵敏度决定的电桥最小可分辨电阻变化量 + + (c) 比例臂 C 视为精确值（标准电阻箱），则: + u(Rx) = C * u(R0) + + (d) 相对不确定度: + u_r(Rx) = u(Rx) / Rx = u(R0) / R0 + + (e) 修约规则: + 不确定度保留一位有效数字，四舍六入五凑偶 + Rx 末位与不确定度对齐 +""" + +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_stats(data_list, delta_inst): + """ + 计算一组数据的平均值、A类、B类及合成不确定度 + u(X) = sqrt( u_A^2 + u_B^2 ) + u_A = s / sqrt(n) + u_B = Δ_inst / sqrt(3) + """ + n = len(data_list) + mean = sum(data_list) / Decimal(n) + + if n > 1: + variance = sum((x - mean)**2 for x in data_list) / Decimal(n - 1) + s = Decimal(str(math.sqrt(float(variance)))) + u_a = s / Decimal(str(math.sqrt(n))) + else: + u_a = Decimal("0") + + u_b = Decimal(delta_inst) / Decimal(str(math.sqrt(3))) + u_combined = Decimal(str(math.sqrt(float(u_a**2 + u_b**2)))) + + return mean, u_a, u_b, u_combined + +def input_decimal(prompt, default=None): + """输入一个 Decimal 值，支持默认值""" + while True: + try: + suffix = f" (默认 {default})" if default else "" + val = input(f"{prompt}{suffix}: ").strip() + if not val and default is not None: + return Decimal(str(default)) + if not val: + continue + return Decimal(val) + except Exception: + print("输入无效，请输入数字。") + +def input_int(prompt, default=None): + """输入一个整数，支持默认值""" + while True: + try: + suffix = f" (默认 {default})" if default else "" + val = input(f"{prompt}{suffix}: ").strip() + if not val and default is not None: + return default + return int(val) + except Exception: + print("输入无效，请输入整数。") + +def measure_one_ratio(ratio_label, C_value, R_nominal, n_measurements): + """ + 对单一比例臂进行测量与分析 + + 参数: + ratio_label: 比例臂标签，如 "1:1" + C_value: 比例臂比值 C = R1/R2 + R_nominal: 标定电阻标称值 (Ω) + n_measurements: 测量次数 + """ + C = Decimal(str(C_value)) + + print(f"\n{'─'*50}") + print(f" 比例臂 Rx:R0 = {ratio_label}，比例因子 C = {C_value}") + print(f" 标定电阻标称值: {R_nominal} Ω") + print(f"{'─'*50}") + + # --- 输入平衡时的 R0 读数 --- + print(f"\n请输入 {n_measurements} 次平衡时的 R0 读数 (Ω):") + R0_data = [] + for i in range(n_measurements): + while True: + try: + val = input(f" 第 {i+1} 次 R0 = ").strip() + if not val: continue + R0_data.append(Decimal(val)) + break + except Exception: + print(" 输入无效，请输入数字。") + + # --- 输入偏差 5 小格时的 R0' --- + print(f"\n请输入偏转 5 小格时的 R0' 读数 (Ω):") + R0_prime = input_decimal(" R0'") + delta_n = Decimal("5") # 检流计偏转格数 + + # --- 计算 R0 的平均值 --- + R0_mean = sum(R0_data) / Decimal(len(R0_data)) + + # --- 计算灵敏度 --- + # S = Δn * R0 / |R0' - R0| + delta_R0 = abs(R0_prime - R0_mean) + if delta_R0 == 0: + print(" ⚠ 警告: R0' = R0，无法计算灵敏度，请检查数据！") + S = Decimal("Infinity") + delta_inst = Decimal("0") + else: + S = delta_n * R0_mean / delta_R0 + # 仪器误差限 (由灵敏度决定) + # Δ_inst = R0 / S = |R0' - R0| / Δn = ΔR0 / 5 + delta_inst = R0_mean / S + + # --- 不确定度分析 --- + R0_mean_final, u_a_R0, u_b_R0, u_R0 = calculate_stats(R0_data, str(delta_inst)) + + # Rx = C * R0 + Rx_mean = C * R0_mean_final + u_Rx = C * u_R0 + + # 相对不确定度 + if Rx_mean != 0: + u_r_Rx = u_Rx / Rx_mean + else: + u_r_Rx = Decimal("0") + + # --- 修约 --- + Rx_final, u_Rx_final = scientific_round(Rx_mean, u_Rx) + + # --- 输出结果 --- + print(f"\n{'═'*50}") + print(f" 结果 (比例臂 {ratio_label}, 标定 {R_nominal} Ω)") + print(f"{'═'*50}") + + # 打印 R0 测量数据表 + print(f"\n R0 测量数据:") + for i, val in enumerate(R0_data): + print(f" 第 {i+1} 次: {val} Ω") + + print(f"\n R0 平均值 = {R0_mean_final:.4f} Ω") + print(f" R0' = {R0_prime} Ω") + print(f" ΔR0 = |R0'-R0| = {delta_R0:.4f} Ω") + print(f" 偏转格数 Δn = {delta_n} 格") + print(f"{'─'*50}") + print(f" 灵敏度 S = Δn·R0/ΔR0 = {S:.2f}") + print(f" 仪器误差限 Δ_inst = R0/S = {delta_inst:.4f} Ω") + print(f"{'─'*50}") + print(f" u_A(R0) = {u_a_R0:.6f} Ω") + print(f" u_B(R0) = {u_b_R0:.6f} Ω") + print(f" u(R0) = {u_R0:.6f} Ω") + print(f"{'─'*50}") + print(f" Rx = C × R0 = {C} × {R0_mean_final:.4f} = {Rx_mean:.4f} Ω") + print(f" u(Rx) = C × u(R0) = {u_Rx:.6f} Ω") + print(f" 相对不确定度 u_r = {u_r_Rx*100:.2f}%") + print(f"{'─'*50}") + print(f" ★ 最终结果: Rx = {Rx_final} ± {u_Rx_final} Ω") + print(f"{'═'*50}") + + return { + "ratio_label": ratio_label, + "C": C, + "R_nominal": R_nominal, + "R0_data": R0_data, + "R0_mean": R0_mean_final, + "R0_prime": R0_prime, + "delta_R0": delta_R0, + "S": S, + "delta_inst": delta_inst, + "u_a_R0": u_a_R0, + "u_b_R0": u_b_R0, + "u_R0": u_R0, + "Rx_mean": Rx_mean, + "u_Rx": u_Rx, + "u_r_Rx": u_r_Rx, + "Rx_final": Rx_final, + "u_Rx_final": u_Rx_final + } + +def main(): + print("════════════════════════════════════════════════") + print(" 电桥法测量电阻计算器 (Wheatstone Bridge)") + print("════════════════════════════════════════════════") + print() + print("实验原理:") + print(" 惠斯通电桥平衡条件: Rx = R0 × C (C = R1/R2)") + print(" 灵敏度: S = Δn × R0 / ΔR0") + print(" 不确定度: u(Rx) = C × u(R0), u(R0) = √(u_A² + u_B²)") + print() + + # --- 1. 基本参数 --- + n_meas = input_int("每组测量次数 n", default=5) + + # --- 2. 定义测量方案 --- + # 标定电阻: 10kΩ 和 2kΩ + # 比例臂: 1:1, 1:10, 10:1 + # 注意: 2kΩ 不会有 1:10 的比例臂 (因为 R0 = Rx/(C) 会超出电阻箱量程) + + calibration_plan = [ + # (标定电阻标称值, 比例臂标签, C 值) + # ---- 10kΩ 标定 ---- + ("10000", "1:1", 1), + ("10000", "10:1", 10), + # ---- 2kΩ 标定 ---- + ("2000", "1:1", 1), + ("2000", "10:1", 10), + # ---- 200Ω 标定 ---- + ("200", "1:1", 1), + ("200", "1:10", 0.1), + ("200", "10:1", 10), + ] + + print(f"\n{'━'*55}") + print(" 测量方案如下:") + print(f"{'━'*55}") + print(f" {'序号':<6}{'标定电阻':>10}{'比例臂 Rx:R0':>14}{'C = R1/R2':>12}") + print(f" {'─'*46}") + for idx, (R_nom, ratio, C_val) in enumerate(calibration_plan, 1): + print(f" {idx:<6}{R_nom + ' Ω':>10}{ratio:>14}{C_val:>12}") + print(f"{'━'*55}") + print() + print(" 注: 2kΩ 标定不进行 1:10 比例臂测量") + print(" (R0 = 2000/0.1 = 20000 Ω，超出电阻箱量程)") + print() + + # --- 3. 逐项测量 --- + all_results = [] + + for idx, (R_nom, ratio_label, C_val) in enumerate(calibration_plan, 1): + print(f"\n{'▓'*55}") + print(f" 第 {idx} 组: 标定 {R_nom} Ω，比例臂 Rx:R0 = {ratio_label}") + print(f"{'▓'*55}") + + result = measure_one_ratio(ratio_label, C_val, R_nom, n_meas) + all_results.append(result) + + if idx < len(calibration_plan): + input("\n按回车继续下一组测量...") + + # --- 4. 汇总表 --- + print("\n") + print("╔══════════════════════════════════════════════════════════════════════════════╗") + print("║ 测 量 结 果 汇 总 表 ║") + print("╠══════════════════════════════════════════════════════════════════════════════╣") + + header = (f" {'标定':>6} {'比例臂':>8} {'R0均值(Ω)':>12} " + f"{'S(灵敏度)':>10} {'Δ_inst(Ω)':>10} " + f"{'Rx(Ω)':>12} {'u(Rx)(Ω)':>10} {'u_r(%)':>8}") + print(header) + print(f" {'─'*74}") + + for res in all_results: + S_str = f"{res['S']:.1f}" if res['S'] != Decimal("Infinity") else "∞" + line = (f" {res['R_nominal']:>6} {res['ratio_label']:>8} " + f"{res['R0_mean']:>12.4f} {S_str:>10} " + f"{res['delta_inst']:>10.4f} " + f"{res['Rx_mean']:>12.4f} {res['u_Rx']:>10.6f} " + f"{res['u_r_Rx']*100:>7.2f}%") + print(line) + + print("╚══════════════════════════════════════════════════════════════════════════════╝") + + # --- 5. 灵敏度分析总结 --- + print("\n") + print("┌──────────────────────────────────────────────┐") + print("│ 灵 敏 度 分 析 总 结 │") + print("├──────────────────────────────────────────────┤") + + for res in all_results: + S_str = f"{res['S']:.2f}" if res['S'] != Decimal("Infinity") else "∞" + print(f"│ 标定 {res['R_nominal']:>5}Ω 比例臂 {res['ratio_label']:>4} " + f"S = {S_str:>8} │") + + print("├──────────────────────────────────────────────┤") + + # 找灵敏度最大和最小 + finite_results = [r for r in all_results if r['S'] != Decimal("Infinity")] + if finite_results: + best = max(finite_results, key=lambda r: r['S']) + worst = min(finite_results, key=lambda r: r['S']) + print(f"│ 最高灵敏度: 标定{best['R_nominal']:>5}Ω 比例臂{best['ratio_label']:>4} │") + print(f"│ S = {best['S']:.2f} │") + print(f"│ 最低灵敏度: 标定{worst['R_nominal']:>5}Ω 比例臂{worst['ratio_label']:>4} │") + print(f"│ S = {worst['S']:.2f} │") + + print("├──────────────────────────────────────────────┤") + print("│ 分析: │") + print("│ · 灵敏度 S 越大，测量精度越高 │") + print("│ · S = Δn·R0/ΔR0 │") + print("│ · 比例臂比值影响灵敏度 │") + print("│ · 合理选择比例臂可优化测量精度 │") + print("└──────────────────────────────────────────────┘") + + # --- 6. 最终修约结果汇总 --- + print("\n") + print("╔═══════════════════════════════════════════════════╗") + print("║ 最 终 结 果 (修约后) ║") + print("╠═══════════════════════════════════════════════════╣") + + for res in all_results: + print(f"║ 标定 {res['R_nominal']:>5}Ω 比例臂 {res['ratio_label']:>4}: ║") + print(f"║ Rx = {res['Rx_final']} ± {res['u_Rx_final']} Ω ║") + print(f"║ (相对不确定度 {res['u_r_Rx']*100:.2f}%) ║") + print(f"╟───────────────────────────────────────────────────╢") + + print("╚═══════════════════════════════════════════════════╝") + +if __name__ == "__main__": + main() diff --git a/Newton's rings.py b/Newton's rings.py new file mode 100644 index 0000000..4db7654 --- /dev/null +++ b/Newton's rings.py @@ -0,0 +1,140 @@ +import math +from decimal import Decimal, ROUND_HALF_EVEN + +""" +Newton's Rings Curvature Radius Calculator (Successive Difference Method) +牛顿环曲率半径计算器 (逐差法) + +参考 lenth.py 的高精度与不确定度修约逻辑。 +""" + +def get_decimal_places(data_strings): + """获取输入数据中小数点后的保留位数，用于自适应展示""" + max_places = 0 + for s in data_strings: + if '.' in s: + places = len(s.split('.')[1]) + max_places = max(max_places, places) + return max_places + +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 main(): + # --- 1. 参数与常量初始化 --- + WAVELENGTH = Decimal("0.0005893") # 钠光波长默认 589.3 nm = 0.0005893 mm + M = 10 # 逐差法的差值项数 (30-20=10) + + print("========================================") + print(" 牛顿环曲率半径计算 (逐差法) ") + print("========================================") + print(f"默认波长 λ = {WAVELENGTH} mm (钠光)") + print("参考环数: 第 11 环 到 第 30 环 (共 20 个数据)") + + # --- 2. 数据输入 --- + d_squared_vals = [] + data_strs = [] + print("\n[第一步] 请输入暗环直径平方 D^2 (单位: mm^2)") + print("支持粘贴或逐个输入，精确到小数点后三位。") + + for i in range(11, 31): + while True: + try: + user_input = input(f"D_{i:02d}^2 = ").strip() + if not user_input: + continue + val = Decimal(user_input) + d_squared_vals.append(val) + data_strs.append(user_input) + break + except Exception: + print("无效输入！请输入正确的直径平方值（如 12.345）。") + + print("\n[第二步] 仪器参数输入") + res_input = input("请输入仪器分度值 (默认 0.01 mm): ").strip() or "0.01" + delta_inst_input = input("请输入仪器误差 Δ_inst (默认 0.005 mm): ").strip() or "0.005" + + res = Decimal(res_input) + delta_inst = Decimal(delta_inst_input) + + # --- 3. 逐差计算 --- + # Δ_i = D_{i+M}^2 - D_i^2 (i = 11...20) + deltas = [] + sum_d_pairs = Decimal("0") # 用于 B 类不确定度传播: sum(D_low^2 + D_high^2) + + print("\n--- 逐差过程 ---") + for i in range(10): + d_low_sq = d_squared_vals[i] + d_high_sq = d_squared_vals[i + 10] + diff = d_high_sq - d_low_sq + deltas.append(diff) + sum_d_pairs += (d_low_sq + d_high_sq) + print(f"Δ_{i+11:02d} = D_{i+21}^2 - D_{i+11}^2 = {diff} mm^2") + + mean_delta = sum(deltas) / Decimal(str(M)) + + # --- 4. 曲率半径计算 --- + # R = Δ_avg / (4 * M * λ) + # M = 10, 所以分母是 40 * λ + denominator = Decimal("4") * Decimal(str(M)) * WAVELENGTH + r_val = mean_delta / denominator + + # --- 5. 不确定度计算 --- + # (1) A 类不确定度 u_A (针对逐差后的 Δ) + variance_delta = sum((x - mean_delta)**2 for x 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))) + + # (2) B 类不确定度 u_B (传播自读数显微镜的仪器误差) + # 推导过程: + # u(x) = delta_inst / sqrt(3) + # u(D^2) = sqrt( (2D u_x)^2 + (2D u_x)^2 ) = 2 * sqrt(2) * D * u_x + # u(D^2)^2 = 8 * D^2 * u_x^2 + # u(Δ)^2 = u(D_high^2)^2 + u(D_low^2)^2 = 8 * u_x^2 * (D_high^2 + D_low^2) + # u(Δ_avg)^2 = (1/M^2) * sum( u(Δ_i)^2 ) = (8 * u_x^2 / M^2) * sum_d_pairs + u_x = delta_inst / Decimal(str(math.sqrt(3))) + u_b_delta_sq = (Decimal("8") * (u_x**2) * sum_d_pairs) / (Decimal(str(M))**2) + u_b_delta = Decimal(str(math.sqrt(float(u_b_delta_sq)))) + + # (3) 合成不确定度 u(Δ) 并在 R 上传播 + u_delta_combined = Decimal(str(math.sqrt(float(u_a_delta**2 + u_b_delta**2)))) + u_r = u_delta_combined / denominator + + # --- 6. 修约与展示 --- + r_final, u_final = scientific_round(r_val, u_r) + + print("\n--- 中间统计结果 (高精度) ---") + print(f"平均差值 Δ_avg: {mean_delta:.6f} mm^2") + print(f"Δ 的标准差 s: {s_delta:.6f}") + print(f"Δ 的 A 类不确定度 u_A(Δ): {u_a_delta:.6f}") + print(f"Δ 的 B 类不确定度 u_B(Δ): {u_b_delta:.6f}") + + print("\n--- 最终计算结果 (R) ---") + print(f"曲率半径 R (原始): {r_val:.4f} mm") + print(f"R 的合成不确定度 u(R): {u_r:.4f} mm") + + print("\n========================================") + print(" 最终修约结果 (四舍六入五凑偶)") + print(f" R = {r_final} ± {u_final} mm") + print("========================================") + +if __name__ == "__main__": + main() diff --git a/README.md b/README.md new file mode 100644 index 0000000..98f4314 Binary files /dev/null and b/README.md differ diff --git a/Speed of sound.py b/Speed of sound.py new file mode 100644 index 0000000..6fc46b3 --- /dev/null +++ b/Speed of sound.py @@ -0,0 +1,160 @@ +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() diff --git a/Surface tension coefficient.py b/Surface tension coefficient.py new file mode 100644 index 0000000..21b998b --- /dev/null +++ b/Surface tension coefficient.py @@ -0,0 +1,257 @@ +import math +from decimal import Decimal, ROUND_HALF_EVEN + +""" +Surface Tension Coefficient Calculator (Pull-off Method) +拉脱法测液体表面张力系数计算器 + +参考了杨氏模量及粘滞系数代码的高精度与不确定度修约逻辑。 + +公式说明: +1. 表面张力系数计算公式: + α = F / [π * (d1 + d2)] + 其中: + F : 拉脱力，通过传感器测得 F = U / K (K为传感器灵敏度) + d1 : 铁环内径 + d2 : 铁环外径 + +2. 传感器灵敏度 K (通过砝码标定): + 测得 2 组不同砝码个数 i (0~5) 下的电压 U1_i 与 U2_i (通常为增重与减重过程)。 + 取其平均值 U_avg_i = (U1_i + U2_i) / 2。 + 设单个砝码质量为 m，则受力 F_i = i * m * g。 + 使用最小二乘法对 (F_i, U_avg_i) 进行线性回归，斜率即为灵敏度 K (V/N)。 + +3. 不确定度传递公式 (相对不确定度): + 假定 K 的不确定度 u_K 通过线性回归标准差求得。 + + 拉脱力 F 的平均值: F_mean = U_pull_mean / K + 拉脱力 F 的不确定度 (相对不确定度): + u_r(F) = sqrt( (u_U / U_pull_mean)^2 + (u_K / K)^2 ) + + 内外径和(d1+d2)的不确定度: + u_d_sum = sqrt( u_d1^2 + u_d2^2 ) + u_r(d_sum) = u_d_sum / (d1_mean + d2_mean) + + 总相对不确定度: + u_r(α) = u_α / α = sqrt( u_r(F)^2 + u_r(d_sum)^2 ) + + 绝对不确定度: + u_α = α * u_r(α) + + 对于多次测量量 X 的不确定度: + u_X = sqrt( u_A,X^2 + u_B,X^2 ) + u_A,X = s_X / sqrt(n) + u_B,X = Δ_inst / sqrt(3) +""" + +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_stats(data_list, delta_inst): + """计算一组数据的平均值、A类、B类及合成不确定度""" + n = len(data_list) + mean = sum(data_list) / Decimal(n) + + if n > 1: + variance = sum((x - mean)**2 for x in data_list) / Decimal(n - 1) + s = Decimal(str(math.sqrt(float(variance)))) + u_a = s / Decimal(str(math.sqrt(n))) + else: + u_a = Decimal("0") + + u_b = Decimal(delta_inst) / Decimal(str(math.sqrt(3))) + u_combined = Decimal(str(math.sqrt(float(u_a**2 + u_b**2)))) + + return mean, u_combined + +def linear_regression(x_list, y_list): + """ + 对 (x, y) 列表进行线性回归，返回斜率 k, 截距 b, 及斜率的标准不确定度 u_k + """ + n = len(x_list) + x_mean = sum(x_list) / Decimal(n) + y_mean = sum(y_list) / Decimal(n) + + L_xx = sum((x - x_mean)**2 for x in x_list) + L_yy = sum((y - y_mean)**2 for y in y_list) + L_xy = sum((x_list[i] - x_mean) * (y_list[i] - y_mean) for i in range(n)) + + # 斜率 + k = L_xy / L_xx + # 截距 + b = y_mean - k * x_mean + + # 相关系数 r + r = L_xy / Decimal(str(math.sqrt(float(L_xx * L_yy)))) + + # 斜率的标准不确定度 u_k + # u_k = k * sqrt( (1/r^2 - 1) / (n - 2) ) + r_float = float(r) + if abs(r_float) < 1.0: + u_k_float = float(k) * math.sqrt( (1.0 / (r_float**2) - 1.0) / (n - 2) ) + else: + u_k_float = 0.0 + u_k = Decimal(str(abs(u_k_float))) + + return k, b, r, u_k + +def input_data_group(name, count, unit=""): + """输入一组重复测量的数值""" + prompt_unit = f" (单位: {unit})" if unit else "" + print(f"\n请输入 {name} 的 {count} 次测量值{prompt_unit}:") + data = [] + for i in range(count): + while True: + try: + val = input(f"第 {i+1} 次 = ").strip() + if not val: continue + data.append(Decimal(val)) + break + except Exception: + print("输入无效，请输入数字。") + return data + +def main(): + print("========================================") + print(" 拉脱法测水表面张力系数计算器 ") + print("========================================") + + # --- 1. 常量设置 --- + G = Decimal("9.80") # 重力加速度 m/s^2 + PI = Decimal(str(math.pi)) + print(f"重力加速度 g = {G} m/s^2") + print(f"圆周率 π ≈ {PI:.5f}") + + print("\n--- 砝码参数与标定 ---") + m_default = "0.5" + m_input = input(f"单个小砝码质量 m (g, 默认 {m_default}): ").strip() or m_default + m_mass = Decimal(m_input) / Decimal("1000") # 转为 kg + + print("\n请输入力敏传感器标定电压 (共 2 组，0~5 个小砝码，例如：第 1 组为增重，第 2 组为减重):") + + print("\n--- 输入第 1 组标定电压 (U1) ---") + U_calib_1 = [] + for i in range(6): + while True: + try: + val = input(f"第1组 - 放 {i} 个砝码时的电压 U1 (V或mV, 请统一单位) = ").strip() + if not val: continue + U_calib_1.append(Decimal(val)) + break + except Exception: + print("输入无效，请输入数字。") + + print("\n--- 输入第 2 组标定电压 (U2) ---") + U_calib_2 = [] + for i in range(6): + while True: + try: + val = input(f"第2组 - 放 {i} 个砝码时的电压 U2 (V或mV, 请统一单位) = ").strip() + if not val: continue + U_calib_2.append(Decimal(val)) + break + except Exception: + print("输入无效，请输入数字。") + + # 计算两组的平均电压并生成对应的受力 + U_calib_avg = [] + F_calib = [] + for i in range(6): + avg_val = (U_calib_1[i] + U_calib_2[i]) / Decimal("2") + U_calib_avg.append(avg_val) + F_calib.append(Decimal(str(i)) * m_mass * G) + + # 打印标定数据对照表 + print("\n" + "-"*65) + print(f"{'砝码数':<6} | {'受力 F (N)':<12} | {'第1组 U1':<10} | {'第2组 U2':<10} | {'平均 U_avg':<10}") + print("-"*65) + for i in range(6): + print(f"{i:<8} | {F_calib[i]:<12.6f} | {U_calib_1[i]:<10.4f} | {U_calib_2[i]:<10.4f} | {U_calib_avg[i]:<10.4f}") + print("-"*65) + + # 计算灵敏度 K + K, K_b, r, u_K = linear_regression(F_calib, U_calib_avg) + print(f"\n[标定结果] 传感器灵敏度 K = {K:.6f}, 截距 = {K_b:.6f}") + print(f"线性相关系数 r = {r:.6f}, 斜率不确定度 u_K = {u_K:.6e}") + + # --- 2. 拉脱测量与仪器误差输入 --- + print("\n--- 输入仪器误差 Δ_inst ---") + delta_U = Decimal(input("电压 U 的仪器误差 (默认 0.001): ").strip() or "0.001") + delta_d = Decimal(input("游标卡尺测量内外径的仪器误差 (mm, 默认 0.02): ").strip() or "0.02") + + # 输入拉脱电压 U_pull (6次) + U_pull_vals = input_data_group("拉脱时的电压差/峰值电压 U_pull", 6, "同标定电压单位") + U_pull_mean, U_pull_u = calculate_stats(U_pull_vals, delta_U) + + # 输入内径 d1 (6次) + d1_vals = input_data_group("铁环内径 d1", 6, "mm") + d1_mean, d1_u = calculate_stats(d1_vals, delta_d) + + # 输入外径 d2 (6次) + d2_vals = input_data_group("铁环外径 d2", 6, "mm") + d2_mean, d2_u = calculate_stats(d2_vals, delta_d) + + # --- 3. 表面张力系数计算 --- + # F = U / K + F_pull_mean = U_pull_mean / K + + # 统一转化内径、外径为国际标准单位 m + d1_m = d1_mean / Decimal("1000") + d2_m = d2_mean / Decimal("1000") + + d1_u_m = d1_u / Decimal("1000") + d2_u_m = d2_u / Decimal("1000") + + # α = F / [π * (d1 + d2)] + d_sum = d1_m + d2_m + alpha_val = F_pull_mean / (PI * d_sum) + + # --- 4. 不确定度传播 --- + # u_r(F) = sqrt( (u_U / U_pull_mean)^2 + (u_K / K)^2 ) + u_r_F_sq = (U_pull_u / U_pull_mean)**2 + (u_K / K)**2 + u_r_F = Decimal(str(math.sqrt(float(u_r_F_sq)))) + + # u(d1 + d2) + u_d_sum = Decimal(str(math.sqrt(float(d1_u_m**2 + d2_u_m**2)))) + u_r_d_sum = u_d_sum / d_sum + + # u_r(α) = sqrt( u_r(F)^2 + u_r(d_sum)^2 ) + u_r_alpha = Decimal(str(math.sqrt(float(u_r_F**2 + u_r_d_sum**2)))) + + # 绝对不确定度 + alpha_u = alpha_val * u_r_alpha + + # --- 5. 修约与输出 --- + alpha_final, u_final = scientific_round(alpha_val, alpha_u) + + print("\n" + "="*40) + print(" 计算结果") + print("="*40) + print(f"拉脱电压 U_pull = {U_pull_mean:.4f} ± {U_pull_u:.4f}") + print(f"拉脱力 F = {F_pull_mean:.6f} N") + print(f"铁环内径 d1 = {d1_mean:.3f} ± {d1_u:.3f} mm") + print(f"铁环外径 d2 = {d2_mean:.3f} ± {d2_u:.3f} mm") + print("-" * 40) + print(f"相对不确定度 u_r = {u_r_alpha*100:.2f}%") + print(f"绝对不确定度 u_α = {alpha_u:.6f} N/m") + print("-" * 40) + print(f"最终表面张力系数: α = {alpha_final} ± {u_final} N/m") + print("="*40) + +if __name__ == "__main__": + main() diff --git a/Viscosity coefficient.py b/Viscosity coefficient.py new file mode 100644 index 0000000..e4a6e63 --- /dev/null +++ b/Viscosity coefficient.py @@ -0,0 +1,188 @@ +import math +from decimal import Decimal, ROUND_HALF_EVEN + +""" +Viscosity Coefficient Calculator (Stokes' Law) +蓖麻油粘滞系数计算器 (斯托克斯定律) + +参考了杨氏模量代码的高精度与不确定度修约逻辑。 +（注：由于实验提供的是同一距离20cm的5次重复时间测量，此处直接求平均计算，不适用逐差法。逐差法通常用于等间距的多组累积物理量测量。） + +公式说明: +1. 粘滞系数计算公式 (考虑管壁修正): + η = (ρ - ρ0) * g * d^2 * t / [18 * l * (1 + 2.4 * d/D)] + 其中: + ρ : 小球密度 + ρ0 : 液体(蓖麻油)密度 + d : 小球直径 + D : 玻璃管内径 + l : 下落距离 + t : 下落平均时间 + +2. 雷诺数 (Re) 判断和修正公式: + Re = ρ0 * v * d / η + 其中 v = l / t。 + 当 Re < 0.1 时，斯托克斯定律完全适用。 + 当 0.1 <= Re <= 1.0 时，可使用奥辛(Oseen)修正: + η_corrected = η / (1 + 3/16 * Re) + +3. 不确定度传递公式 (相对不确定度): + 假定 ρ, ρ0, g 的不确定度可忽略，则: + u_r(η) = u_η / η = sqrt( (2 * u_d / d)^2 + (u_t / t)^2 + (u_l / l)^2 + (u_D / D)^2 ) + + 绝对不确定度: + u_η = η * u_r(η) + + 其中各个量 X 的不确定度 u_X = sqrt( u_A,X^2 + u_B,X^2 ) + u_A,X = s_X / sqrt(n) + u_B,X = Δ_inst / sqrt(3) +""" + +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_stats(data_list, delta_inst): + """计算一组数据的平均值、A类、B类及合成不确定度""" + n = len(data_list) + mean = sum(data_list) / Decimal(n) + + if n > 1: + variance = sum((x - mean)**2 for x in data_list) / Decimal(n - 1) + s = Decimal(str(math.sqrt(float(variance)))) + u_a = s / Decimal(str(math.sqrt(n))) + else: + u_a = Decimal("0") + + u_b = Decimal(delta_inst) / Decimal(str(math.sqrt(3))) + u_combined = Decimal(str(math.sqrt(float(u_a**2 + u_b**2)))) + + return mean, u_combined + +def input_data_group(name, count, unit=""): + """输入一组重复测量的数值""" + print(f"\n请输入 {name} 的 {count} 次测量值" + (f" (单位: {unit})" if unit else "") + ":") + data = [] + for i in range(count): + while True: + try: + val = input(f"第 {i+1} 次 = ").strip() + if not val: continue + data.append(Decimal(val)) + break + except Exception: + print("输入无效，请输入数字。") + return data + +def main(): + print("========================================") + print(" 蓖麻油粘滞系数计算 (斯托克斯定律) ") + print("========================================") + + # --- 1. 常量设置 --- + G = Decimal("9.80") # 重力加速度 m/s^2 + RHO_BALL = Decimal("7800") # 钢球密度 kg/m^3 (默认为钢球，通常实验室用此值) + RHO_OIL = Decimal("960") # 蓖麻油密度 kg/m^3 (随温度变化，常温附近约 960) + + print(f"重力加速度 g = {G} m/s^2") + print(f"小球密度 ρ = {RHO_BALL} kg/m^3") + print(f"蓖麻油密度 ρ0 = {RHO_OIL} kg/m^3") + + # --- 2. 测量数据输入 --- + print("\n--- 仪器与实验参数设定 ---") + D_mean = Decimal("20.0") + d_mean = Decimal("1.0") + l_mean = Decimal("20.0") # 默认单位 cm，后续转化为 mm 或 m 计算 + + ans = input("使用默认参数 (玻璃管直径D=20mm, 小球直径d=1mm, 下落距离l=20cm) 吗？[Y/n]: ").strip().lower() + if ans == 'n': + D_mean = Decimal(input("玻璃管直径 D (mm): ").strip() or "20.0") + d_mean = Decimal(input("小球直径 d (mm): ").strip() or "1.0") + l_mean = Decimal(input("下落距离 l (cm): ").strip() or "20.0") + + l_mean_mm = l_mean * Decimal("10") # 统一 cm 转为 mm + + print("\n--- 输入仪器误差 Δ_inst ---") + delta_D = Decimal(input("D 的仪器误差 (mm, 默认 0): ").strip() or "0") + delta_d = Decimal(input("d 的仪器误差 (mm, 默认 0): ").strip() or "0") + delta_l = Decimal(input("l 的仪器误差 (mm, 默认 0): ").strip() or "0") + delta_t = Decimal(input("t(时间)的仪器误差 (s, 默认 0.01): ").strip() or "0.01") + + # D, d, l 我们假设只进行了单次测量或作为已知常数处理，所以只算B类不确定度 + D_u = delta_D / Decimal(str(math.sqrt(3))) + d_u = delta_d / Decimal(str(math.sqrt(3))) + l_u = delta_l / Decimal(str(math.sqrt(3))) + + # 输入5个时间数据 + t_vals = input_data_group("小球下落时间 t", 5, "s") + t_mean, t_u = calculate_stats(t_vals, delta_t) + + # --- 3. 粘滞系数计算 --- + # 统一转化为国际标准单位 (m, kg, s) + D_m = D_mean / Decimal("1000") + d_m = d_mean / Decimal("1000") + l_m = l_mean_mm / Decimal("1000") + + # η = (ρ - ρ0) * g * d^2 * t / [18 * l * (1 + 2.4 * d/D)] + correction_factor = Decimal("1") + Decimal("2.4") * (d_m / D_m) + numerator = (RHO_BALL - RHO_OIL) * G * (d_m**2) * t_mean + denominator = Decimal("18") * l_m * correction_factor + eta_val = numerator / denominator + + # --- 4. 雷诺数判断与修正 --- + # v = l / t + v_val = l_m / t_mean + # Re = ρ0 * v * d / η + Re = RHO_OIL * v_val * d_m / eta_val + + eta_corrected = eta_val + if Re > Decimal("0.1"): + # Oseen 修正 + eta_corrected = eta_val / (Decimal("1") + Decimal("3") / Decimal("16") * Re) + + # --- 5. 不确定度传播 --- + # 相对不确定度平方和 + rel_u_sq = (Decimal("2") * d_u / d_mean)**2 + (t_u / t_mean)**2 + (l_u / l_mean_mm)**2 + (D_u / D_mean)**2 + rel_u = Decimal(str(math.sqrt(float(rel_u_sq)))) + + # 绝对不确定度 + eta_u = eta_corrected * rel_u + + # --- 6. 修约与输出 --- + eta_final, u_final = scientific_round(eta_corrected, eta_u) + + print("\n" + "="*40) + print(" 计算结果") + print("="*40) + print(f"下落时间 t = {t_mean:.3f} ± {t_u:.3f} s") + print(f"下落速度 v = {v_val:.5f} m/s") + print("-" * 40) + print(f"雷诺数 Re = {Re:.4f}") + if Re <= Decimal("0.1"): + print("Re <= 0.1，满足斯托克斯定律条件，无需 Oseen 修正。") + else: + print("Re > 0.1，已应用 Oseen 修正。") + print(f"修正前 η = {eta_val:.4f} Pa·s") + print("-" * 40) + print(f"相对不确定度 u_r = {rel_u*100:.2f}%") + print(f"绝对不确定度 u_η = {eta_u:.4f} Pa·s") + print("-" * 40) + print(f"最终结果: η = {eta_final} ± {u_final} Pa·s") + print("="*40) + +if __name__ == "__main__": + main() diff --git a/Young's modulus.py b/Young's modulus.py new file mode 100644 index 0000000..383e49a --- /dev/null +++ b/Young's modulus.py @@ -0,0 +1,201 @@ +import math +from decimal import Decimal, ROUND_HALF_EVEN + +""" +Young's Modulus Calculator (Optical Lever & Successive Difference Method) +杨氏模量计算器 (光杠杆 & 逐差法) + +不确定度传递公式: +1. 杨氏模量公式: + E = (32 * m * g * L * D) / (π * d^2 * b * Δn_avg) + 其中 Δn_avg 是 4 个砝码对应的平均刻度位移（由 8 个刻度值通过逐差法求得）。 + +2. 相对不确定度传递公式: + u_r(E) = u_E / E = sqrt( (u_L/L)^2 + (u_D/D)^2 + (2 * u_d/d)^2 + (u_b/b)^2 + (u_Δn/Δn_avg)^2 ) + +3. 合成不确定度: + u_E = E * u_r(E) + 其中各物理量 X 的不确定度 u_X = sqrt( u_A,X^2 + u_B,X^2 ) + u_A,X = s_X / sqrt(n) + u_B,X = Δ_inst / sqrt(3) +""" + +def scientific_round(value, uncertainty): + """ + 按照“四舍六入五凑偶”修约： + 1. 不确定度 u 保留一位有效数字（若首位是1或2，有时保留两位，此处统一保留一位） + 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_stats(data_list, delta_inst): + """计算一组数据的平均值、A类、B类及合成不确定度""" + n = len(data_list) + mean = sum(data_list) / Decimal(n) + + # A类不确定度 + if n > 1: + variance = sum((x - mean)**2 for x in data_list) / Decimal(n - 1) + s = Decimal(str(math.sqrt(float(variance)))) + u_a = s / Decimal(str(math.sqrt(n))) + else: + u_a = Decimal("0") + + # B类不确定度 + u_b = Decimal(delta_inst) / Decimal(str(math.sqrt(3))) + + # 合成 + u_combined = Decimal(str(math.sqrt(float(u_a**2 + u_b**2)))) + + return mean, u_combined + +def input_data_group(name, count, unit="mm"): + """输入一组重复测量的数值""" + print(f"\n请输入 {name} 的 {count} 次测量值 (单位: {unit}):") + data = [] + for i in range(count): + while True: + try: + val = input(f"{name}_{i+1} = ").strip() + if not val: continue + data.append(Decimal(val)) + break + except Exception: + print("输入无效，请输入数字。") + return data + +def main(): + # --- 1. 常量设置 --- + M_WEIGHT = Decimal("0.350") # 单个砝码质量 350g = 0.350kg + G = Decimal("9.80") # 重力加速度 + PI = Decimal(str(math.pi)) + + print("========================================") + print(" 杨氏模量计算 (光杠杆 & 逐差法) ") + print("========================================") + print(f"单个砝码质量 m = {M_WEIGHT} kg") + print(f"重力加速度 g = {G} m/s^2") + + # --- 2. 几何长度输入 (各6次测量) --- + # L: 钢丝长度, D: 镜面到刻度尺距离, d: 钢丝直径, b: 光杠杆臂长 + L_vals = input_data_group("钢丝长度 L", 6, "mm") + D_vals = input_data_group("镜面到尺距离 D", 6, "mm") + d_vals = input_data_group("钢丝直径 d", 6, "mm") + b_vals = input_data_group("光杠杆臂长 b", 6, "mm") + + print("\n--- 输入仪器误差 Δ_inst ---") + delta_L = Decimal(input("L 的仪器误差 (mm, 默认 1.0): ").strip() or "1.0") + delta_D = Decimal(input("D 的仪器误差 (mm, 默认 1.0): ").strip() or "1.0") + delta_d = Decimal(input("d 的仪器误差 (mm, 默认 0.005): ").strip() or "0.005") + delta_b = Decimal(input("b 的仪器误差 (mm, 默认 0.02): ").strip() or "0.02") + + L_mean, L_u = calculate_stats(L_vals, delta_L) + D_mean, D_u = calculate_stats(D_vals, delta_D) + d_mean, d_u = calculate_stats(d_vals, delta_d) + b_mean, b_u = calculate_stats(b_vals, delta_b) + + # --- 3. 刻度值输入 (0-7个砝码, 加重+减重) --- + print("\n--- 输入刻度尺读数 n (单位: mm) ---") + n_up = [] + n_down = [] + for i in range(8): + while True: + try: + line = input(f"砝码数量 {i} (加重, 减重): ").replace(',', ' ').split() + if len(line) != 2: + print("请输入两个数值（空格或逗号隔开）") + continue + n_up.append(Decimal(line[0])) + n_down.append(Decimal(line[1])) + break + except Exception: + print("输入无效！") + + delta_n_inst = Decimal(input("刻度尺仪器误差 (mm, 默认 0.5): ").strip() or "0.5") + + # 计算平均读数 n_i + n_means = [(u + d) / 2 for u, d in zip(n_up, n_down)] + + print("\n--- 中间量 (加重与减重的平均值 n_i) ---") + for i, val in enumerate(n_means): + print(f"n_{i} = {val} mm") + + # 逐差法计算 Δn (隔 4 项相减) + # Δn_j = n_{j+4} - n_j, j = 0,1,2,3 + dn_list = [] + print("\n--- 逐差计算过程 ---") + for j in range(4): + diff = n_means[j+4] - n_means[j] + dn_list.append(diff) + print(f"Δn_{j} = n_{j+4} - n_{j} = {diff} mm") + + # Δn 的平均值及不确定度 + dn_mean = sum(dn_list) / Decimal("4") + # A类: 针对这 4 个差值 + variance_dn = sum((x - dn_mean)**2 for x in dn_list) / Decimal("3") + s_dn = Decimal(str(math.sqrt(float(variance_dn)))) + u_a_dn = s_dn / Decimal(str(math.sqrt(4))) + + # B类: 刻度尺误差 + # 因为 n_i = (n_up + n_down)/2, Δn = n_high - n_low + # 传播得 u_B(Δn) = Δ_inst / sqrt(3) + u_b_dn = delta_n_inst / Decimal(str(math.sqrt(3))) + dn_u = Decimal(str(math.sqrt(float(u_a_dn**2 + u_b_dn**2)))) + + # --- 4. 杨氏模量计算 --- + # 全部转为标准单位 (m) + L_m = L_mean / 1000 + D_m = D_mean / 1000 + d_m = d_mean / 1000 + b_m = b_mean / 1000 + dn_m = dn_mean / 1000 + + # E = (32 * m * g * L * D) / (pi * d^2 * b * Δn) + # 这里的 Δn 对应 4 个砝码的位移 + numerator = Decimal("32") * M_WEIGHT * G * L_m * D_m + denominator = PI * (d_m**2) * b_m * dn_m + E_val = numerator / denominator + + # --- 5. 不确定度传播 --- + # rel_u = sqrt( (uL/L)^2 + (uD/D)^2 + (2*ud/d)^2 + (ub/b)^2 + (udn/dn)^2 ) + rel_u_sq = (L_u/L_mean)**2 + (D_u/D_mean)**2 + (Decimal("2")*d_u/d_mean)**2 + (b_u/b_mean)**2 + (dn_u/dn_mean)**2 + rel_u = Decimal(str(math.sqrt(float(rel_u_sq)))) + E_u = E_val * rel_u + + # --- 6. 修约与输出 --- + # 通常用 10^11 Pa 表示 + E_11 = E_val / Decimal("1e11") + E_u_11 = E_u / Decimal("1e11") + + E_final, u_final = scientific_round(E_11, E_u_11) + + print("\n" + "="*40) + print(" 计算结果") + print("="*40) + print(f"L = {L_mean:.2f} ± {L_u:.2f} mm") + print(f"D = {D_mean:.2f} ± {D_u:.2f} mm") + print(f"d = {d_mean:.4f} ± {d_u:.4f} mm") + print(f"b = {b_mean:.3f} ± {b_u:.3f} mm") + print(f"Δn = {dn_mean:.3f} ± {dn_u:.3f} mm (4个砝码位移)") + print("-" * 40) + print(f"杨氏模量 E (原始): {E_val:.4e} Pa") + print(f"相对不确定度 u_r: {rel_u*100:.2f}%") + print("-" * 40) + print(f"最终结果: E = {E_final} ± {u_final} × 10^11 Pa") + print("="*40) + +if __name__ == "__main__": + main() diff --git a/__pycache__/Bridge resistance measurement.cpython-313.pyc b/__pycache__/Bridge resistance measurement.cpython-313.pyc new file mode 100644 index 0000000..4963430 Binary files /dev/null and b/__pycache__/Bridge resistance measurement.cpython-313.pyc differ diff --git a/__pycache__/Surface tension coefficient.cpython-313.pyc b/__pycache__/Surface tension coefficient.cpython-313.pyc new file mode 100644 index 0000000..b5a909f Binary files /dev/null and b/__pycache__/Surface tension coefficient.cpython-313.pyc differ diff --git a/degree.py b/degree.py new file mode 100644 index 0000000..44ccb23 --- /dev/null +++ b/degree.py @@ -0,0 +1,130 @@ +import math +from decimal import Decimal, ROUND_HALF_EVEN + +def parse_angle_to_minutes(s): + """ + 将输入解析为内部统一单位：分钟（Decimal）。 + 支持格式： + 1. 'D.M' -> D度M分 (例如 '120.15' -> 120*60 + 15 = 7215分) + 2. 'D' (纯数字无点) -> D度 (例如 '120' -> 120*60 = 7200分) + 3. 'D.M.S' -> D度M分S秒 (例如 '120.15.30' -> 7215.5分) + """ + parts = s.split('.') + if len(parts) == 1: + # 视为度 + return Decimal(parts[0]) * 60 + elif len(parts) == 2: + # 度.分 + return Decimal(parts[0]) * 60 + Decimal(parts[1]) + else: + # 度.分.秒 或更多 + res = Decimal(parts[0]) * 60 + Decimal(parts[1]) + if len(parts) >= 3: + res += Decimal(parts[2]) / 60 + return res + +def format_minutes_as_angle(minutes, prec=None): + """ + 将总分钟数转换为 D° M' 格式，支持保持对齐精度。 + """ + neg = minutes < 0 + minutes = abs(minutes) + + d = int(minutes // 60) + m = minutes % 60 + + # 进位处理：如果 m 由于位宽限制进位到 60 + if m >= 60: + d += 1 + m -= 60 + + if prec is not None and prec < 1: + # 计算小数位数 + places = abs(prec.as_tuple().exponent) + m_str = format(m, f'.{places}f') + else: + # 整数或自动小数位 + if m == int(m): + m_str = str(int(m)) + else: + m_str = format(m, 'f').rstrip('0').rstrip('.') + + return f"{'-' if neg else ''}{d}°{m_str}'" + +def calc_uncertainty(data_strs, res_str="1"): + # 转换为分钟进行高精度计算 + data = [parse_angle_to_minutes(s) for s in data_strs] + res = Decimal(res_str) # 仪器误差，单位：分 + n = len(data) + + # 1. 计算平均值 + sum_val = sum(data) + mean = sum_val / n + + # 2. 计算标准差 s + if n > 1: + variance = sum((x - mean) ** 2 for x in data) / (n - 1) + s = Decimal(str(math.sqrt(float(variance)))) + else: + s = Decimal('0') + + # 3. A类不确定度 u_A + u_A = s / Decimal(str(math.sqrt(n))) + + # 4. B类不确定度 u_B (仪器误差) + u_B = res / Decimal(str(math.sqrt(3))) + + # 5. 合成不确定度 u + u = Decimal(str(math.sqrt(float(u_A**2 + u_B**2)))) + + return mean, s, u_A, u_B, u + +def format_final_result(mean_minutes, u_minutes): + """ + 按照“四舍六入五凑偶”修约： + 1. 不确定度 u 保留一位有效数字 + 2. 平均值末位与 u 对齐 + """ + u_float = float(u_minutes) + if u_float == 0: + return f"θ = {format_minutes_as_angle(mean_minutes)} ± 0'" + + # 找到第一位非零数字的位置 + first_digit_pos = math.floor(math.log10(u_float)) + prec = Decimal('1e' + str(first_digit_pos)) + + # 修约 + u_final = u_minutes.quantize(prec, rounding=ROUND_HALF_EVEN) + mean_final_val = mean_minutes.quantize(prec, rounding=ROUND_HALF_EVEN) + + # 转换为角度字符串输出 + mean_str = format_minutes_as_angle(mean_final_val, prec) + u_str = f"{u_final}'" + + return f"θ = {mean_str} ± {u_str}" + +if __name__ == "__main__": + print("=== 角度不确定度计算器 ===") + print("输入格式: '度.分' (如 120.15 表示 120°15')，多个数据用空格隔开。") + + try: + user_input = input("请输入数据：").split() + if not user_input: + print("未输入数据！") + else: + res_input = input("请输入仪器分度值 (单位: 分，默认 1)：").strip() or "1" + + mean_v, s, u_A, u_B, u = calc_uncertainty(user_input, res_input) + + print("\n--- 中间计算过程 ---") + print(f"平均值: {format_minutes_as_angle(mean_v)} ({mean_v:.4f}')") + print(f"标准差 s: {s:.6f}'") + print(f"A类不确定度 u_A: {u_A:.6f}'") + print(f"B类不确定度 u_B: {u_B:.6f}'") + print(f"合成不确定度 u: {u:.6f}'") + + print("\n--- 最终修约结果 (四舍六入五凑偶) ---") + print(format_final_result(mean_v, u)) + + except Exception as e: + print(f"错误: {e}") \ No newline at end of file diff --git a/lenth.py b/lenth.py new file mode 100644 index 0000000..21c9496 --- /dev/null +++ b/lenth.py @@ -0,0 +1,78 @@ +import math +from decimal import Decimal, ROUND_HALF_EVEN + +def get_decimal_places(data_strings): + """获取输入数据中小数点后的位数""" + max_places = 0 + for s in data_strings: + if '.' in s: + places = len(s.split('.')[1]) + max_places = max(max_places, places) + return max_places + +def calc_uncertainty(data_strs, res_str="0.05"): + # 转换为 Decimal 进行高精度计算 + data = [Decimal(s) for s in data_strs] + res = Decimal(res_str) + n = len(data) + + # 1. 计算平均值 + places = get_decimal_places(data_strs) + sum_val = sum(data) + mean = sum_val / n + # 修改：平均值修约使用“四舍六入五凑偶” + mean_final = mean.quantize(Decimal('1.' + '0' * places), rounding=ROUND_HALF_EVEN) + + # 2. 计算标准差 s + variance = sum((x - mean) ** 2 for x in data) / (n - 1) + s = Decimal(str(math.sqrt(float(variance)))) + + # 3. A类不确定度 u_A + u_A = s / Decimal(str(math.sqrt(n))) + + # 4. B类不确定度 u_B (仪器误差) + u_B = res / Decimal(str(math.sqrt(3))) + + # 5. 合成不确定度 u + u = Decimal(str(math.sqrt(float(u_A**2 + u_B**2)))) + + return mean_final, s, u_A, u_B, u + +def format_final_result(mean, u): + """ + 按照“四舍六入五凑偶”修约： + 1. 不确定度 u 保留一位有效数字 + 2. 平均值末位与 u 对齐 + """ + u_float = float(u) + # 找到第一位非零数字的位置 + first_digit_pos = math.floor(math.log10(u_float)) + prec = Decimal('1e' + str(first_digit_pos)) + + # 修改：不确定度 u 的修约改为 ROUND_HALF_EVEN + u_final = u.quantize(prec, rounding=ROUND_HALF_EVEN) + + # 修改：平均值修约保持一致 + mean_final = mean.quantize(prec, rounding=ROUND_HALF_EVEN) + + return f"f = {mean_final} ± {u_final}" + +if __name__ == "__main__": + # ... (后续输入输出逻辑保持不变) + print("请输入数据，用空格分隔：") + raw_input = input().split() + + print("请输入仪器分度值（默认 0.05）：") + res_input = input().strip() or "0.05" + + mean_f, s, u_A, u_B, u = calc_uncertainty(raw_input, res_input) + + print("\n--- 中间计算过程 (高精度展示) ---") + print(f"平均值 (原始精度): {mean_f}") + print(f"标准差 s: {s:.6f}") + print(f"A类不确定度 u_A: {u_A:.6f}") + print(f"B类不确定度 u_B: {u_B:.6f}") + print(f"合成不确定度 u: {u:.6f}") + + print("\n--- 最终修约结果 (四舍六入五凑偶) ---") + print(format_final_result(mean_f, u)) \ No newline at end of file