光电效应及普朗克常数的测定数据处理 Python实现

内容仅供参考，如有错误，欢迎指正，如有疑问，欢迎交流。 

因为我不会Excel所以只能用Python来处理

祝大家早日摆脱物理实验的苦海

用到的一些方法

PCHIP （分段三次埃尔米特插值多项式）

因为实验时记录的数据比较少，记录的也比较随便，直接将这些数据连起来的话画出的图像十分的直，而伏安特性曲线顾名思义应该是一条曲线，因此我们需要进行数据拟合，但是呢我们发现直接进行数据拟合会出现下面的情况，即在某些区间内，电压上升电流反而下降了，这很不合理，所以我让AI帮我想到了PCHIP插值法进行数据点的扩展，具体是啥我也不知道，反正用就完了。

利用导数找拐点

看了示范报告，一开始把零点当成拐点算了，爆炸了

找拐点的方法我试了好多种，二阶导等于零、一阶导的极值、一阶导作差、一阶导和各种阈值比较等等，最后我用的是区间斜率与整体斜率的比较，也算是一阶导和阈值比较的一种吧。

这些方法的核心思想都是在伏安特性曲线上找到电流变化陡峭的拐点，并将其对应的电压作为截止电压，关键是怎样才算陡峭。本文用到的方法中，先计算电流随电压的平均变化率（即平均斜率），作为基准值 threshold =(I_max - I_min) / (U_max - U_min)，随后遍历所有插值后的数据点，利用五点差分法计算其局部斜率，如果这个局部斜率是基准斜率的1.5倍，就将其认为是拐点，记录并退出循环。但是因为找到了第一个就退出循环了并且“1.5倍”不一定是每组数据的最佳选择，所以结果可能会出现较大误差。

threshold = (max(y_dense) - min(y_dense)) / (x_max - x_min)  # 导数阈值
                
for i in range(nums_data):
    if (y_dense[i + 5] - y_dense[i]) / (x_dense[i + 5] - x_dense[i]) > threshold * 1.5:
        cutoff_voltage = x_dense[i]
        break

代码

import openpyxl
import matplotlib.pyplot as plt
import numpy as np
import os
from matplotlib.ticker import FormatStrFormatter
from scipy.interpolate import PchipInterpolator
from docx import Document
from docx.shared import Inches

# 设置支持中文的字体和正常显示负号
plt.rcParams['font.sans-serif'] = ['SimHei']  # 或者 ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False

# ------------------------------------------控制台--------------------------------------------------
N = 5  # 数据组数

cutoff_voltages = [] # 截止电压 (V)
calculate = 1 # 数据计算+图片输出
save_to_word = 0 # 将图片保存到 Word 文档

output_dir = r"C:\Users\86138\Desktop\这是文件夹\这是文件夹里的文件夹" # 指定图片保存目录
excel_file = r"C:\Users\86138\Desktop\这是文件夹\这是xlsx.xlsx"  # Excel 文件路径 (用于加载原始数据和储存计算后的数据)
word_file = r"C:\Users\86138\Desktop\这是文件夹\这是docx.docx"  # Word 文件路径（用于读取图片并保存）

# --------------------------------------------------------------------------------------------------

auto = 0 if cutoff_voltages else 1
os.makedirs(output_dir, exist_ok=True)  # 确保目标文件夹存在

