基于数值积分法利用玻尔-索末菲量子化条件计算一维谐振子的能谱并绘制对应的位移动量(r-p)图。（不能用现有的能级公式直接得出结论）

为了基于数值积分方法和玻尔-索末菲（Bohr-Sommerfeld）量子化条件计算一维谐振子的能谱，我们需要首先了解几个关键概念：

1. **玻尔-索末菲量子化条件**：这是经典力学中的一个近似，用来描述量子粒子在势场中的能量状态。在量子力学中，它对应于角动量取离散值的情况。

2. **一维谐振子**：这是一个简单的物理模型，代表了一个粒子在一个弹性势阱中的运动，其势能函数为 \( V(x) = \frac{1}{2}kx^2 \)，其中 \( k \) 是刚度常数。

3. **数值积分法**：我们将使用数值方法来求解微分方程或积分，因为解析解对于复杂的量子条件可能并不适用。

下面是步骤概述：

### 步骤1: 定义量子化条件
玻尔-索末菲量子化条件要求粒子的行动量 \( p \) 和位置 \( x \) 的乘积在每个能量态上满足整数倍的普朗克常数 \( \hbar \)：
\[ \oint pdx = n\hbar \]

### 步骤2: 将连续的角动量积分转换成离散的能量
你需要对动能 \( T = \frac{p^2}{2m} \) 和势能 \( V(x) \) 的差 \( E = T - V \) 进行积分，然后找到使得积分结果等于整数倍的 \( \hbar \) 的 \( n \) 值。这通常涉及到数值积分（例如辛普森法则、梯形法则等）。

### 步骤3: 编写MATLAB代码
创建一个函数来执行上述积分，比如使用`quadgk`函数来进行高精度数值积分。定义一个范围，比如从 \( -L \) 到 \( L \)，并将步长缩小到足够精确地捕捉量子跃迁。

function energy = bs_quantum_condition(k, hbar, L, N_points)
    m = 1; % 一维谐振子的质量假设为1
    delta_x = L / N_points;
    px = @(x) sqrt(2 * k * (L^2 - x.^2)); % 动量
    energies = zeros(1, N_states); % 初始化存储能级

    for n = 1:N_states % 迭代寻找符合条件的n值
        integral = 0;
        x_values = -L : delta_x : L;
        for i = 1:length(x_values)
            integral = integral + quadgk(@(x) px(x).*dx(x_values(i)), -L, L);
        end
        if abs(integral - n*hbar) < tolerance % 检查是否接近整数倍的hbar
            energies(n) = integral / (2*m); % 计算能量
        else
            break; % 如果找不到，跳出循环
        end
    end
    energy = energies;
end

### 步骤4: 绘制 \( r-p \) 图
为了得到 \( r-p \) 图，你需要先计算粒子的位置和动量分布，然后根据量子化的能级把这些点画出来。由于这里没有具体的 \( r \) 变量（粒子在谐振子上的位置），我们通常讨论的是 \( q-p \) 对应的波函数平方，表示概率密度。

[r, p] = meshgrid(-L, -sqrt(2*k*L), 'uniform'); % 假设r和p均匀分布
psi_squared = exp(-(r.^2 + p.^2)/(2*hbar^2)); % 波函数的幅度

figure
surf(r, p, psi_squared, 'FaceColor', 'interp')
xlabel('Position (q)')
ylabel('Momentum (p)')
title('Probability Density in the (q,p) Plane')

% 根据能量水平，选取对应于不同量子态的ψ的值并在图表中标记
for i = 1:length(energy)
    plot([0, p(end)], [r(end), r(end)], 'k:', 'LineWidth', 1.5);
    hold on
    text(p(end), r(end), num2str(energy(i)), 'HorizontalAlignment', 'right');
end
hold off

### 相关问题--
1. 如何处理边界条件影响量子化结果？
2. 数值积分时如何选择合适的误差容忍度（tolerance）？
3. 如何改进这个算法以提高计算效率？
4. 如何理解\( r \)-\( p \) 图和它与量子态的关系？