跳转至

大物实验数据处理脚本

[!SUMMARY]

本人大学物理实验处理数据的脚本，仅供学习参考（尤其是数据！不要抄袭！），如出现错误脚本提供者不承担任何责任，使用时希望大家顺便检查是否有误，可在评论区留言；同时欢迎留言 / 私戳提供自己的脚本（希望有一定的中 / 英文注释），方式包括但不限于直接提供脚本 / 提供对应的链接。

[!ATTENTION]

鉴于不是所有人都有 python(3) 运行环境，以及有时代码过长影响观感，代码将放在我的 LeetCode/playground 中，在那里可以直接修改数据运行，或者是 fork 到你自己的 playground 中。当然，由于是在线环境，部分需要导入的模块缺失，可能导致运行错误。

或者看 github.com/darstib/physics_exp，文件名纯机翻。

声速的测定

LeetCode

measurementOfTheSpeedOfSound.py
# 驻波法/相位差法 读数
read1 = [11.335, 15.540, 19.805, 24.110, 28.550, 33.002, 37.625, 42.065]
read2 = [38.805, 43.065, 47.390, 51.625, 56.000, 60.305, 64.655, 69.085]

# 频率
f = 40.1
# 温度
t = 25.2
# 不确定度 ub
ub = 0.030

# 获得半波长
def get_lambda_2(read):
    l = len(read)
    m = l // 2
    # print(m)
    sum1 = sum(read[:m])
    sum2 = sum(read[m:])
    # print(sum1, sum2)
    return (sum2 - sum1) / (m**2)


# 获得测量声速
def get_v(read, f):
    lambda_2 = get_lambda_2(read)
    # print(lambda_2)
    return 2 * lambda_2 * f

# 两种方法的声速
v1 = get_v(read1, f)
v2 = get_v(read2, f)

print("测得声速：", v1, v2)

from math import sqrt

# 获得温度为t时的声速
def get_v(t):
    return 331.45 * sqrt(t / 273.15 + 1)

vt = get_v(t)

# 测量相对误差
def get_err(v, vt):
    return abs(v - vt) / vt

print("实际声速：", vt)
err1 = get_err(v1, vt)
err2 = get_err(v2, vt)

print("相对误差：", err1, err2)

# 计算平均值
def get_overline(read):
    return sum(read) / len(read)

# 计算不确定度 ua
def get_ua(read):
    read_overline = get_overline(read)
    n = len(read)
    return sqrt(sum([(x - read_overline) ** 2 for x in read]) / ((n - 1) * n))

# print("测得 U_a ：", get_ua(read1, read1_overline), get_ua(read2, read2_overline))

# 获得多次测量的 lambda
def get_lams(read):
    l = len(read)
    m = l // 2
    lams = []
    for i in range(m):
        lams.append((read[i + m] - read[i]) * 2)
    return lams

# A类不确定度
ua1 = get_ua(get_lams(read1))
ua2 = get_ua(get_lams(read2))

print("测得 U_a ：", ua1, ua2)


def get_u(ua, ub):
    return sqrt(ua**2 + ub**2)


u1 = get_u(ua1, ub)
u2 = get_u(ua2, ub)

print("不确定度：", u1, u2)

光速测量

LeetCode

measurementOfTheSpeedOfLight.py
# 测量数据（使用空格是因为打字空格比逗号方便）
## 第一组 方波形 
s1_1 = "0.04 0.05 0.07 0.10 0.15 0.20"
S1_1 = [float(i) for i in s1_1.split()]
s2_1 = "0.44 0.45 0.37 0.50 0.45 0.40"
S2_1 = [float(i) for i in s2_1.split()]
dt_1 = "58.5 58.0 44.5 58.5 45.0 29.5"
dt_1s = [float(i) for i in dt_1.split()]
## 第二组 正弦波形
s1_1 = "0.01 0.03 0.05 0.07 0.09 0.11"
S1_2 = [float(i) for i in s1_1.split()]
s2_1 = "0.51 0.43 0.45 0.47 0.49 0.51"
S2_2 = [float(i) for i in s2_1.split()]
dt_2 = "73.0 59.5 57.5 59.0 59.5 58.0"
dt_2s = [float(i) for i in dt_2.split()]

