johnson.py 文件源码

def amplitude(wf, R, edash, mu):
    # Mies    F ~ JA + NB       J ~ sin(kR)/kR
    # normalization sqrt(2 mu/pu hbar^2) = zz
    zz = np.sqrt(2*mu/const.pi)/const.hbar

    oo, n, nopen = wf.shape

    # two asymptotic points on wavefunction wf[:, j]
    i1 = oo-5
    i2 = i1-1
    x1 = R[i1]*1.0e-10
    x2 = R[i2]*1.0e-10

    A = np.zeros((nopen, nopen))
    B = np.zeros((nopen, nopen))
    oc = 0
    for j in range(n):
        if edash[j] < 0:
            continue
        # open channel
        ke = np.sqrt(2*mu*edash[j]*const.e)/const.hbar
        rtk = np.sqrt(ke)
        kex1 = ke*x1
        kex2 = ke*x2

        j1 = sph_jn(0, kex1)[0]*x1*rtk*zz
        y1 = sph_yn(0, kex1)[0]*x1*rtk*zz

        j2 = sph_jn(0, kex2)[0]*x2*rtk*zz
        y2 = sph_yn(0, kex2)[0]*x2*rtk*zz

        det = j1*y2 - j2*y1

        for k in range(nopen):
            A[oc, k] = (y2*wf[i1, j, k] - y1*wf[i2, j, k])/det
            B[oc, k] = (j1*wf[i2, j, k] - j2*wf[i1, j, k])/det

        oc += 1

    AI = linalg.inv(A)
    K = B @ AI

    return K, AI, B