def Htlang():
    # 物理常数
    c = 3.0e8  # 光速，单位 m/s

    # 1. 加载 Excel 文件
    wb = openpyxl.load_workbook(excel_file)
    sheet = wb.active

    # 2. 确定数据组的行间隔（第一组从第 1 行开始，每组占 2 行）
    group_start_rows = [1 + i * 2 for i in range(N)]  # 生成 [1, 3, 5, 7, 9]

    # 3. 存储所有数据以便后续绘图
    wavelengths = []     # 波长 (nm)
    frequencies = []     # 频率 (10^-14 Hz)

    # 4. 处理每一组数据
    for group_idx, start_row in enumerate(group_start_rows, start=1):
        # 读取表格名称、坐标轴标签
        table_name = sheet.cell(row=start_row, column=1).value  # A列存储波长信息（例如 "365nm"）
        x_label = sheet.cell(row=start_row, column=2).value       # B列第一行: 横坐标标签（电压）
        y_label = sheet.cell(row=start_row + 1, column=2).value     # B列第二行: 纵坐标标签（电流）

        # **处理 y_label，将乘方变为上标**
        if y_label:
            y_label = y_label.replace("10-10", r"$10^{-10}$") \
                            .replace("10-11", r"$10^{-11}$") \
                            .replace("10-12", r"$10^{-12}$")

        # 读取电压与电流数据（从 C 到 L 列）
        voltage_data = []
        current_data = []
        for col in range(3, 13):  # C(3)到L(12)
            v = sheet.cell(row=start_row, column=col).value  # 电压数据
            i = sheet.cell(row=start_row + 1, column=col).value  # 电流数据
            voltage_data.append(v)
            current_data.append(i)

        # 计算频率
        # 从 table_name 中去除末尾两个字符（例如 "nm"）转换为 int 得到波长
        wavelength_nm = int(table_name[:-2])
        wavelength_m = wavelength_nm * 1e-9  # 转换为 m
        frequency = c / wavelength_m  # 计算频率 (Hz)
        frequency_14 = frequency / 1e14  # 以 10^-14 Hz 为单位

        wavelengths.append(wavelength_nm)
        frequencies.append(frequency_14)

        # 过滤并排序有效数据
        valid_indices = [i for i, (v, i_val) in enumerate(zip(voltage_data, current_data))
                        if v is not None and i_val is not None]
        x = np.array([voltage_data[i] for i in valid_indices])
        y = np.array([current_data[i] for i in valid_indices])
        
        # 按电压升序排列
        sort_idx = np.argsort(x)
        x_sorted = x[sort_idx]
        y_sorted = y[sort_idx]
        
        if auto:
            cutoff_voltage = None
        else:
            cutoff_voltage = cutoff_voltages[group_idx - 1]
        try:
            # 使用保形插值（保证单调性）
            pchip = PchipInterpolator(x_sorted, y_sorted)
            
            nums_data = 300
            # 生成密集采样点，用于绘制曲线（扩展范围为原数据范围两侧各扩展10%）
            x_min, x_max = np.min(x_sorted), np.max(x_sorted)
            x_range = x_max - x_min
            x_dense = np.linspace(x_min - 0.1*x_range, x_max + 0.1*x_range, nums_data)
            y_dense = pchip(x_dense)
            
            if auto:
                threshold = (max(y_dense) - min(y_dense)) / (x_max - x_min)  # 导数阈值
                
                for i in range(nums_data):
                    if (y_dense[i + 5] - y_dense[i]) / (x_dense[i + 5] - x_dense[i]) > threshold * 1.5:
                        cutoff_voltage = x_dense[i]
                        break

        except Exception as e:
            print(f"[{group_idx}] 插值失败，使用线性连接: {str(e)}")
            cutoff_voltage = None
            x_dense = x_sorted
            y_dense = y_sorted

        if auto:
            if cutoff_voltage is not None:
                print(f"[{group_idx}] 拐点确定的截止电压: {cutoff_voltage:.3f} V")
            else:
                print(f"[{group_idx}] 未检测到拐点，截止电压无法确定")
            cutoff_voltages.append(cutoff_voltage)

        # 绘制伏安特性曲线
        plt.figure(figsize=(10, 7))
        plt.plot(x_dense, y_dense, 'b-', linewidth=2, label="拟合曲线")
        plt.scatter(x_sorted, y_sorted, c='red', s=50, marker='o',
                    edgecolors='k', label="实验数据")

        if cutoff_voltage is not None:
            plt.axvline(cutoff_voltage, color='red', linestyle='--', label=f"截止电压: {cutoff_voltage:.2f} V")
        
        plt.xlabel(x_label)
        plt.ylabel(y_label)
        plt.title(f"伏安特性曲线，波长 {table_name}")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()

        output_path = os.path.join(output_dir, f"output_{group_idx}.png")
        plt.savefig(output_path, dpi=300)
        plt.close()
        print(f"[{group_idx}] 图像已保存至: {output_path}")

    # 5. 线性拟合：截止电压 vs 频率（计算拟合参数及其不确定度）
    freq_array = np.array(frequencies)
    volt_array = np.array(cutoff_voltages, dtype=float)
    
    if len(freq_array) > 2:
        # 数据点多于两个，正常计算协方差矩阵
        coeffs, cov = np.polyfit(freq_array, volt_array, 1, cov=True)
        err_k = np.sqrt(cov[0, 0])  # 斜率 k 的标准差
    else:
        # 数据点只有两个，直接拟合，不计算协方差矩阵
        coeffs = np.polyfit(freq_array, volt_array, 1)
        err_k = 0.0  # 或者用其他方式表示不确定度为零
    k, b = coeffs

    # 6. 绘制截止电压 vs 频率 图（增加数据点连接线，并设置刻度格式）
    plt.figure(figsize=(8, 6))
    plt.scatter(freq_array, volt_array, color='blue', label="实验数据")
    plt.plot(freq_array, volt_array, 'b-', label="数据连接线")
    plt.plot(freq_array, k * freq_array + b, 'r--', label=f"拟合直线: $U_s = {k:.4f}f + {b:.4f}$\n斜率不确定度: ±{err_k:.4f}")
    plt.xlabel(r"频率 ($10^{-14}$ Hz)")
    plt.ylabel("截止电压 (V)")
    plt.title("截止电压-频率图线")
    plt.grid(True)
    plt.legend()

    ax = plt.gca()
    ax.xaxis.set_major_formatter(FormatStrFormatter('%.3f'))  # 横轴保留三位小数
    ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))  # 纵轴保留两位小数
    plt.tight_layout()

    output_path = os.path.join(output_dir, f"output_{N + 1}.png")
    plt.savefig(output_path, dpi=300)
    plt.close()
    print(f"截止电压-频率图已保存至: {output_path}")

    # 7. 计算普朗克常数
    h0 = 6.64 * 10**(-34)  # 公认普朗克常数 (J·s)
    e = 1.6 * 10**(-19)  # 电子电荷量 (C)
    h = e * (-k) * 1e-14  # 计算出的普朗克常数 (J·s)
    _h = e * err_k * 1e-14  # 计算出的普朗克常数的不确定度 (J·s)
    Ep = (h - h0) * 100 / h0

    # 8. 写入 Excel
    # 将标题写入 A15, A16, A17
    sheet["A15"] = "波长 (nm)"
    sheet["A16"] = "频率 (10^-14 Hz)"
    sheet["A17"] = "截止电压 (V)"

    # 假设数据组数为 5，对应写入 B 到 F 列
    for i, (w, f, u) in enumerate(zip(wavelengths, frequencies, cutoff_voltages)):
        col_letter = chr(ord('B') + i)  # B, C, D, E, F
        sheet[f"{col_letter}15"] = f"{w:.0f}"
        sheet[f"{col_letter}16"] = f"{f:.2f}"
        sheet[f"{col_letter}17"] = f"{u:.2f}"

    # 将计算结果写入 Excel（从 A19 开始）
    sheet["A19"] = "h"
    sheet["B19"] = f"{h * 10 ** 34:.2f} × 10^(-34) J·s"
    sheet["A20"] = "Δh"
    sheet["B20"] = f"{_h * 10 ** 34:.2f} × 10^(-34) J·s"
    sheet["A21"] = "h±Δh"
    sheet["B21"] = f"{h * 10 ** 34:.2f} ± {_h * 10 ** 34:.2f} × 10^(-34) J·s"
    sheet["A22"] = "Ep"
    sheet["B22"] = f"{Ep:.2f}%"

    wb.save(excel_file)
    print("计算结果已写入 Excel")
    print(f"{Ep:.2f}%")

    if save_to_word:
        from docx import Document
    from docx.shared import Inches
    
