Improving speed of multidimensional Gaussian quadrature

Question

I am trying to do a numerical integration using Gauss-Legendre quadrature over rectangular regions. The equations used in the code below is given as:

$$ Q_{1}= \int_{c_{i}}^{d_{i}}\int_{a_{i}}^{b_{i}} I_{fn}(x_{o}, y_{o},x_{i},y_{i},x_{p},y_{p},z_{p},z_{o},z_{d},p_D,p_A) \ \mathrm{d} x_{i} \ \mathrm{d} y_{i} \\ \approx \frac{(d_i - c_i)}{2} \frac{(b_i - a_i)}{2} \sum_{u=1}^{n_{i}} \sum_{v=1}^{n_{i}}w_uw_v \ I_{fn}(x_{o}, y_{o},\frac{(a_{i} + b_{i})}{2} + \frac{(b_{i} - a_{i})}{2}x_{u},\frac{(c_{i} + d_{i})}{2} + \frac{(d_{i} - c_{i})}{2}y_{v},x_{p},y_{p},z_{p},z_{o},z_{d},p_D,p_A) \\ Q_{2}= \int_{c_{o}}^{d_{o}}\int_{a_{o}}^{b_{o}} Q_{1}(x_{o}, y_{o},x_{p},y_{p},z_{p},z_{o},z_{d},p_D,p_A) \ \mathrm{d} x_{o} \ \mathrm{d} y_{o} \\ \approx \frac{(d_o - c_o)}{2} \frac{(b_o - a_o)}{2} \sum_{u=1}^{n_{o}} \sum_{v=1}^{n_{o}}w_uw_v \ Q_1(\frac{(a_{o} + b_{o})}{2} + \frac{(b_{o} - a_{o})}{2}x_{u},\frac{(c_{o} + d_{o})}{2} + \frac{(d_{o} - c_{o})}{2}y_{v},x_{p},y_{p},z_{p},z_{o},z_{d},p_D,p_A) $$

My code segment is given here:

Needs["NumericalDifferentialEquationAnalysis`"];
Remove[Result, GaussWeightsRin, GaussWeightsRout, FirstTerm, Ifn, Q1, Q2]
GaussWeightsRin = GaussianQuadratureWeights[15, -1, 1];
GaussWeightsRout = GaussianQuadratureWeights[2, -1, 1];

FirstTerm[xi_, yi_, xp_, yp_, zp_, Zo_, Zd_, pD_, pA_] := 
  Piecewise[{{(Abs[Zo - zp Cos[pA] + Zd Cos[pA - pD] + xp Sin[pA] - xi Sin[pA + ArcTan[-Cos[pD], Sin[pD]]]]^3 (Zd + Zo Cos[pA - pD] - zp Cos[pD] + xp Sin[pD]) (Zd Cos[pA - pD] - xi Sin[pA + ArcTan[-Cos[pD], Sin[pD]]]))/(Sqrt[Abs[Cos[pD]]^2 +  Abs[Sin[pD]]^2] (Abs[yi - yp]^2 +  Abs[zp - Zo Cos[pA] - Zd Cos[pD] + xi Sin[pD]]^2 +  Abs[xp + xi Cos[pD] + Zo Sin[pA] + Zd Sin[pD]]^2)^(3/2) Abs[Zd Cos[pA - pD] - xi Sin[pA + ArcTan[-Cos[pD], Sin[pD]]]]^3 (Zo - zp Cos[pA] + Zd Cos[pA - pD] + xp Sin[pA] - xi Sin[pA + ArcTan[-Cos[pD], Sin[pD]]])), 
            Cos[pD] != 0 || Sin[pD] != 0}}, ComplexInfinity];

Ifn = Compile[{{xo, _Real, 0}, {yo, _Real, 0}, {xi, _Real, 0}, {yi, _Real, 0}, 
               {xp, _Real, 0}, {yp, _Real, 0}, {zp, _Real, 0}, {Zo, _Real, 0}, 
               {Zd, _Real, 0}, {pD, _Real, 0}, {pA, _Real, 0}}, 

               FirstTerm[xi, yi, xp, yp, zp, Zo, Zd, pD, pA] 2 Exp[-(((xo - xi - 0)/5)^2 + ((yo - yi - 0)/5)^2)], 

              RuntimeAttributes -> {Listable}, CompilationTarget -> "C", Parallelization -> True, RuntimeOptions -> "Speed"];

