I am trying to find the definite integral of $$I=\int_{0}^{\infty}dx\left(\frac{\sqrt{n^{2l-1}}M_{n,l+\frac{1}{2}}(\frac{2x}{n})}{\Gamma(2l+2)}\right)x\left(\frac{(2r)^{l_f+1} {}_0F_1(2l_f+2;-2r)}{\sqrt{2}\Gamma(2l_f+2)}\right)$$.
Where $M_{n,l+\frac{1}{2}}(\frac{2x}{n})$ is the WhittakerM function, ${}_0F_1(2l_f+1;-2r)$ is the Hypergeometric0F1 function, $r=a_0*x$, where $a_0$ is the bohr radius in Rydberg units ($a_0=\frac{h^2}{4\pi^2m_ee_0^2}$). I am trying to compute this integral for when $n=2,l=0$ and $l_f=1$ and my input is this:
h=UnitConvert[Quantity["PlanckConstant"]];
m=UnitConvert[Quantity["ElectronMass"]];
q=UnitConvert[Quantity["ElectronCharge"]];
a=h^2/(4*Pi^2*m*q^2)
Integrate[
x (Sqrt[2^(2*0 - 1)]*WhittakerM[2, 0 + 1/2, (2*x)/2])/
Gamma[2*0 + 2] ((2*a*x)^(1 + 1) Hypergeometric0F1[2*1 + 2, -2*a*x])/(
Sqrt[2] Gamma[2*1 + 2]), {x, 0, \[Infinity]}]
But I get
Integrate:Integrate was unable to determine the units of quantities that appear
in the input"
And if I omit the bohr radius the output I get is the input, any ideas on what's going wrong?