(Regular) Coulomb wave function

Question

I'm looking for a way to implement the regular Coulomb wave function in python. This function is a solution to

\begin{align} \frac{\text{d}^2\,u}{\text{d}z^2}+\left(1-\frac{2\eta}{z}-\frac{\ell(\ell+1)}{z^2} \right)u(z)=0 \end{align}

The regular Coulomb wave function is given by

\begin{align} F_\ell(\eta,z) = C_\ell(\eta)z^{l+1}e^{-iz} \mathstrut_1 F_1(\ell+1-i\eta,2\ell+2,2iz), \end{align}

which is the function I'm interested in implementing in Python. In a physics context this is Schödinger's equation with a Coulomb potential; $z$ is the radius, $\ell$ is a quantum number and $\eta=Zm\alpha/(\hbar k)$. More specifically I'm interested in implementing the regular Coulomb function for the repulsive interaction with $\ell=1$.

I've looked at mpmath but I'm not sure how to use the mpf class with np.array and scipy's curve_fit

Answer 1

I ended up using the integral representation of the regular Coulomb wave function given by \begin{align} F_\ell(\eta,\rho) = \frac{\rho^{\ell+1}2^\ell e^{i\rho-(\pi\eta/2)}}{|\Gamma(\ell+1+i\eta)|} \int_0^1 e^{-2i\rho t}t^{\ell+i\eta}(1-t)^{\ell-i\eta} \, \text{d}t \end{align} Something like this

    def complex_quadrature(func, a, b, **kwargs):
        def real_func(x):
            return np.real(func(x))
        def imag_func(x):
            return np.imag(func(x))
        real_integral = quad(real_func, a, b, **kwargs)
        imag_integral = quad(imag_func, a, b, **kwargs)
        return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

    def RegularCoulomb(l,eta,rho):
        First = rho**(l+1)*2**l*np.exp(1j*rho-(np.pi*eta/2),dtype='complex_')/(abs(gamma(l+1+1j*eta)))
        integral = complex_quadrature(lambda t: np.exp(-2*1j*rho*t,dtype='complex_')*t**(l+1j*eta)*(1-t)**(l-1j*eta),0,1)[0]
        return np.array(First*integral,dtype='complex_')

    def C(l,eta):
        return 2**l*np.exp(-np.pi*eta/2)*(abs(gamma(l+1+1j*eta))/(factorial(2*l+1)))
    #Compare to mpmath
    xes = np.linspace(0,10,100)
    coulombfwave = [RegularCoulomb(1,-2,i) for i in xes]

    F1 = lambda x: mpmath.coulombf(1,-2,x)
    mpmath.plot([F1], [0,10])
    plt.plot(xes,coulombfwave)