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

  • $\begingroup$ Scipy has the confluent hypergeometric function, is that the main thing you need? $\endgroup$ Commented Oct 10, 2022 at 15:58
  • $\begingroup$ Unfortunately Scipy's function requires real a,b parameters $\endgroup$
    – mmikkelsen
    Commented Oct 10, 2022 at 18:53
  • $\begingroup$ Ah gotcha, sympy has this function which could also work. If efficiency and compatibility with numpy / scipy is important then you can probably lambdify it. $\endgroup$ Commented Oct 10, 2022 at 19:01
  • $\begingroup$ That might work! I will add some context to the question to better explain why I need numpy / scipy $\endgroup$
    – mmikkelsen
    Commented Oct 10, 2022 at 19:17

1 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])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.