How to calculate this complementary Bessel function? I am trying to calculate this complementary Bessel function
$$\Psi(a,b,\gamma)=\int_0^\infty\Phi({a\over \sqrt{u}}+b\sqrt{u}){u^{\gamma-1}e^{-u}\over \Gamma(\gamma)}du$$
where $\Phi$ is the standard normal distribution function and $\Gamma$ is the gamma function.
I am trying to find a way to calculate a result based upon this equation.  From my limited understanding of mathematics, this is impossible to do but can be approximated possibly with an improper integral.
Please show me how to calculate this complementary Bessel function.  I don't think I could even realistically solve this over the next decade.
 A: It seems that you are somehow referring to integral representation of the Macdonald function $K_{\gamma}(r)$, also called modified Bessel function of the 2nd kind:
\begin{align}
K_{\gamma}(r)&=\int_0^{\infty}e^{-r\cosh\theta}\cosh\gamma\theta\,d\theta=\\
&=\frac12\int_{-\infty}^{\infty}e^{-r\cosh\theta}\cosh\gamma\theta\,d\theta=\\
&=\frac12\int_{-\infty}^{\infty}e^{-r\cosh\theta}e^{\gamma\theta}\,d\theta=\\
&=\frac12\int_0^{\infty}e^{-r(v+v^{-1})/2}v^{\gamma-1}dv.\tag{1}
\end{align}
(At the last step, I made the change of variables $v=e^{\theta}$.)
To reduce your integral to (1), it suffices to make a linear change of variables:
\begin{align}
\Psi(a,b,\gamma)&=\frac{1}{\Gamma(\gamma)}\int_0^{\infty}\frac{\exp\left\{-\frac12\left(\frac{a}{\sqrt{u}}+b\sqrt{u}\right)^2\right\}}{\sqrt{2\pi}}u^{\gamma-1}e^{-u}du=\\
&=\frac{e^{-ab}}{\Gamma(\gamma)\sqrt{2\pi}}\int_{0}^{\infty}e^{-\left((b^2+2)u+a^2u^{-1}\right)/2}u^{\gamma-1}du=\\
&=\frac{e^{-ab}}{\Gamma(\gamma)\sqrt{2\pi}}\left(\frac{a}{\sqrt{b^2+2}}\right)^{\gamma}\int_{0}^{\infty}e^{-a\sqrt{b^2+2}\left(v+v^{-1}\right)/2}v^{\gamma-1}dv
=\\
&=\frac{2e^{-ab}}{\Gamma(\gamma)\sqrt{2\pi}}\left(\frac{a}{\sqrt{b^2+2}}\right)^{\gamma}
K_{\gamma}\left(a\sqrt{b^2+2}\right).\tag{2}
\end{align}
The continuation really depends on what you want to do. One can e.g. write $K_{\gamma}$ as a linear combination of the modified Bessel functions $I_{\pm\gamma}$ of the 1st kind, which have explicit series representations. This can be convenient for small arguments. For large arguments, there are known asymptotic expansions. I can develop the details if you explain what you need.

Added: (in the process of discussion with OP the answer for the cumulative density function was found in Madan et al, European Finance Review 2, (1998)).
Introduce
\begin{align}
\Phi(x)&=\int_{-\infty}^{x}\frac{e^{-t^2/2}}{\sqrt{2\pi}}dx=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{x}{\sqrt{2}}\right)\right],\\
\Psi(a,b,\gamma)&=\frac{1}{\Gamma(\gamma)}\int_0^{\infty}\Phi\left(\frac{a}{\sqrt{u}}+b\sqrt{u}\right)u^{\gamma-1}e^{-u}du.
\end{align}
Also define $c=|a|\sqrt{2+b^2}$, $d=\displaystyle \frac{b}{\sqrt{2+b^2}}$, $s_a=\mathrm{sign}(a)$ 
  and introduce degenerate hypergeometric function of two variables:
  $$F(\alpha,\beta,\gamma;x,y)=\frac{\Gamma(\gamma)}{\Gamma(\alpha)\Gamma(\gamma-\alpha)}
  \int_0^1\frac{u^{\alpha-1}(1-u)^{\gamma-\alpha-1}}{(1-u x)^{\beta}}e^{uy}du.
  $$
Then, according to the mentioned paper (+ confirmed by numerics):
  \begin{align}
  \Psi(a,b,\gamma)=\frac{c^{\gamma+\frac12}e^{s_a c}(1+d)^{\gamma}}{\sqrt{2\pi}\Gamma(\gamma)}\biggl[\frac{1}{\gamma}K_{\gamma+\frac12}(c)
  F\left(\gamma,1-\gamma,1+\gamma;\frac{1+d}{2},-s_a c(1+u)\right)&\\
  -\frac{s_a(1+d)}{1+\gamma}K_{\gamma-\frac12}(c)
  F\left(1+\gamma,1-\gamma,2+\gamma;\frac{1+d}{2},-s_a c(1+u)\right)&\\
+  \frac{s_a}{\gamma}K_{\gamma-\frac12}(c)
  F\left(\gamma,1-\gamma,1+\gamma;\frac{1+d}{2},-s_a c(1+u)\right)\biggr].
  \end{align}