def Htlanglanglang():
    doc = Document()
    
    # 创建表格（3 行 2 列），确保六张图片能排在一页上
    table = doc.add_table(rows=3, cols=2)
    table.autofit = True  # 自动调整列宽
    
    img_index = 0  # 计数器

    for i in range(1, N + 2):  # 遍历所有图片（包括截止电压-频率图）
        image_path = os.path.join(output_dir, f"output_{i}.png")
        if os.path.exists(image_path):
            row = img_index // 2  # 计算行号
            col = img_index % 2   # 计算列号
            
            cell = table.cell(row, col)  # 获取单元格
            paragraph = cell.paragraphs[0]
            run = paragraph.add_run()
            run.add_picture(image_path, width=Inches(3))  # 调整宽度，确保 3 张图片一行
            
            img_index += 1

    doc.save(word_file)
    print(f"实验结果已保存至 Word: {word_file}")
    
if calculate:
    Htlang()
if save_to_word:
    Htlanglanglang()

 操作步骤

第0步 创建所需文件

在电脑某个位置新建一个文件夹，然后在新建的文件夹里新建一个.xlsx文件、一个.py文件、一个.docx文件以及一个文件夹，完成后大概长下面这样

第1步 填入数据、粘贴代码

打开excel文件，将原始数据输入，格式如下

注意：A列的几个单元格要合并一下不然会报错（我懒的调了），电流的单位中的-11次方这些不需要在excel中就给它改为上标不然不知道会不会报错（我懒的试了）

打开Python文件，将代码粘贴进去，并在代码中注释了控制台的部分对文件路径进行修改

第2步 安装必要的Python扩展库

打开cmd，输入以下代码（二选一）

pip install scipy
pip install python-docx
pip install openpyxl
pip install matplotlib
pip install numpy

pip install openpyxl matplotlib numpy scipy python-docx

第3步 运行代码

一切顺利的话运行结果大致如下

 此外，打开excel，会发现多了一些数据

文件夹里也会多出6张图片

第4步 发现问题并稍加调整

因为找拐点的方法比较粗略，所以可能会存在较大的误差，导致算出来的百分差严重超标，所以我们需要进行手动调整，然后将五个截止电压在代码中事先输入，如下图

第5步 再次运行代码

我们需要不断的调整和尝试，直到把百分差降到10%以下

注意：运行时记得把excel关掉否则会报错

第6步 把图片保存到word中

经过不断的调整，将百分差控制在10以内之后，将代码中的save_to_word这个变量设置为1（当然一开始就设置为1也不是不行，开始教学之前我忘记改了），再次运行代码（此时可以将变量calculate设置为0，不过无所谓，数据量不大，运行也耗不了多长时间），会发现输出多了一句话

然后就可以打开word，按自己的需求调整图片的大小（代码已经帮我们初步调好了），然后去打印店打印结果了。

第7步 定量分析

重新建一个excel，数据填进前4行，把文件路径复制到代码相应位置，同时还要改一下存图片的文件夹地址以及.docx文件的地址（这两个地址如果不存在的话代码会帮你自动建立），不然会将原先保存的东西覆盖，然后把代码注释中控制台第一行的数据组数改成2即可。随后就是同样的运行代码，先得到大概的截止电压，然后自己进行微调，调好后得出h就行。