Solving an integral over gaussian function in spherical coordinates (or an intermediate function involving BesselI)

Question

I am trying to solve an integral from the following function (int) and set of assumptions (as):

int = Exp[(-a)*(r^2 + r1^2 + 2*r*r1*(Cos[f]*Cos[f1] + 
                   Cos[t - t1]*Sin[f]*Sin[f1]))]*Sin[f]*r^2; 
as = Assumptions -> {Element[{a, r, r1, f, f1, t, t1, a}, 
           Reals], a >= 0, r >= 0, r1 >= 0, t >= 0, t <= Pi, 
         t >= 0, t1 <= Pi, f >= 0, f1 <= 2*Pi}; 

The integral is:

Integrate[int, {t, 0, 2*Pi}, {f, 0, Pi}, {r, 0, s}, as]

which seems hopeless to obtain by brute force. I have no problem in obtaining the integral numerically, but it is very costly as I need to compute it many times. Do you have ideas on how to solve it (at least w.r.t. 2 variables would be great and then treat the third variable numerically)? I could obtain analytical expressions w.r.t. one of each of the three variables. E.g.,

intt = Integrate[int, {t, 0, 2*Pi}, as]

yields:

intt = (2*Pi*r^2*BesselI[0, 2*r*r1*a*Sin[f]*Sin[f1]]*Sin[f])/
   E^(a*(r^2 + r1^2 + 2*r*r1*Cos[f]*Cos[f1]))

But then, I get stuck if I try to integrate w.r.t. the remaining variables (including assumptions above does not help):

inttr = Integrate[intt, {r, 0, s}]
inttf = Integrate[intt, {f, 0, Pi}, as]

Answer 1

OP's integral is

ii[s_,r1_,a_,t1_,f1_] := NIntegrate[Exp[(-a)*(r^2+r1^2+2*r*r1*(
   Cos[f]*Cos[f1]+Cos[t-t1]*Sin[f]*Sin[f1]))]*Sin[f]*r^2,{t,0,2*Pi},{f,0,Pi},{r,0,s}]

The short answer is that this is actually independent of t1 and f1 and equal to

iii[s_,r1_,a_] := (Pi (-((E^(-a (r1+s)^2) (-1+E^(4 a r1 s)))/r1)+
   Sqrt[a] Sqrt[Pi] (-Erf[Sqrt[a] (r1-s)]+Erf[Sqrt[a] (r1+s)])))/(2 a^2)

Example:

ii[0.3,0.4,0.7,0.2,0.4]
(* 0.0976496 *)

iii[0.3,0.4,0.7]
(* 0.0976496 *)

To check this claim, first check that ii does not depend on t1, which is clear because of the integral over t. Set t1=0 for convenience. At this point, interpret the combined integral over t and f as an integral over the unit sphere in 3 dimensions, and note that different choices of f1 correspond to different choices of the three coordinate axes (...) hence the integral does not depend on f1. Set f1=0 for convenience. At this point, the integral over t is trivial and simply gives a factor $2\pi$, one can then use Mathematica to explicitly integrate over f, then over r, and this gives the expression in iii.