Does an an analytic solution exist to these integrals? For a quantum mechanical simulation
I have these 3 functions of euclidian position (x,y) and constants (p,q,r,s,t,u):
$A1(x,y,p,q)=-\cos \left(\frac{\sqrt{3} \pi  x (p+q)}{A}\right) \sin \left(\frac{\pi  y (p-q)}{A}\right)+\cos \left(\frac{\sqrt{3} \pi  q x}{A}\right) \sin \left(\frac{\pi  y (2
   p+q)}{A}\right)-\cos \left(\frac{\sqrt{3} \pi  p x}{A}\right) \sin \left(\frac{\pi  y (p+2 q)}{A}\right)$
$A2(x,y,r,s)=\sin \left(\frac{\sqrt{3} \pi  x (r+s)}{A}\right) \sin \left(\frac{\pi  y (r-s)}{A}\right)+\sin \left(\frac{\sqrt{3} \pi  s x}{A}\right) \sin \left(\frac{\pi  y (2
   r+s)}{A}\right)-\sin \left(\frac{\sqrt{3} \pi  r x}{A}\right) \sin \left(\frac{\pi  y (r+2 s)}{A}\right)$
Where
$q=0,1,2,3... \quad p=q+1,q+2,q+3...$
$s=0,1,2,3..., \quad r=s+1,s+2,s+3...$
and 
$E(x,y,t,u)=A2(x,y,t,u)+iA1(x,y,t,u)$
i.e.
$E(x,y,t,u)=\sin \left(\frac{\sqrt{3} \pi  x (t+u)}{A}\right) \sin \left(\frac{\pi  y (t-u)}{A}\right)+\sin \left(\frac{\sqrt{3} \pi  u x}{A}\right) \sin \left(\frac{\pi  y (2
   t+u)}{A}\right)-\sin \left(\frac{\sqrt{3} \pi  t x}{A}\right) \sin \left(\frac{\pi  y (t+2 u)}{A}\right)+i \Biggl(-\cos \left(\frac{\sqrt{3} \pi  x (t+u)}{A}\right) \sin
   \left(\frac{\pi  y (t-u)}{A}\right)+\cos \left(\frac{\sqrt{3} \pi  u x}{A}\right) \sin \left(\frac{\pi  y (2 t+u)}{A}\right)-\cos \left(\frac{\sqrt{3} \pi  t x}{A}\right) \sin
   \left(\frac{\pi  y (t+2 u)}{A}\right)\Biggr)$
Where
$u=0,\frac{1}{3},\frac{2}{3},\frac{3}{3}...\quad t=q+1,q+2,q+3...$
I wish to calculate the integrals of products of pairs of different combinations of these functions in equilateral triangle shaped regions of the plane, i.e. I want to find expressions for
$\int\int A1(x,y,p,q)A1(x,y,r,s) dxdy$ 
$\int\int A1(x,y,p,q)A2(x,y,r,s) dxdy$
$\int\int A2(x,y,p,q)A2(x,y,r,s) dxdy$
$\int\int A1(x,y,p,q)E(x,y,t,u) dxdy$
$\int\int A2(x,y,p,q)E(x,y,t,u) dxdy$
$\int \int E(x,y,t,u)E(x,y,v,w) dxdy$     
Where the integral is over the surface of an equilateral triangle. I have tried transforming the problem, like this:
https://math.stackexchange.com/a/955188/441529
But have not yet managed to find an analytic solution to the integral.
 A: If I am right, all the terms are of the form
