判断氢原子非零跃迁矩阵元的代码

上一个代码的补充（似乎发不了评论区）

通过数值计算与电偶极、磁偶极、电四极算子有相同表示的算子（如x、yz、Ly（或按照化学书籍常用的写法 Ry）等）的矩阵元，查找非零的跃迁矩阵元。

考虑了E1、E2、M1跃迁，没有考虑自旋等因素。

Dim = Num^3; % 维数

dx = 2*L/(Num-1); % 格点间隔

NormalL = 1/(dx)^3

% 位矢算子

x1 = spdiags(reshape(X,[Dim,1]),0,Dim,Dim);

y1 = spdiags(reshape(Y,[Dim,1]),0,Dim,Dim);

z1 = spdiags(reshape(Z,[Dim,1]),0,Dim,Dim);

% 微分算子和轨道角动量算子

e = ones(Num,1);

Dif1 = spdiags([-e 0*e e],[-1 0 1],Num,Num)/(2*dx);

I1 = speye(Num,Num);

PDx = kron(kron(I1,Dif1),I1);

PDy = kron(kron(I1,I1),Dif1);

PDz = kron(kron(Dif1,I1),I1);

Lz = -1i*(x1*PDy-y1*PDx);

Ly = -1i*(z1*PDx-x1*PDz);

Lx = -1i*(y1*PDz-z1*PDy);

% 态矢量

psi01=reshape(psi1,[Dim,1]);

psi02=reshape(psi2,[Dim,1]);

%检验归一化（波函数显著非零的区域可能因为“轨道半径”比较大而超出边界，造成区域内的总概率明显小于1）

psi01*psi01/NormalL

psi02*psi02/NormalL

% 电偶极检验

psi01*x1*psi02/NormalL

psi01*y1*psi02/NormalL

psi01*z1*psi02/NormalL

% 磁偶极检验

psi01*Lx*psi02/NormalL

psi01*Ly*psi02/NormalL

psi01*Lz*psi02/NormalL

% 电四极检验

psi01*x1*y1*psi02/NormalL

psi01*y1*z1*psi02/NormalL

psi01*x1*z1*psi02/NormalL

psi01*x1*x1*psi02/NormalL

psi01*y1*y1*psi02/NormalL

psi01*z1*z1*psi02/NormalL

% 特别小的数（比如10^-15量级）可以当做0