I want to numerically integrate $$\iint_P e^{ik(y-y^\prime)^2}e^{ik(x-x^\prime)^2} \,dx\,dy$$
where $P$ is the region defined by a general $N$-gon with an arbitrary number of sides (a polygon has $N=5$).
For simple shapes such as a triangle and square, I can derive the limits/functional limits and then split it into two integrals that can be solved via Simpson's rule or another method. However, for a generalised $N$-gon, this is much harder to derive general limits. I was thinking of using Green's theorem in conjunction with a list of co-ordinates of the vertex and then split it into $N-1$ $1$-dimensional integrals. This, however, requires solving $$\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}= e^{ik(y-y^\prime)^2}e^{ik(x-x^\prime)^2}$$ for $Q$ and $P.$ I have no idea how to do this, and for that matter if it is even possible.
Is there a way to solve this problem (doesn't have to be via Green's theorem) that can be computationally implemented, the simpler, the better.