# 频率
v = 1 * 10**8
v_ = 4.55 * 10**5 # v_ 相当于 v'


# 计算光速
def get_cs(S1, S2, dt_s):
    def get_c(ds, dt_):
        return (2 * ds * v) / (dt_ * v_)

    cs = []
    for s1, s2, dt_ in zip(S1, S2, dt_s):
        cs.append(get_c(s2 - s1, dt_))
    return cs


cs1 = get_cs(S1_1, S2_1, dt_1s)
cs2 = get_cs(S1_2, S2_2, dt_2s)
print("第一组：", cs1)
print("第二组：", cs2)

# 平均光速
def get_avg(cs):
    return sum(cs) / len(cs)

c1_avg = get_avg(cs1)
c2_avg = get_avg(cs2)
print("平均光速", c1_avg, c2_avg)

# 测量相对误差
def get_err(c_avg):
    c = 3  # * 10**8
    return abs(c - c_avg) / c

err1 = get_err(c1_avg)
err2 = get_err(c2_avg)
print("相对误差：", err1, err2)

# 计算不确定度 ua
def get_ua(cs):
    from math import sqrt

    c_avg = get_avg(cs)
    n = len(cs)
    return sqrt(sum([(c - c_avg) ** 2 for c in cs]) / ((n - 1) * n))

ua1 = get_ua(cs1)
ua2 = get_ua(cs2)
print("测得 U_a ：", ua1, ua2)

# 本题我们老师没说 u_b 在哪，我参考的报告也没写，我也就没算了……

金属材料杨氏模量的测定

数据处理部分理解有问题，但是可能只是 hkk 老师这样要求？就不放出来了。

鹤翔万里的处理代码

动态法测量的杨氏模量

Leetcode

4_youngSModulus2.py
# ZJU 动态法测量材料杨氏模量 数据处理脚本
# 处理 f
def get_f(fs):

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.optimize import curve_fit

    xs = [5 * i for i in range(1, 11) if i != 7]

    # Define the quadratic function for fitting
    def quadratic(x, a, b, c):
        return a * x**2 + b * x + c

    # Perform curve fitting
    popt, _ = curve_fit(quadratic, xs, fs)
    # Generate points for smooth curve
    x_smooth = np.linspace(min(xs), max(xs), 100)
    y_smooth = quadratic(x_smooth, *popt)
    # Find the minimum point
    x_min = -popt[1] / (2 * popt[0])
    y_min = quadratic(x_min, *popt)
    # Plotting
    plt.figure(figsize=(10, 6))
    plt.scatter(xs, fs, color="red", label="Data points")
    plt.plot(x_smooth, y_smooth, color="blue", label="Fitted curve")
    plt.plot(x_min, y_min, "go", label="Minimum point")
    plt.xlabel("x (mm)")
    plt.ylabel("Frequency (Hz)")
    plt.title("Frequency vs Distance")
    plt.legend()
    plt.grid(True)
    # Annotate the minimum point
    plt.annotate(
        f"Min: ({x_min:.2f}, {y_min:.2f})",
        xy=(x_min, y_min),
        xytext=(5, 5),
        textcoords="offset points",
    )
    plt.savefig("4_youngSModulus2.png")
    print(f"Minimum point: x = {x_min:.2f}, f = {y_min:.2f}")
    return y_min


fs = [736.4, 734.8, 733.8, 733.1, 732.5, 731.8, 732.2, 732.3, 732.5]
# f_fix = get_f(fs)
f_fix = 731.86

