def solveCSE(en, rot, mu, R, VT):
n, m, oo = VT.shape
V = np.transpose(VT)
# find channels that are open, as defined by E' > 0
edash = en - np.diag(VT[:, :, -1])
openchann = edash > 0
nopen = edash[openchann].size
mx = matching_point(en, rot, V, R, mu)
if mx < oo-5:
out = leastsq(eigen, (en, ), args=(rot, mx, V, R, mu), xtol=0.01)
en = float(out[0])
# solve CSE according to Johnson renormalized Numerov method
WI = WImat(en, rot, V, R, mu)
RI = RImat(WI, mx)
wf = []
if nopen > 0:
oc = 0
for j, ed in enumerate(edash):
if ed > 0:
f = fmat(j, RI, WI, mx)
wf.append(wavefunction(WI, oc, f))
oc += 1
else:
f = fmat(0, RI, WI, mx)
wf.append(wavefunction(WI, nopen, f))
wf = np.array(wf)
wf = np.transpose(wf)
if nopen == 0:
wf = normalize(wf, R)
else:
K, AI, B = amplitude(wf, R, edash, mu) # shape = nopen x nopen
# K = BA-1 = U tan xi UT
eig, U = linalg.eig(K)
# form A^-1 U cos xi exp(i xi) UT
I = np.identity(nopen, dtype=complex)
xi = np.arctan(eig)*I
cosxi = np.cos(xi)*I
expxi = np.exp((0+1j)*xi)*I
expxiUT = expxi @ np.transpose(U)
cosxiexpxiUT = cosxi @ expxiUT
UcosxiexpxiUT = U @ cosxiexpxiUT
Norm = AI @ UcosxiexpxiUT # == (cu, su) complex
# complex wavefunction array oo x n x nopen
wf = wf @ Norm
return wf, en
评论列表
文章目录