# (Regular) Coulomb wave function

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

• Scipy has the confluent hypergeometric function, is that the main thing you need? Oct 10, 2022 at 15:58
• Unfortunately Scipy's function requires real a,b parameters Oct 10, 2022 at 18:53
• 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. Oct 10, 2022 at 19:01
• That might work! I will add some context to the question to better explain why I need numpy / scipy Oct 10, 2022 at 19:17

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 + 1j*imag_integral, 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)))