rkr.py 文件源码

def turning_points(mu, vv, Gv, Bv, dv=0.1):
    DD = np.sqrt(const.h/(8*mu*const.m_u*const.c*100))*1.0e10/np.pi
    # Gv spline
    gsp = splrep(vv, Gv, s=0)
    # Bv spline
    bsp = splrep(vv, Bv, s=0)

    # vibrational QN at which to evaluate turning points
    V = np.arange(dv-1/2, vv[-1], dv) 
    # add a point close to v=-0.5, the bottom of the well
    V = np.append([-1/2+0.0001], V)
    Rmin = []; Rmax = []; E = []
    # compute turning points using RKR method
    print(u"RKR: v   Rmin(A)  Rmax(A)  E(cm-1)")
    for vib in V:
        E.append(G(vib, gsp))    # energy of vibrational level
        ff = fg_integral(vib, gsp, bsp, lambda x, y: 1)
        gg = fg_integral(vib, gsp, bsp, B)
        fg = np.sqrt(ff/gg + ff**2)
        Rmin.append((fg - ff)*DD) # turning points
        Rmax.append((fg + ff)*DD)
        frac, integ = np.modf(vib)
        if frac > 0 and frac < dv: 
            print(u"     {:d}   {:6.4f}    {:6.4f}    {:6.2f}".
                  format(int(vib), Rmin[-1], Rmax[-1], np.float(E[-1])))

    return Rmin, Rmax, E