##### 处理读数 #####
d0 = -0.015  # 螺旋测微器零点读数 单位: mm
ds = [5.950, 5.945, 5.948, 5.955, 5.952]  # 单位: mm
Ls = [160.5, 160.1, 160.2, 160.1, 160.1]  # 单位: mm
ms = [37.338, 37.337, 37.336, 37.337, 37.337]  # 单位: g，可以只测一次，自己设置 m_avg

from math import sqrt


def get_avg(li):
    return sum(li) / len(li)


d_avg = get_avg(ds) - d0
L_avg = get_avg(Ls)
m_avg = get_avg(ms)
print("d_avg:", d_avg, "L_avg:", L_avg, "m_avg:", m_avg)


# dl 为仪器误差
def get_U(li, dl):
    # get U_a
    n = len(li)
    i_avg = get_avg(li)
    U_a = sqrt(sum([(i - i_avg) ** 2 for i in li]) / (n * (n - 1)))
    U_b = dl / sqrt(3)
    U = sqrt(U_a**2 + U_b**2)
    return U_a, U_b, U


# L 不确定度
dL = 0.2
dd = 0.004
dm = 0.001
df = 0.1
L_a, L_b, U_L = get_U(Ls, dL)
print("L_a:", L_a, "L_b:", L_b, "U_L:", U_L)
# d 不确定度
d_a, d_b, U_d = get_U(ds, dd)
print("d_a:", d_a, "d_b:", d_b, "U_d:", U_d)

E = [
    1.6067 * (L**3) * m * (f_fix**2) / (d**4) / (1e11) for L, m, d in zip(Ls, ms, ds)
]  # 单位: 10**11 N/m^2

E_avg = get_avg(E)
print("E_avg: ", E_avg)

E = 1.6067 * (L_avg**3) * m_avg * (f_fix**2) / (d_avg**4) / (1e11)
print("E:", E)
# 计算相对不确定度
dE = E_avg * sqrt(
    (
        (3 * dL / L_avg) ** 2
        + (4 * dd / d_avg) ** 2
        + (dm / m_avg) ** 2
        + (2 * df / f_fix) ** 2
    )
)

print("dE:", dE)

普朗克常量（Matlab）

phylab-plangke 来自 @ 吉水飞云

空气密度测量

Leetcode

5_measurementOfAirDensity
# ZJU 空气密度测量 数据处理脚本
# p 可能是 “气压”，也可能是 “密度”，请依据注释区分
t = 22.0  # ℃
pw0 = 2643.38  # Pa ，饱和水蒸气压
n = 0.615  # 相对湿度
pw = pw0 * n  # Pa ，水蒸气压
print("pw:", pw)

V = 159.326  # cm^3 玻璃泡的体积
dms = [0.1893, 0.1892, 0.1891]  # g m1-m0 空气质量
ps = [m / V for m in dms]  # g/cm^3 空气密度
print("ps:", ps)
avg_p = sum(ps) / len(ps) * 1e3  # g/cm^3 => kg/m^3 平均空气密度
print("平均密度：", avg_p, "kg/m^3")
p0 = 101325  # Pa ，标准大气压
p_ = 102510  # Pa ，气压
p = p_ * (1 - 0.000163)  # Pa 修正后的实验气压
print("p: ", p_)
a = 1 / 273.15
T0 = 273.15  # K
p_gan = avg_p * (p0 / p) * (a * t + 1) * ((3 * pw) / (8 * p) + 1)
print("p_gan:", p_gan)  # 标况下干燥气体的密度
Ma = 0.02898  # kg/mol ，空气的摩尔质量
R = (p0 * Ma) / (T0 * p_gan)  # J/(kg·K) 普适气体常数
print("R:", R)

print(p / 133.322 / 10)  # cmHg 附录中查找理论温度

抛射体运动的照相法研究

参考 ZJU-General-Physics-Experiment-I 略长，仅放 Leetcode 上。

评论