the four Whittaker functions reads [@gradstejn_table_2009] $$\begin{aligned}
    \int
    &M_{a_1+b_1,b_1-1/2}(q_1 r) M_{a_2+b_2,b_2-1/2}(q_2 r)
      M_{a_3+b_3,b_3-1/2}(q_3 r') M_{a_4+b_4,b_4-1/2}(q_4 r')
      \frac{r_<^l}{r_>^{l+1}} dr dr' \nonumber
    \\
    &=\sum_{i_1=0}^{a_1} \sum_{i_2=0}^{a_2} \sum_{i_3=0}^{a_3}
      \sum_{i_4=0}^{a_4} T_{{\boldsymbol}{a},{\boldsymbol}{b},{\boldsymbol}{q}}({\boldsymbol}{i})
      \left(f_{i_1+i_2+b_1+b_2+l+1}^{i_3+i_4+b_3+b_4-l}
      \left(\frac{q_3+q_4}{2},\frac{q_1+q_2}{2}\right)
      + f_{i_3+i_4+b_3+b_4+l+1}^{i_1+i_2+b_1+b_2-l}
      \left(\frac{q_1+q_2}{2},\frac{q_3+q_4}{2}\right)\right),
  \end{aligned}$$ where $$\begin{aligned}
    T_{{\boldsymbol}{a},{\boldsymbol}{b},{\boldsymbol}{q}}({\boldsymbol}{i}) =
    \prod_{k=1}^4{{a_k}\choose{i_k}}
    \frac{\Gamma(2b_k)}{\Gamma(2b_k+i_k)}
    (-1)^{i_k}q_k^{b_k+i_k}
  \end{aligned}$$ and we have made use of $$\begin{aligned}
    f_{a}^{b}(x,y) = \int_0^\infty \int_r^\infty e^{-\lambda r -
    \lambda' r'} r^{a-1}
    {r'}^{b-1} dr' dr=\frac{\Gamma(a+b)}{a{\lambda'}^{a+b}}
    {_2F_1}\left(a, a+b, a + 1, -\frac{\lambda}{\lambda'}\right).
  \end{aligned}$$ The properties of $f_{a}^{b}(x,y)$ functions are described in some detail in [@dzikowski_generating_2019]. The bold ${\boldsymbol}a$, ${\boldsymbol}b$, ${\boldsymbol}q$ and ${\boldsymbol}i$ are lists of four values, i.e., ${\boldsymbol}a = \{a_{1}, a_{2}, a_{3}, a_{4}\}$ with similar expressions for ${\boldsymbol}b$, ${\boldsymbol}q$ and ${\boldsymbol}i$.

Now we would like to compute the scattering factors. For this we consider the electronic density, which is averaged over the angular variables. By using the orthogonality of spherical spinors the angular dependence integrates out trivially. Therefore, we get the radial density as $$\begin{aligned}
    r^2 \rho_{n_r,\kappa} (r,Z^*) =
    |G_{n_r,\kappa}(r,Z^*)|^2 + |F_{n_r,\kappa}(r,Z^*)|^2, \label{eq:39}
  \end{aligned}$$ which gives the scattering factors as $$\begin{aligned}
    f_{n_r,\kappa}(q,Z^*) = \int \left(|G_{n_r,\kappa}(r,Z^*)|^2 +
    |F_{n_r,\kappa}(r,Z^*)|^2\right) e^{i {{\boldsymbol}q\cdot{\boldsymbol}r}}
    d{\boldsymbol}{r}. \label{eq:40}
  \end{aligned}$$ Expanding Whittaker functions in a finite series in Eq. (\[eq:30\]) and using [@gradstejn_table_2009] $$\begin{aligned}
    \int e^{-\alpha r}r^{n-2} e^{i {\boldsymbol}{q} \cdot {\boldsymbol}{r}} d{\boldsymbol}{r} =
    