Legendre-Gauss Quadrature on an Ugly Integral

I'm trying to perform nth-Order Legendre-Gauss Quadrature on an integral of mine. The algorithm I've written for the quadrature works on simple integrals, like the following.

$\int \int \int dxdydz \sin (x) \sin (y) \sin (z)$

However, the problem I'd like to solve is particularly ugly and the results from quadrature don't converge no matter how high in n I go (I've gone to 256th order!). The results are (seemingly) random across several orders of magnitude. I'm wondering if anyone has any insight into this issue and might have ideas on how to attack it.

The integral is used to calculate the intensity in x-ray diffraction of a lattice of atoms. The lattice is perfect except for a collection of dislocations. There are five types of dislocations detailed below, indexed by $\alpha$, which have an associated Burger's vector, $\vec{B}$, and gives a displacement vector on the lattice. Only the z-component of the displacement vector, here defined as $u$ is needed. Note that $\nu$ is a physical parameter (a constant for the lattice) called Poisson's ratio, $a \mbox{ and } c$ are also constants.

Screw Dislocations
Type 1: $\vec{B}^{(1)} = <0 \ 0 \ c>$
Type 2: $\vec{B}^{(2)} = <0 \ 0 \ -c>$
$u_{s} (\vec{r}) = \frac{B^{(\alpha)}_{z}}{2 \pi} \arctan{\frac{y}{x}}$

Edge Dislocations
Type 3: $\vec{B}^{(3)} = < a \ 0 \ 0>$
Type 4: $\vec{B}^{(4)} = < \frac{-a}{2} \ \frac{\sqrt{3}a}{2} \ 0>$
Type 5: $\vec{B}^{(5)} = < \frac{-a}{2} \ \frac{-\sqrt{3}a}{2} \ 0>$
$u_{e} (\vec{r}) = \frac{\nu}{1 - \nu} \frac{\vec{r} \cdot \vec{B}^{(\alpha)}}{2 \pi} \frac{z+2r(\nu - 1)}{r(r - z)}$

Note that $\vec{r} = x \hat{x} + y \hat{y} + z \hat{z} \ , \ r = \sqrt{x^2 + y^2 + z^2}$
Now, the total (z-component of the) displacement vector in the lattice is

$u = \sum_{\alpha = 1}^{5} \sum_{j}^{N^{(\alpha)}} u^{(\alpha)}(x - x_{j}^{(\alpha)}, \ y - y_{j}^{(\alpha)}, \ z)$

$N^{(\alpha)}$ is the number of $\alpha$-type disolcations, and $(x_{j}^{(\alpha)}, \ y_{j}^{(\alpha)})$ is the intersection of the $j$th $\alpha$-type dislocation with the sample surface.

Finally, we can get to the integral. $\vec{h} = h \hat{z}$ and $\vec{q}$ are vectors that are constant in space, $T$ is the sample's thickness, and $S$ is the coherent area, which is simplified as a rectangle bounded by $[x_{0},x_{1},y_{0},y_{1}]$.

$I = \int \int_{S} dxdy \int^{0}_{-T} dz e^{-ihu(\vec{r})}e^{-i \vec{q} \cdot \vec{r}}$

Why doesn't the quadrature work on my integral?!! I've outlined Gauss Quadrature below.

Quadrature
An integral of the function $f(x,y,z)$ can (in theory) be approximated by the following. $\{\theta\}$ are the roots of the nth Legendre polynomial, and $\{w\}$ are the associated weights.

$I = \int_{x_{0}}^{x_{1}} \int_{y_{0}}^{y_{1}} \int{z_{0}}^{z_{1}} dxdydz f(x,y,z)$
$I \approx x_{a}y_{a}z_{a} \sum_{i,j,k}^{n} w_{i}w_{j}w_{k} f(x_{a}\theta_{i} + x_{b}, y_{a}\theta_{j} + y_{b}, z_{a}\theta_{k} + z_{b})$

$x_{a} = \frac{x_{1} - x_{0}}{2} \ , \ x_{b} = \frac{x_{1} + x_{0}}{2} \ , \ \mbox{etc.}$