Can this integral involving Gaussian and sine functions be written in terms of special functions? Recently, I encountered the following integral. Given that it is relatively of a simple form, I suspect it is written in terms of known special functions. I know for $\phi = \pi/2$, it corresponds to Craig formula. 
$$
I(a, \phi) = \int_{0}^{\phi}  \exp \left[ -\frac{a^2}{(\sin{\theta})^2}    \right] d\theta
$$
Can anyone enlighten me regarding this integral?
 A: It is embarrassing to write my own answer to my question, but I do so for future readers seeking for a solution of similar integrals. 
The integral $I(a,\theta)$, defined in the question turns out to be written in terms of error functions and Owen's T function, (i.e., it implies that the integral has something to do with bivariate normal distribution). 
The answer is 
$$
I(a,\theta) = 2 \pi  T\left(\sqrt{2} a \cot (\theta ),\tan (\theta )\right)+\frac{1}{2}
   \pi  \text{erf}(a) (\text{erf}(a \cot (\theta ))-1),
$$
where $T$ is Owen's T function and $a$ is assumed to be positive. 
Edit: 
Let me explain how it was derived.
This identity is evaluated first by calculating $\frac{\partial I}{\partial a}$ which is 
$$
\int -2 a e^{-a^2 \csc ^2(\theta )} \csc ^2(\theta ) \, d\theta = \sqrt{\pi } e^{-a^2} \text{erf}(a \cot (\theta )),
$$
by making the change of variable $\cot \theta = y$. 
Next, $I$ is recovered by integrating the RHS of the above equation with regard to $a$ using the fact that this integral amounts to calculating the CDF of skew-normal distribution](https://en.wikipedia.org/wiki/Skew_normal_distribution).
