計算 Rabi 振蕩的振蕩頻率需要知道偶極躍遷矩陣元，即形如

的矩陣元，分別對應(yīng) ?的情形。而這可以用CG系數(shù)表示為（櫻井純的課后題）

其中 R 代表徑向波函數(shù)。如果引入等效的一維本征波函數(shù)

那么就是

可以方便地用 matlab 數(shù)值計算第一個括號內(nèi)的積分，

后面部分查CG系數(shù)（）的表再帶入就行了

至于含時勢能情形的Schrodinger方程數(shù)值求解，只用在之前的代碼上稍作修改：

1. 修改 UMmt

UMmt = exp(-1i*dT*(Kx.^2+Ky.^2+Kz.^2)/4); % exp(-i*(-0.5*Laplacian)*dT/2), 處于動量空間

2. 球?qū)ΨQ勢能 Up 固定，總勢能修改為

Zr = R.*cos(Tt);

Xr = R.*sin(Tt).*cos(P);

Yr = R.*sin(Tt).*sin(P);

F0 = 1.5;%按需修改

Om = 6.25; %按需修改，應(yīng)當(dāng)是能級間距+失諧

Utm = @(z,t) F0*sin(Om*t)*z; % q=0

% Utm = @(x,y,t) F0*(sin(Om*t)*x+cos(Om*t)*y)/sqrt(2); % |q|=1

Unow = Up+ Utm(Zr,0);

% Unow = Up+ Utm(Xr,Yr,0);

3. for 循環(huán)中時間演化的部分修改為

Unow = Up+Utm(Zr,t*dT);

% Unow = Up+Utm(Xr,Yr,t*dT);

UVh = exp(-1i*Unow*dT);

%exp(-i*(-0.5*Laplacian)*dT/2)

psik = UMmt.*fftshift(fftn(psi3));

psi3 = ifftn(ifftshift(psik));

% exp(-i*U*dT)

psi3 = UVh.*psi3;

%exp(-i*(-0.5*Laplacian)*dT/2)

psik = UMmt.*fftshift(fftn(psi3));

psi3 = ifftn(ifftshift(psik));

至于Rabi振蕩的角頻率，零失諧時就是

M 就是之前提到的矩陣元。如果還有非零的失諧，那就是 計算失諧與上式的平方和，再開根號。

矩陣元 M 中積分計算的一個例子：

N1 = 4000;

x1 = linspace(1e-3,L,N1)';

l = 0;?l2 = 1;?

m = 0;?m2 = 0;?

nl1 = 1;nl2 = 1;

[En1,Vl] = RadEig(x1,SphP(x1),l);

[En2,Vl2] = RadEig(x1,SphP(x1),l2);

I=Vl2(:,nl2)'*((x1).*Vl(:,nl1))

I 就是想要的積分（的近似值）了。

自定義函數(shù) RadEig 見上一個文章