I'm trying to numerically integrate the function $$I(\mu) =\exp{(-\tau/\beta)} I_0[1-u_1 (1 - \mu)-u_2 (1-\mu)^2], $$ $I_0$, $u_1$ , $u_2$, $\tau$ and $\beta$ are known, and $$\mu = \sqrt{1-(x^2+z^2)}$$
Now it comes the difficult part. The nasty non-rectangular region where I want the integral is: $$ \left\{\begin{array}{ll} R_{in} ^2 < \left(\frac{x-X_p-(z-Z_p)\cdot\sin \phi}{\cos \theta}\right)^2 + (\frac{z-Z_p}{\sin\theta})^2<R_{out}^2 \\ x^2+z^2<1 \\ (x-X_p)^2+(z-Z_p)^2>R_p^2 \end{array}\right. $$
$R_p<R_{in}<R_{out}$, $R_{p}<1$, $X_p$ and $Z_p$ all are known.
I tried to use software to numerically solve it, such as MATLAB, Python or Mathematica. But the problem I have is that I can't parametrize this region in order to be able to integrate it. So, what can I use to numerically solve it? Do I need to parametrize it?.
Pd: another thing I tried is integrate the function in a piecewise function (where the function is 0 outside the region), using the trapeze method. But it doesn't work since I have to deal with discontinuities. So, what should I do?