I want to solve the following expression (used to obtain an analytic solution to a current distribution inside a workpiece):
$$a_{mn} = -\frac{\frac{4}{ab} \int_0^a \int_0^b f(x',y')\sin(px')\sin(qy')\mathrm{d}x'\mathrm{d}y'}{t\sinh(tc)}$$
Here. $a$,$b$ are scalar constants and $p = \frac{m\pi}{a}$, $q = \frac{n\pi}{b}$, $t=\sqrt{p^2+q^2}$.
The function $f$ is a Gaussian distribution: $f(x,y) = \frac{I_0d}{\pi\sigma^2} \exp(-\frac{r^2d}{\sigma^2})$
I was wondering if using the Simpson's method is the smartest way to solve the double integral. Are there alternative solutions which would be more computationally efficient?
Additional information: The resulting matrix $a_{mn}$ is used in the following infinite series to obtain the final result: $\sum_{m=1}^\infty \sum_{n=1}^\infty a_{mn}p\cos(px)\sin(qy)\cosh(tz)$