With[{Ifn = Ifn}, 
  Q1 = Compile[{{xo, _Real, 0}, {yo, _Real, 0}, {ai, _Real, 0}, {bi, _Real, 0},
                {ci, _Real, 0}, {di, _Real, 0}, {ni, _Real, 0}, {xp, _Real, 0}, 
                {yp, _Real, 0}, {zp, _Real, 0}, {Zo, _Real, 0}, {Zd, _Real, 0}, 
                {pD, _Real, 0}, {pA, _Real, 0}},

               (* Quadrature Q1 *)
               (di - ci)/2*(bi - ai)/2*
               Sum[Sum[GaussWeightsRin[[u, 2]] GaussWeightsRin[[v, 2]] 
                       Ifn[xo, yo, (ai + bi)/ 2 + (bi - ai)/2 GaussWeightsRin[[u, 1]], (ci + di)/ 2 + (di - ci)/2 GaussWeightsRin[[v, 1]], xp, yp, zp, Zo, Zd, pD, pA], {u, 1, ni}], {v, 1, ni}], 

 RuntimeAttributes -> {Listable}, CompilationTarget -> "C", Parallelization -> True, RuntimeOptions -> "Speed"]];

With[{Q1 = Q1}, 
  Q2 = Compile[{{ma, _Integer, 0}, {mb, _Integer, 0}, {no, _Integer, 0}, {ni, _Integer, 0}, 
                {xp, _Real, 0}, {yp, _Real, 0}, {zp, _Real, 0}, {Zo, _Real, 0}, 
                {Zd, _Real, 0}, {pD, _Real, 0}, {pA, _Real, 0}}, 

               Module[{ao = ma*0.048046875 (1 - 1/(2 Abs[ma])) - 0.048/2, 
                       bo = ma*0.048046875 (1 - 1/(2 Abs[ma])) + 0.048/2, 
                       co = mb*0.048046875 (1 - 1/(2 Abs[mb])) - 0.048/2, 
                       do = mb*0.048046875 (1 - 1/(2 Abs[mb])) + 0.048/2, 
                       ai = ma*0.048046875 (1 - 1/(2 Abs[ma])) - 8, 
                       bi = ma*0.048046875 (1 - 1/(2 Abs[ma])) + 8, 
                       ci = mb*0.048046875 (1 - 1/(2 Abs[mb])) - 8, 
                       di = mb*0.048046875 (1 - 1/(2 Abs[mb])) + 8 },

               (* Quadrature Q2 *)
               (do - co)/2*(bo - ao)/2*
               Sum[Sum[GaussWeightsRout[[u, 2]] GaussWeightsRout[[v, 2]] Q1[(ao + bo)/2 + (bo - ao)/2 GaussWeightsRout[[u, 1]], (co + do)/2 + (do - co)/2 GaussWeightsRout[[v, 1]], ai, bi, ci, di, ni, xp, yp, zp, Zo, Zd, pD, pA], {u, 1, no}], {v, 1, no}]], 

               RuntimeAttributes -> {Listable}, CompilationTarget -> "C", Parallelization -> True, RuntimeOptions -> "Speed"]];

With[{Q2 = Q2}, 
  Result = Compile[{{M, _Integer, 1}, {no, _Integer, 0}, {ni, _Integer, 0}, {xp, _Real, 0}, {yp, _Real, 0}, {zp, _Real, 0}, {Zo, _Real, 0}, {Zd, _Real, 0}, {pD, _Real, 0}, {pA, _Real, 0}}, 

                  Q2[Compile`GetElement[M, 1], Compile`GetElement[M, 2], no, ni, xp, yp, zp, Zo, Zd, pD, pA],

                  RuntimeAttributes -> {Listable}, CompilationTarget -> "C", Parallelization -> True, RuntimeOptions -> "Speed"]];

I want to run this code (Result[]) for a set of different points {ma,mb} and for that, I use Outer.

Num = 10;  
CordsX = DeleteCases[Range[Num/2, -Num/2, -1], 0, Infinity];
CordsY = {1};
Outer[Result[{#2, #1}, 2, 15, 10, 10, 10, 100, 200, \[Pi], 0] &, CordsY, CordsX, 1] // AbsoluteTiming

This is very slow for my application as I need to run this over more points (i.e. for higher values of Num). Have I correctly used compile in my code? and how do I speed up my code?.