$$\int_D\text{ soc}(ax)\text{ soc}(by)\text{ soc}(cx)\text{ soc}(dy)\,dx\,dy$$
where $\text{soc}$ denotes a sine or a cosine, and the integration domain is a triangle.
By repeated application of the product-to-sum formulas, you can turn that to a sum of
$$\int_D \text{ soc}(a'x+b'y)\text{ soc}(c'x+d'y)\,dx\,dy,$$
then
$$\int_D \text{ soc}(a''x+b''y)\,dx\,dy.$$
You can decompose the domain as an algebraic sum of trapezoids and evaluate
$$\int_{x=x_0}^{x_1}\int_{y=0}^{px+q}\text{ soc}(a''x+b''y)\,dx\,dy.$$
The integration on $y$ will yield terms in $\text{soc}(a'''x+b''')$, so, yes, there is an analytic solution.
A: Notice functions $A_1, A_2$ and $E$ are linear combinations of traveling waves $e^{ip\cdot r} = e^{i(p_1 x + p_2 y)}$. So do integrals of their product. All integrals you mentioned are linear combinations of the form 
$$\int_\Delta e^{i p\cdot r} dxdy$$
where $\Delta$ is your triangle.
Instead of triangle, consider generalization of this sort of integral over
any simple polygon.
Let $\Omega$ be any simple polygonal region with vertices $v_1 = (x_1,y_1), v_2  = (x_2,y_2), \ldots, v_{n} = (x_n,y_n)$ ordered in counter-clockwise orientation. For convenience, let $v_{n+1} = v_1$. 
Let $p = (p_1,p_2)$ be any non-zero vector such that all $p \cdot v_k \ne p \cdot v_{k+1}$.
Embed the plane $\mathbb{R}^2$ into $\mathbb{R}^3$, lift the polygon $\Omega$
to a polygonal prism $\Omega \times [0,1]$ in $\mathbb{R}^3$ and apply divergence theorem, we have
$$\begin{align}
\int_\Omega e^{ip\cdot r} dxdy
&= \int_{\Omega \times [0,1]} e^{ip\cdot r}dxdydz
= -\frac{i}{|p|^2} \int_{\Omega \times [0,1]} \nabla \cdot (p e^{ip\cdot r}) dxdydz\\
&= -\frac{i}{|p|^2} \int_{\partial(\Omega \times [0,1])} e^{ip\cdot r} p \cdot dS
= -\frac{i}{|p|^2} \int_{\partial\Omega} e^{ip\cdot r} p\cdot(ds \times \hat{z})\\
&= -\frac{i}{|p|^2} \int_{\partial\Omega} e^{ip\cdot r} \hat{z}\cdot (p \times ds)
\end{align}
$$
where $dS$ is the surface element of $\partial(\Omega \times [0,1])$ and $ds$ is the line element for $\partial\Omega$.  
Break $\partial\Omega$ into $n$ line segments $v_iv_{i+1}$, integrate over each line segment and sum, we obtain
$$
\int_\Omega e^{ip\cdot r} dxdy
= -\frac{1}{|p|^2}\sum_{k=1}^n 
(e^{ip\cdot v_{k+1}} - e^{ip\cdot v_k})\frac{ \hat{z}\cdot( p \times (v_{k+1} - v_k))}{p\cdot(v_{k+1}-v_k)}
\tag{*1}
$$
Expanding all factors out, RHS becomes following mess:
$$
-\frac{1}{p_1^2+p_2^2}\sum_{k=1}^n 
\big(e^{i(p_1x_{k+1} + p_2y_{k+1})} - e^{i(p_1x_{k} + p_2y_{k})}\big)
\frac{p_1(y_{k+1}-y_k) - p_2(x_{k+1}-x_k)}{p_1(x_{k+1}-x_k) + p_2(y_{k+1}-y_k)}
$$
In general, the integral of $e^{ip\cdot r}$ over polygon $\Omega$
will be a linear combination of the value of $e^{ip\cdot r}$ at the vertices $v_k$ weighted by rational function of coordinates of $p$ and $v_k$.
The formula $(*1)$ is valid only when all $p\cdot v_k \ne p \cdot v_{k+1}$. In the case when this fails, one could use the fact the integral depends on $p$ continuously and extract the desired coefficient in front of $e^{ip\cdot v_k}$ by taking appropriate limit in $(*1)$. Introduce two set of new variables
$$r_k = \frac12(v_{k+1} + v_{k})\quad\text{ and }\quad \epsilon_k = \frac12(v_{k+1} - v_k)$$
We obtain following expression which works as long as $p \ne 0$.
$$\int_\Omega e^{ip\cdot r} dxdy
= \frac{2}{i|p|^2}\sum_{k=1}^n e^{ip\cdot r_k} {\rm sinc}(p\cdot \epsilon_k) 
 ( \hat{z} \cdot p \times \epsilon_k ),\quad
{\rm sinc}(x) = \begin{cases}
\frac{\sin x}{x}, & x \ne 0\\
1, & x = 0
\end{cases}
$$
I worked out above formulas but this stuff has been known for a long time.
A literature search return a recent article
Form factor (Fourier shape transform) of polygon and polyhedron by Joachim Wuttke. Look at that and
references there for more details.
